From a0edbee7d2d950db6d4f9d8068f0b1e09004e51d Mon Sep 17 00:00:00 2001 From: "Alan D. Snow" Date: Fri, 27 Nov 2020 04:15:26 -0600 Subject: [PATCH] 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 * pytest.mark.skipif Co-authored-by: Martin Fleischmann * Update geopandas/geodataframe.py Co-authored-by: Martin Fleischmann * more skipif * use pyproj to transform bounds Co-authored-by: Martin Fleischmann --- doc/source/docs/reference/geodataframe.rst | 1 + doc/source/docs/reference/geoseries.rst | 1 + geopandas/geodataframe.py | 40 ++++++++++++ geopandas/geoseries.py | 73 ++++++++++++++++++++++ geopandas/tests/test_geodataframe.py | 14 +++++ geopandas/tests/test_geoseries.py | 32 ++++++++++ 6 files changed, 161 insertions(+) diff --git a/doc/source/docs/reference/geodataframe.rst b/doc/source/docs/reference/geodataframe.rst index 05ae651..71c1a85 100644 --- a/doc/source/docs/reference/geodataframe.rst +++ b/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 ------------------------ diff --git a/doc/source/docs/reference/geoseries.rst b/doc/source/docs/reference/geoseries.rst index fb545b2..a1499d3 100644 --- a/doc/source/docs/reference/geoseries.rst +++ b/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 -------------- diff --git a/geopandas/geodataframe.py b/geopandas/geodataframe.py index 0257919..e9a2047 100644 --- a/geopandas/geodataframe.py +++ b/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 + + 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 diff --git a/geopandas/geoseries.py b/geopandas/geoseries.py index b12ae0a..02e5f11 100644 --- a/geopandas/geoseries.py +++ b/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 + + 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. diff --git a/geopandas/tests/test_geodataframe.py b/geopandas/tests/test_geodataframe.py index 07d011e..4c6bc69 100644 --- a/geopandas/tests/test_geodataframe.py +++ b/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) diff --git a/geopandas/tests/test_geoseries.py b/geopandas/tests/test_geoseries.py index 72b6f23..9bcc763 100644 --- a/geopandas/tests/test_geoseries.py +++ b/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()