ENH: Added rio.estimate_utm_crs

.travis.yml

@@ -54,7 +54,7 @@ before_install:
 - conda config --add channels conda-forge
 - conda config --set channel_priority strict
 # Create conda environment
-- conda create -n test python=$PYTHON_VERSION rasterio scipy xarray netcdf4 dask pandoc
+- conda create -n test python=$PYTHON_VERSION rasterio xarray scipy pyproj netcdf4 dask pandoc
 - source $HOME/miniconda/bin/activate test
 
 install:

docs/getting_started/crs_management.ipynb

@@ -35,7 +35,8 @@
     "#### API Documentation\n",
     "\n",
     "- [rio.write_crs()](../rioxarray.rst#rioxarray.rioxarray.XRasterBase.write_crs)\n",
-    "- [rio.crs](../rioxarray.rst#rioxarray.rioxarray.XRasterBase.crs)\n"
+    "- [rio.crs](../rioxarray.rst#rioxarray.rioxarray.XRasterBase.crs)\n",
+    "- [rio.estimate_utm_crs()](../rioxarray.rst#rioxarray.rioxarray.XRasterBase.estimate_utm_crs)"
    ]
   },
   {
@@ -327,7 +328,7 @@
    "name": "python",
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
-   "version": "3.8.3"
+   "version": "3.8.6"
   }
  },
  "nbformat": 4,

docs/history.rst

@@ -3,7 +3,7 @@ History
 
 Latest
 ------
-
+- ENH: Added `rio.estimate_utm_crs` (issue #181)
 
 0.1.1
 ------

rioxarray/rioxarray.py

@@ -496,6 +496,54 @@ def write_crs(self, input_crs=None, grid_mapping_name=None, inplace=False):
             grid_mapping_name=grid_mapping_name, inplace=True
         )
 
+    def estimate_utm_crs(self, datum_name="WGS 84"):
+        """Returns the estimated UTM CRS based on the bounds of the dataset.
+
+        .. versionadded:: 0.2
+
+        .. note:: Requires pyproj 3+
+
+        Parameters
+        ----------
+        datum_name : str, optional
+            The name of the datum to use in the query. Default is WGS 84.
+
+        Returns
+        -------
+        rasterio.crs.CRS
+        """
+        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 self.crs is None:
+            raise RuntimeError("crs must be set to estimate UTM CRS.")
+
+        # ensure using geographic coordinates
+        if self.crs.is_geographic:
+            minx, miny, maxx, maxy = self.bounds(recalc=True)
+        else:
+            minx, miny, maxx, maxy = self.transform_bounds("EPSG:4326", recalc=True)
+
+        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 _cached_transform(self):
         """
         Get the transform from:

test/integration/test_integration_rioxarray.py

@@ -5,6 +5,7 @@
 from functools import partial
 
 import numpy
+import pyproj
 import pytest
 import rasterio
 import scipy
@@ -31,6 +32,8 @@
     _assert_xarrays_equal,
 )
 
+PYPROJ_LT_3 = LooseVersion(pyproj.__version__) < LooseVersion("3")
+
 
 @pytest.fixture(params=[xarray.open_dataset, xarray.open_dataarray])
 def modis_reproject(request):
@@ -2164,3 +2167,42 @@ def test_add_spatial_ref_warning():
 def test_add_xy_grid_meta_warning():
     with pytest.raises(RuntimeError):
         add_xy_grid_meta({})
+
+
+def test_estimate_utm_crs():
+    xds = rioxarray.open_rasterio(
+        os.path.join(TEST_INPUT_DATA_DIR, "cog.tif"),
+    )
+    if PYPROJ_LT_3:
+        with pytest.raises(RuntimeError, match=r"pyproj 3\+ required"):
+            xds.rio.estimate_utm_crs()
+    else:
+        assert xds.rio.estimate_utm_crs() == CRS.from_epsg(32618)
+        assert xds.rio.reproject("EPSG:4326").rio.estimate_utm_crs() == CRS.from_epsg(
+            32618
+        )
+        assert xds.rio.estimate_utm_crs("NAD83") == CRS.from_epsg(26918)
+
+
+@pytest.mark.skipif(PYPROJ_LT_3, reason="pyproj 3+ required")
+def test_estimate_utm_crs__missing_crs():
+    with pytest.raises(RuntimeError, match=r"crs must be set to estimate UTM CRS"):
+        xarray.Dataset().rio.estimate_utm_crs("NAD83")
+
+
+def test_estimate_utm_crs__out_of_bounds():
+    xds = xarray.DataArray(
+        numpy.zeros((2, 2)),
+        dims=("latitude", "longitude"),
+        coords={
+            "latitude": [-90.0, -90.0],
+            "longitude": [-5.0, 5.0],
+        },
+    )
+    xds.rio.write_crs("EPSG:4326", inplace=True)
+    if PYPROJ_LT_3:
+        with pytest.raises(RuntimeError, match=r"pyproj 3\+ required"):
+            xds.rio.estimate_utm_crs()
+    else:
+        with pytest.raises(RuntimeError, match=r"Unable to determine UTM CRS"):
+            xds.rio.estimate_utm_crs()