From c90ce139124d2012155fe14ab6d33f0da8a6b2aa Mon Sep 17 00:00:00 2001 From: "Alan D. Snow" Date: Fri, 24 Sep 2021 13:47:52 -0500 Subject: [PATCH] ENH: Update estimate_utm_crs to handle crossing the antimeridian with pyproj 3.1+ (#2049) --- ci/envs/37-pd10.yaml | 2 +- geopandas/_compat.py | 1 + geopandas/array.py | 34 ++++++++++++++++++++++++++-------- geopandas/tests/test_array.py | 11 +++++++++++ 4 files changed, 39 insertions(+), 9 deletions(-) diff --git a/ci/envs/37-pd10.yaml b/ci/envs/37-pd10.yaml index 3d1d5ff..a42aee3 100644 --- a/ci/envs/37-pd10.yaml +++ b/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 diff --git a/geopandas/_compat.py b/geopandas/_compat.py index a93139a..dbe5311 100644 --- a/geopandas/_compat.py +++ b/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") diff --git a/geopandas/array.py b/geopandas/array.py index bf09528..f4924c8 100644 --- a/geopandas/array.py +++ b/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, diff --git a/geopandas/tests/test_array.py b/geopandas/tests/test_array.py index 2ddc667..885f92d 100644 --- a/geopandas/tests/test_array.py +++ b/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"):