ENH: Add estimate_utm_crs method to GeoSeries and GeoDataFrame (#1646)

* ENH: Add estimate_tum_crs method to GeoSeries and GeoDataFrame

* Add examples & expect release in 0.9

* undo changes in CI

* Update geopandas/geoseries.py

Co-authored-by: Martin Fleischmann <martin@martinfleischmann.net>

* pytest.mark.skipif

Co-authored-by: Martin Fleischmann <martin@martinfleischmann.net>

* Update geopandas/geodataframe.py

Co-authored-by: Martin Fleischmann <martin@martinfleischmann.net>

* more skipif

* use pyproj to transform bounds

Co-authored-by: Martin Fleischmann <martin@martinfleischmann.net>

doc/source/docs/reference/geodataframe.rst

@@ -37,6 +37,7 @@ Projection handling
    GeoDataFrame.crs
    GeoDataFrame.set_crs
    GeoDataFrame.to_crs
+   GeoDataFrame.estimate_utm_crs
 
 Active geometry handling
 ------------------------

doc/source/docs/reference/geoseries.rst

@@ -123,6 +123,7 @@ Projection handling
    GeoSeries.crs
    GeoSeries.set_crs
    GeoSeries.to_crs
+   GeoSeries.estimate_utm_crs
 
 Missing values
 --------------

geopandas/geodataframe.py

@@ -1075,6 +1075,46 @@ box': (2.0, 1.0, 2.0, 1.0)}], 'bbox': (1.0, 1.0, 2.0, 2.0)}
         if not inplace:
             return df
 
+    def estimate_utm_crs(self, datum_name="WGS 84"):
+        """Returns the estimated UTM CRS based on the bounds of the dataset.
+
+        .. versionadded:: 0.9
+
+        .. note:: Requires pyproj 3+
+
+        Parameters
+        ----------
+        datum_name : str, optional
+            The name of the datum to use in the query. Default is WGS 84.
+
+        Returns
+        -------
+        pyproj.CRS
+
+        Examples
+        --------
+        >>> world = geopandas.read_file(
+        ...     geopandas.datasets.get_path("naturalearth_lowres")
+        ... )
+        >>> germany = world.loc[world.name == "Germany"]
+        >>> germany.estimate_utm_crs()  # doctest: +SKIP
+        <Projected CRS: EPSG:32632>
+        Name: WGS 84 / UTM zone 32N
+        Axis Info [cartesian]:
+        - E[east]: Easting (metre)
+        - N[north]: Northing (metre)
+        Area of Use:
+        - name: World - N hemisphere - 6°E to 12°E - by country
+        - bounds: (6.0, 0.0, 12.0, 84.0)
+        Coordinate Operation:
+        - name: UTM zone 32N
+        - method: Transverse Mercator
+        Datum: World Geodetic System 1984
+        - Ellipsoid: WGS 84
+        - Prime Meridian: Greenwich
+        """
+        return self.geometry.estimate_utm_crs(datum_name=datum_name)
+
     def __getitem__(self, key):
         """
         If the result is a column containing only 'geometry', return a

geopandas/geoseries.py

@@ -777,6 +777,79 @@ class GeoSeries(GeoPandasBase, Series):
             GeometryArray(new_data), crs=crs, index=self.index, name=self.name
         )
 
+    def estimate_utm_crs(self, datum_name="WGS 84"):
+        """Returns the estimated UTM CRS based on the bounds of the dataset.
+
+        .. versionadded:: 0.9
+
+        .. note:: Requires pyproj 3+
+
+        Parameters
+        ----------
+        datum_name : str, optional
+            The name of the datum to use in the query. Default is WGS 84.
+
+        Returns
+        -------
+        pyproj.CRS
+
+        Examples
+        --------
+        >>> world = geopandas.read_file(
+        ...     geopandas.datasets.get_path("naturalearth_lowres")
+        ... )
+        >>> germany = world.loc[world.name == "Germany"]
+        >>> germany.geometry.estimate_utm_crs()  # doctest: +SKIP
+        <Projected CRS: EPSG:32632>
+        Name: WGS 84 / UTM zone 32N
+        Axis Info [cartesian]:
+        - E[east]: Easting (metre)
+        - N[north]: Northing (metre)
+        Area of Use:
+        - name: World - N hemisphere - 6°E to 12°E - by country
+        - bounds: (6.0, 0.0, 12.0, 84.0)
+        Coordinate Operation:
+        - name: UTM zone 32N
+        - method: Transverse Mercator
+        Datum: World Geodetic System 1984
+        - Ellipsoid: WGS 84
+        - Prime Meridian: Greenwich
+        """
+        try:
+            from pyproj.aoi import AreaOfInterest
+            from pyproj.database import query_utm_crs_info
+        except ImportError:
+            raise RuntimeError("pyproj 3+ required for estimate_utm_crs.")
+
+        if not self.crs:
+            raise RuntimeError("crs must be set to estimate UTM CRS.")
+
+        minx, miny, maxx, maxy = self.total_bounds
+        # ensure using geographic coordinates
+        if not self.crs.is_geographic:
+            lon, lat = Transformer.from_crs(
+                self.crs, "EPSG:4326", always_xy=True
+            ).transform((minx, maxx, minx, maxx), (miny, miny, maxy, maxy))
+            x_center = np.mean(lon)
+            y_center = np.mean(lat)
+        else:
+            x_center = np.mean([minx, maxx])
+            y_center = np.mean([miny, maxy])
+
+        utm_crs_list = query_utm_crs_info(
+            datum_name=datum_name,
+            area_of_interest=AreaOfInterest(
+                west_lon_degree=x_center,
+                south_lat_degree=y_center,
+                east_lon_degree=x_center,
+                north_lat_degree=y_center,
+            ),
+        )
+        try:
+            return CRS.from_epsg(utm_crs_list[0].code)
+        except IndexError:
+            raise RuntimeError("Unable to determine UTM CRS")
+
     def to_json(self, **kwargs):
         """
         Returns a GeoJSON string representation of the GeoSeries.

geopandas/tests/test_geodataframe.py

@@ -2,11 +2,14 @@ import json
 import os
 import shutil
 import tempfile
+from distutils.version import LooseVersion
 
 import numpy as np
 import pandas as pd
 
 import fiona
+import pyproj
+from pyproj import CRS
 from pyproj.exceptions import CRSError
 from shapely.geometry import Point
 
@@ -20,6 +23,9 @@ from pandas.testing import assert_frame_equal, assert_index_equal, assert_series
 import pytest
 
 
+PYPROJ_LT_3 = LooseVersion(pyproj.__version__) < LooseVersion("3")
+
+
 class TestDataFrame:
     def setup_method(self):
         N = 10
@@ -659,6 +665,14 @@ class TestDataFrame:
         assert_frame_equal(self.df, unpickled)
         assert self.df.crs == unpickled.crs
 
+    def test_estimate_utm_crs(self):
+        if PYPROJ_LT_3:
+            with pytest.raises(RuntimeError, match=r"pyproj 3\+ required"):
+                self.df.estimate_utm_crs()
+        else:
+            assert self.df.estimate_utm_crs() == CRS("EPSG:32618")
+            assert self.df.estimate_utm_crs("NAD83") == CRS("EPSG:26918")
+
 
 def check_geodataframe(df, geometry_column="geometry"):
     assert isinstance(df, GeoDataFrame)

geopandas/tests/test_geoseries.py

@@ -1,3 +1,4 @@
+from distutils.version import LooseVersion
 import json
 import os
 import random
@@ -8,6 +9,7 @@ import numpy as np
 from numpy.testing import assert_array_equal
 import pandas as pd
 
+from pyproj import CRS
 from shapely.geometry import (
     LineString,
     MultiLineString,
@@ -17,6 +19,7 @@ from shapely.geometry import (
     Polygon,
 )
 from shapely.geometry.base import BaseGeometry
+import pyproj
 
 from geopandas import GeoSeries, GeoDataFrame
 from geopandas.array import GeometryArray, GeometryDtype
@@ -26,6 +29,9 @@ from pandas.testing import assert_series_equal
 import pytest
 
 
+PYPROJ_LT_3 = LooseVersion(pyproj.__version__) < LooseVersion("3")
+
+
 class TestSeries:
     def setup_method(self):
         self.tempdir = tempfile.mkdtemp()
@@ -182,6 +188,32 @@ class TestSeries:
         with pytest.raises(ValueError):
             self.landmarks.to_crs(crs=None, epsg=None)
 
+    def test_estimate_utm_crs__geographic(self):
+        if PYPROJ_LT_3:
+            with pytest.raises(RuntimeError, match=r"pyproj 3\+ required"):
+                self.landmarks.estimate_utm_crs()
+        else:
+            assert self.landmarks.estimate_utm_crs() == CRS("EPSG:32618")
+            assert self.landmarks.estimate_utm_crs("NAD83") == CRS("EPSG:26918")
+
+    @pytest.mark.skipif(PYPROJ_LT_3, reason="requires pyproj 3 or higher")
+    def test_estimate_utm_crs__projected(self):
+        assert self.landmarks.to_crs("EPSG:3857").estimate_utm_crs() == CRS(
+            "EPSG:32618"
+        )
+
+    @pytest.mark.skipif(PYPROJ_LT_3, reason="requires pyproj 3 or higher")
+    def test_estimate_utm_crs__out_of_bounds(self):
+        with pytest.raises(RuntimeError, match="Unable to determine UTM CRS"):
+            GeoSeries(
+                [Polygon([(0, 90), (1, 90), (2, 90)])], crs="EPSG:4326"
+            ).estimate_utm_crs()
+
+    @pytest.mark.skipif(PYPROJ_LT_3, reason="requires pyproj 3 or higher")
+    def test_estimate_utm_crs__missing_crs(self):
+        with pytest.raises(RuntimeError, match="crs must be set"):
+            GeoSeries([Polygon([(0, 90), (1, 90), (2, 90)])]).estimate_utm_crs()
+
     def test_fillna(self):
         # default is to fill with empty geometry
         na = self.na_none.fillna()