geopandas/doc/source/docs/user_guide/projections.rst

.. currentmodule:: geopandas

.. ipython:: python
   :suppress:

   import geopandas


Managing Projections
=========================================


Coordinate Reference Systems
-----------------------------

The Coordinate Reference System (CRS) is important because the geometric shapes
in a GeoSeries or GeoDataFrame object are simply a collection of coordinates in
an arbitrary space. A CRS tells Python how those coordinates relate to places on
the Earth.

You can find the codes for most commonly used projections from
`www.spatialreference.org <https://spatialreference.org/>`_.

The same CRS can often be referred to in many ways. For example, one of the most
commonly used CRS is the WGS84 latitude-longitude projection. This can be
referred to using the authority code ``"EPSG:4326"``.

*geopandas* can accept anything accepted by :meth:`pyproj.CRS.from_user_input() <pyproj.crs.CRS.from_user_input>`:

- CRS WKT string
- An authority string (i.e. "epsg:4326")
- An EPSG integer code (i.e. 4326)
- A ``pyproj.CRS``
- An object with a to_wkt method.
- PROJ string
- Dictionary of PROJ parameters
- PROJ keyword arguments for parameters
- JSON string with PROJ parameters

For reference, a few very common projections and their EPSG codes:

* WGS84 Latitude/Longitude: ``"EPSG:4326"``
* UTM Zones (North): ``"EPSG:32633"``
* UTM Zones (South): ``"EPSG:32733"``


What is the best format to store the CRS information?
-----------------------------------------------------

Generally, WKT or SRID's are preferred over PROJ strings as they can contain more information about a given CRS.
Conversions between WKT and PROJ strings will in most cases cause a loss of information, potentially leading to erroneous transformations. If possible WKT2 should be used.

For more details, see https://proj.org/faq.html#what-is-the-best-format-for-describing-coordinate-reference-systems


Setting a Projection
----------------------

There are two relevant operations for projections: setting a projection and re-projecting.

Setting a projection may be necessary when for some reason *geopandas* has coordinate data (x-y values), but no information about how those coordinates refer to locations in the real world. Setting a projection is how one tells *geopandas* how to interpret coordinates. If no CRS is set, *geopandas* geometry operations will still work, but coordinate transformations will not be possible and exported files may not be interpreted correctly by other software.

Be aware that **most of the time** you don't have to set a projection. Data loaded from a reputable source (using the :func:`geopandas.read_file()` command) *should* always include projection information. You can see an objects current CRS through the :attr:`GeoSeries.crs` attribute.

From time to time, however, you may get data that does not include a projection. In this situation, you have to set the CRS so *geopandas* knows how to interpret the coordinates.

For example, if you convert a spreadsheet of latitudes and longitudes into a
GeoSeries by hand, you would set the projection by passing the WGS84
latitude-longitude CRS to the :meth:`GeoSeries.set_crs` method (or by setting
the :attr:`GeoSeries.crs` attribute):

.. sourcecode:: python

    my_geoseries = my_geoseries.set_crs("EPSG:4326")
    my_geoseries = my_geoseries.set_crs(epsg=4326)


Re-Projecting
----------------

Re-projecting is the process of changing the representation of locations from one coordinate system to another. All projections of locations on the Earth into a two-dimensional plane `are distortions <https://en.wikipedia.org/wiki/Map_projection#Which_projection_is_best.3F>`_, the projection that is best for your application may be different from the projection associated with the data you import. In these cases, data can be re-projected using the :meth:`GeoDataFrame.to_crs` command:

.. ipython:: python

    # load example data
    world = geopandas.read_file(geopandas.datasets.get_path('naturalearth_lowres'))

    # Check original projection
    # (it's Platte Carre! x-y are long and lat)
    world.crs

    # Visualize
    ax = world.plot()
    @savefig world_starting.png
    ax.set_title("WGS84 (lat/lon)");

    # Reproject to Mercator (after dropping Antartica)
    world = world[(world.name != "Antarctica") & (world.name != "Fr. S. Antarctic Lands")]
    world = world.to_crs("EPSG:3395") # world.to_crs(epsg=3395) would also work
    ax = world.plot()
    @savefig world_reproj.png
    ax.set_title("Mercator");


Projection for multiple geometry columns
----------------------------------------

GeoPandas 0.8 implements support for different projections assigned to different geometry
columns of the same GeoDataFrame. The projection is now stored together with geometries per column (directly
on the GeometryArray level).

