ENH: Update estimate_utm_crs to handle crossing the antimeridian with pyproj 3.1+ (#2049)

ci/envs/37-pd10.yaml

@@ -23,7 +23,7 @@ dependencies:
   - libspatialite
   - pip
   - pip:
-      - pyproj==2.3.1
+      - pyproj==3.0.1
       - geopy
       - mapclassify==2.4.0
       - pyarrow

geopandas/_compat.py

@@ -223,3 +223,4 @@ except ImportError:
 # -----------------------------------------------------------------------------
 
 PYPROJ_LT_3 = LooseVersion(pyproj.__version__) < LooseVersion("3")
+PYPROJ_GE_31 = LooseVersion(pyproj.__version__) >= LooseVersion("3.1")

geopandas/array.py

@@ -828,16 +828,34 @@ class GeometryArray(ExtensionArray):
             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:
+        if self.crs.is_geographic:
             x_center = np.mean([minx, maxx])
             y_center = np.mean([miny, maxy])
+        # ensure using geographic coordinates
+        else:
+            transformer = Transformer.from_crs(self.crs, "EPSG:4326", always_xy=True)
+            if compat.PYPROJ_GE_31:
+                minx, miny, maxx, maxy = transformer.transform_bounds(
+                    minx, miny, maxx, maxy
+                )
+                y_center = np.mean([miny, maxy])
+                # crossed the antimeridian
+                if minx > maxx:
+                    # shift maxx from [-180,180] to [0,360]
+                    # so both numbers are positive for center calculation
+                    # Example: -175 to 185
+                    maxx += 360
+                    x_center = np.mean([minx, maxx])
+                    # shift back to [-180,180]
+                    x_center = ((x_center + 180) % 360) - 180
+                else:
+                    x_center = np.mean([minx, maxx])
+            else:
+                lon, lat = transformer.transform(
+                    (minx, maxx, minx, maxx), (miny, miny, maxy, maxy)
+                )
+                x_center = np.mean(lon)
+                y_center = np.mean(lat)
 
         utm_crs_list = query_utm_crs_info(
             datum_name=datum_name,

geopandas/tests/test_array.py

@@ -948,6 +948,17 @@ class TestEstimateUtmCrs:
             "EPSG:32618"
         )
 
+    @pytest.mark.skipif(not compat.PYPROJ_GE_31, reason="requires pyproj 3.1 or higher")
+    def test_estimate_utm_crs__antimeridian(self):
+        antimeridian = from_shapely(
+            [
+                shapely.geometry.Point(1722483.900174921, 5228058.6143420935),
+                shapely.geometry.Point(4624385.494808555, 8692574.544944234),
+            ],
+            crs="EPSG:3851",
+        )
+        assert antimeridian.estimate_utm_crs() == CRS("EPSG:32760")
+
     @pytest.mark.skipif(compat.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"):