Note that if GeometryArray has assigned projection, it is preferred over the
projection passed to GeoSeries or GeoDataFrame during the creation:

.. code-block:: python

   >>> array.crs
   <Geographic 2D CRS: EPSG:4326>
   Name: WGS 84
   Axis Info [ellipsoidal]:
   - Lat[north]: Geodetic latitude (degree)
   - Lon[east]: Geodetic longitude (degree)
   ...
   >>> GeoSeries(array, crs=3395).crs  # crs=3395 is ignored as array already has CRS
   FutureWarning: CRS mismatch between CRS of the passed geometries and 'crs'. Use 'GeoDataFrame.set_crs(crs, allow_override=True)' to overwrite CRS or 'GeoDataFrame.to_crs(crs)' to reproject geometries. CRS mismatch will raise an error in the future versions of GeoPandas.
       GeoSeries(array, crs=3395).crs

   <Geographic 2D CRS: EPSG:4326>
   Name: WGS 84
   Axis Info [ellipsoidal]:
   - Lat[north]: Geodetic latitude (degree)
   - Lon[east]: Geodetic longitude (degree)
   ...

If you want to overwrite projection, you can then assign it to the GeoSeries
manually or re-project geometries to the target projection using either
``GeoSeries.set_crs(epsg=3395, allow_override=True)`` or
``GeoSeries.to_crs(epsg=3395)``.

All GeometryArray-based operations preserve projection; however, if you loop over a column
containing geometry, this information might be lost.


Upgrading to GeoPandas 0.7 with pyproj > 2.2 and PROJ > 6
---------------------------------------------------------

Starting with GeoPandas 0.7, the `.crs` attribute of a GeoSeries or GeoDataFrame
stores the CRS information as a ``pyproj.CRS``, and no longer as a proj4 string
or dict.

Before, you might have seen this:

.. code-block:: python

   >>> gdf.crs
   {'init': 'epsg:4326'}

while now you will see something like this:

.. code-block:: python

   >>> gdf.crs
   <Geographic 2D CRS: EPSG:4326>
   Name: WGS 84
   Axis Info [ellipsoidal]:
   - Lat[north]: Geodetic latitude (degree)
   - Lon[east]: Geodetic longitude (degree)
   ...
   >>> type(gdf.crs)
   pyproj.crs.CRS

This gives a better user interface and integrates improvements from pyproj and
PROJ 6, but might also require some changes in your code. See `this blogpost
<https://jorisvandenbossche.github.io/blog/2020/02/11/geopandas-pyproj-crs/>`__
for some more background, and the subsections below cover different possible
migration issues.

See the `pyproj docs <https://pyproj4.github.io/pyproj/stable/>`__ for more on
the ``pyproj.CRS`` object.

Importing data from files
^^^^^^^^^^^^^^^^^^^^^^^^^

When reading geospatial files with :func:`geopandas.read_file`, things should
mostly work out of the box. For example, reading the example countries dataset
yields a proper CRS:

.. ipython:: python

   df = geopandas.read_file(geopandas.datasets.get_path('naturalearth_lowres'))
   df.crs

However, in certain cases (with older CRS formats), the resulting CRS object
might not be fully as expected. See the :ref:`section below <unrecognized-crs-reasons>`
for possible reasons and how to solve it.


Manually specifying the CRS
^^^^^^^^^^^^^^^^^^^^^^^^^^^

When specifying the CRS manually in your code (e.g., because your data has not
yet a CRS, or when converting to another CRS), this might require a change in
your code.

**"init" proj4 strings/dicts**

Currently, a lot of people (and also the GeoPandas docs showed that before)
specify the EPSG code using the "init" proj4 string:

.. code-block:: python

   ## OLD
   GeoDataFrame(..., crs={'init': 'epsg:4326'})
   # or
   gdf.crs = {'init': 'epsg:4326'}
   # or
   gdf.to_crs({'init': 'epsg:4326'})

The above will now raise a deprecation warning from pyproj, and instead of the
"init" proj4 string, you should use only the EPSG code itself as follows:

.. code-block:: python

   ## NEW
   GeoDataFrame(..., crs="EPSG:4326")
   # or
   gdf.crs = "EPSG:4326"
   # or
   gdf.to_crs("EPSG:4326")


**proj4 strings/dicts**

Although a full proj4 string is not deprecated (as opposed to the "init" string
above), it is still recommended to change it with an EPSG code if possible.

For example, instead of:

.. code-block:: python

   gdf.crs = "+proj=laea +lat_0=45 +lon_0=-100 +x_0=0 +y_0=0 +a=6370997 +b=6370997 +units=m +no_defs"

we recommend to do:

.. code-block:: python

   gdf.crs = "EPSG:2163"

*if* you know the EPSG code for the projection you are using.

One possible way to find out the EPSG code is using pyproj for this:

.. code-block:: python

   >>> import pyproj
   >>> crs = pyproj.CRS("+proj=laea +lat_0=45 +lon_0=-100 +x_0=0 +y_0=0 +a=6370997 +b=6370997 +units=m +no_defs")
   >>> crs.to_epsg()
   2163

(you might need to set the ``min_confidence`` keyword of ``to_epsg`` to a lower
value if the match is not perfect)

Further, on websites such as `spatialreference.org <https://spatialreference.org/>`__
and `epsg.io <https://epsg.io/>`__ the descriptions of many CRS can be found
including their EPSG codes and proj4 string definitions.

**Other formats**

Next to the EPSG code mentioned above, there are also other ways to specify the
CRS: an actual ``pyproj.CRS`` object, a WKT string, a PROJ JSON string, etc.
Anything that is accepted by ``pyproj.CRS.from_user_input`` can by specified
to the ``crs`` keyword/attribute in GeoPandas.

Also compatible CRS objects, such as from the ``rasterio`` package, can be
passed directly to GeoPandas.


The axis order of a CRS
^^^^^^^^^^^^^^^^^^^^^^^

Starting with PROJ 6 / pyproj 2, the axis order of the official EPSG definition
is honoured. For example, when using geographic coordinates (degrees of longitude
and latitude) in the standard EPSG:4326, the CRS will look like:

.. code-block:: python

   >>> pyproj.CRS(3EPSG:4326")
   <Geographic 2D CRS: EPSG:4326>
   ...
   Axis Info [ellipsoidal]:
   - Lat[north]: Geodetic latitude (degree)
   - Lon[east]: Geodetic longitude (degree)
   ...

This mentions the order as (lat, lon), as that is the official order of coordinates
in EPSG:4326. In GeoPandas, however, the coordinates are always stored as (x, y),
and thus as (lon, lat) order, regardless of the CRS (i.e. the "traditional" order used
in GIS). When reprojecting, GeoPandas and pyproj will under the hood take care of
this difference in axis order, so the user doesn't need to care about this.

.. _unrecognized-crs-reasons:

Why is it not properly recognizing my CRS?
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

There are many file sources and CRS definitions out there "in the wild" that
might have a CRS description that does not fully conform to the new standards of
PROJ > 6 (proj4 strings, older WKT formats, ...). In such cases, you will get a
``pyproj.CRS`` object that might not be fully what you expected (e.g. not equal
to the expected EPSG code). Below we list a few possible cases.

I get a "Bound CRS"?
~~~~~~~~~~~~~~~~~~~~

Some CRS definitions include a *"towgs84" clause*, which can give problems in
recognizing the actual CRS.

For example, both the proj4 and WKT representation for EPSG:31370 (the local
projection used in Belgium) as can be found at `https://spatialreference.org/ref/epsg/31370/ <https://spatialreference.org/ref/epsg/31370/>`__
include this. When taking one of those definitions from that site, and creating
a CRS object:

.. code-block:: python

   >>> import pyproj
   >>> crs = pyproj.CRS("+proj=lcc +lat_1=51.16666723333333 +lat_2=49.8333339 +lat_0=90 +lon_0=4.367486666666666 +x_0=150000.013 +y_0=5400088.438 +ellps=intl +towgs84=106.869,-52.2978,103.724,-0.33657,0.456955,-1.84218,1 +units=m +no_defs")
   >>> crs
   <Bound CRS: +proj=lcc +lat_1=51.16666723333333 +lat_2=49.83333 ...>
   Name: unknown
   Axis Info [cartesian]:
   - E[east]: Easting (metre)
   - N[north]: Northing (metre)
   Area of Use:
   - undefined
   Coordinate Operation:
   - name: Transformation from unknown to WGS84
   - method: Position Vector transformation (geog2D domain)
   Datum: Unknown based on International 1909 (Hayford) ellipsoid
   - Ellipsoid: International 1909 (Hayford)
   - Prime Meridian: Greenwich
   Source CRS: unknown

You notice that the above is a not a "Projected CRS" as expected, but a "Bound CRS".
This is because it is "bound" to a conversion to WGS84, and will always use this
when reprojecting instead of letting PROJ determine the best conversion.

To get the actual underlying projected CRS, you can use the ``.source_crs`` attribute:

.. code-block:: python

   >>> crs.source_crs
   <Projected CRS: PROJCRS["unknown",BASEGEOGCRS["unknown",DATUM["Unk ...>
   Name: unknown
   ...

Now we have a "Projected CRS", and now it will also recognize the correct EPSG
number:

.. code-block:: python

   >>> crs.to_epsg()

   >>> crs.source_crs.to_epsg()
   31370

I have a different axis order?
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

As mentioned above, pyproj now honours the axis order of the EPSG definition.
However, proj4 strings or older WKT versions don't specify this correctly, which
can be a reason that the CRS object is not equal to the expected EPSG code.

Consider the following example of a Canadian projected CRS "EPSG:2953". When
constructing the CRS object from the WKT string as provided on
`https://epsg.io/2953 <https://epsg.io/2953>`__:

.. code-block:: python

   >>> crs = pyproj.CRS("""PROJCS["NAD83(CSRS) / New Brunswick Stereographic",
   ...     GEOGCS["NAD83(CSRS)",
   ...         DATUM["NAD83_Canadian_Spatial_Reference_System",
   ...             SPHEROID["GRS 1980",6378137,298.257222101,
   ...                 AUTHORITY["EPSG","7019"]],
   ...             AUTHORITY["EPSG","6140"]],
   ...         PRIMEM["Greenwich",0,
   ...             AUTHORITY["EPSG","8901"]],
   ...         UNIT["degree",0.0174532925199433,
   ...             AUTHORITY["EPSG","9122"]],
   ...         AUTHORITY["EPSG","4617"]],
   ...     PROJECTION["Oblique_Stereographic"],
   ...     PARAMETER["latitude_of_origin",46.5],
   ...     PARAMETER["central_meridian",-66.5],
   ...     PARAMETER["scale_factor",0.999912],
   ...     PARAMETER["false_easting",2500000],
   ...     PARAMETER["false_northing",7500000],
   ...     UNIT["metre",1,
   ...         AUTHORITY["EPSG","9001"]],
   ...     AUTHORITY["EPSG","2953"]]""")

   >>> crs
   <Projected CRS: PROJCS["NAD83(CSRS) / New Brunswick Stereographic" ...>
   Name: NAD83(CSRS) / New Brunswick Stereographic
   Axis Info [cartesian]:
   - E[east]: Easting (metre)
   - N[north]: Northing (metre)
   ...

Although this is the WKT string as found online for "EPSG:2953", this CRS object
does not evaluate equal to this EPSG code:

.. code-block:: python

   >>> crs == "EPSG:2953"
   False

If we construct the CRS object from the EPSG code (truncated output):

.. code-block:: python

   >>> pyproj.CRS("EPSG:2953")
   <Projected CRS: EPSG:2953>
   Name: NAD83(CSRS) / New Brunswick Stereographic
   Axis Info [cartesian]:
   - N[north]: Northing (metre)
   - E[east]: Easting (metre)
   ...

You can see that the CRS object constructed from the WKT string has a "Easting,
Northing" (i.e. x, y) axis order, while the CRS object constructed from the EPSG
code has a (Northing, Easting) axis order.

Only having this difference in axis order is no problem when using the CRS in
GeoPandas, since GeoPandas always uses a (x, y) order to store the data
regardless of the CRS definition. But, you might still want to verify it is
equivalent to the expected EPSG code. By lowering the `min_confidence`, the axis
order will be ignored:

.. code-block:: python

   >>> crs.to_epsg()

   >>> crs.to_epsg(min_confidence=20)
   2953


The ``.crs`` attribute is no longer a dict or string
^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^

If you relied on the ``.crs`` object being a dict or a string, such code can
be broken given it is now a ``pyproj.CRS`` object. But this object actually
provides a more robust interface to get information about the CRS.

For example, if you used the following code to get the EPSG code:

.. code-block:: python

   gdf.crs['init']

This will no longer work. To get the EPSG code from a ``crs`` object, you can use
the ``to_epsg()`` method.

Or to check if a CRS was a certain UTM zone:

.. code-block:: python

   '+proj=utm ' in gdf.crs

could be replaced with the more robust check (requires pyproj 2.6+):

.. code-block:: python

   gdf.crs.utm_zone is not None

And there are many other methods available on the ``pyproj.CRS`` class to get
information about the CRS.