BUG: Raise informative error for 2D spatial coordinates and fix rio.shape (#848)

docs/history.rst

@@ -4,6 +4,11 @@ History
 Latest
 ------
 
+- BUG: Raise informative ``CoordinateDimensionError`` when spatial coordinates
+  are not 1-dimensional and fix ``rio.width``/``rio.height`` with 2D
+  coordinates, which also fixes the output size of ``rio.reproject`` with
+  ``src_geoloc_array`` (issue #848)
+
 0.23.0
 ------
 - ENH: Add write support for Zarr spatial and proj conventions

rioxarray/exceptions.py

@@ -39,5 +39,9 @@ class DimensionMissingCoordinateError(RioXarrayError):
     """This is raised when the dimension does not have the supporting coordinate."""
 
 
+class CoordinateDimensionError(RioXarrayError):
+    """This is raised when the coordinate of a spatial dimension is not 1-dimensional."""
+
+
 class MissingCRS(RioXarrayError):
     """Missing the CRS in the dataset."""

rioxarray/rioxarray.py

@@ -41,6 +41,7 @@
 from rioxarray.crs import crs_from_user_input
 from rioxarray.enum import Convention
 from rioxarray.exceptions import (
+    CoordinateDimensionError,
     DimensionError,
     DimensionMissingCoordinateError,
     InvalidDimensionOrder,
@@ -670,15 +671,15 @@ def width(self) -> int:
         """int: Returns the width of the dataset (x dimension size)"""
         if self._width is not None:
             return self._width
-        self._width = self._obj[self.x_dim].size
+        self._width = self._obj.sizes[self.x_dim]
         return self._width
 
     @property
     def height(self) -> int:
         """int: Returns the height of the dataset (y dimension size)"""
         if self._height is not None:
             return self._height
-        self._height = self._obj[self.y_dim].size
+        self._height = self._obj.sizes[self.y_dim]
         return self._height
 
     @property
@@ -736,6 +737,16 @@ def _internal_bounds(self) -> tuple[float, float, float, float]:
             raise DimensionMissingCoordinateError(f"{self.x_dim} missing coordinates.")
         if self.y_dim not in self._obj.coords:
             raise DimensionMissingCoordinateError(f"{self.y_dim} missing coordinates.")
+        for dim in (self.x_dim, self.y_dim):
+            if self._obj[dim].ndim != 1:
+                raise CoordinateDimensionError(
+                    f"The {dim} coordinate must be 1-dimensional to determine "
+                    f"the bounds, resolution, or transform; it has "
+                    f"{self._obj[dim].ndim} dimensions. To reproject data on a "
+                    f"curvilinear grid, pass the 2D coordinates to 'reproject' "
+                    f"using the 'src_geoloc_array' keyword argument."
+                    f"{_get_data_var_message(self._obj)}"
+                )
         try:
             left = float(self._obj[self.x_dim][0])
             right = float(self._obj[self.x_dim][-1])

test/integration/test_integration_rioxarray.py

@@ -18,11 +18,13 @@
 from rasterio import DatasetReader
 from rasterio.control import GroundControlPoint
 from rasterio.crs import CRS
+from rasterio.warp import calculate_default_transform
 from rasterio.windows import Window
 
 import rioxarray
 from rioxarray._spatial_utils import _generate_spatial_coords, _make_coords
 from rioxarray.exceptions import (
+    CoordinateDimensionError,
     DimensionMissingCoordinateError,
     MissingCRS,
     MissingSpatialDimensionError,
@@ -2895,6 +2897,68 @@ def test_shape_order():
         assert rds.air_temperature.rio.shape == (585, 1386)
 
 
+def _two_dimensional_coordinate_dataarray():
+    yy, xx = numpy.meshgrid(
+        numpy.linspace(-30.0, -40.0, 5),
+        numpy.linspace(150.0, 155.0, 6),
+        indexing="ij",
+    )
+    return xarray.DataArray(
+        numpy.zeros((5, 6)),
+        dims=("y", "x"),
+        coords={"x": (("y", "x"), xx), "y": (("y", "x"), yy)},
+    ).rio.write_crs("EPSG:4326")
+
+
+def test_bounds__2d_coordinates():
+    xds = _two_dimensional_coordinate_dataarray()
+    with pytest.raises(CoordinateDimensionError, match="must be 1-dimensional"):
+        xds.rio.bounds()
+
+
+def test_resolution__2d_coordinates():
+    xds = _two_dimensional_coordinate_dataarray()
+    with pytest.raises(CoordinateDimensionError, match="must be 1-dimensional"):
+        xds.rio.resolution()
+
+
+def test_transform__2d_coordinates():
+    xds = _two_dimensional_coordinate_dataarray()
+    with pytest.raises(CoordinateDimensionError, match="must be 1-dimensional"):
+        xds.rio.transform(recalc=True)
+
+
+def test_reproject__2d_coordinates():
+    xds = _two_dimensional_coordinate_dataarray()
+    with pytest.raises(CoordinateDimensionError, match="must be 1-dimensional"):
+        xds.rio.reproject("EPSG:2193")
+
+
+def test_reproject__2d_coordinates__src_geoloc_array():
+    xds = _two_dimensional_coordinate_dataarray()
+    src_geoloc_array = (xds.x.values, xds.y.values)
+    _, width, height = calculate_default_transform(
+        "EPSG:4326",
+        "EPSG:2193",
+        xds.rio.width,
+        xds.rio.height,
+        src_geoloc_array=src_geoloc_array,
+    )
+    reprojected = xds.rio.reproject("EPSG:2193", src_geoloc_array=src_geoloc_array)
+    assert reprojected.rio.shape == (height, width)
+
+
+def test_shape__2d_coordinates():
+    xds = _two_dimensional_coordinate_dataarray()
+    assert xds.rio.shape == (5, 6)
+
+
+def test_bounds__2d_coordinates__dataset():
+    xds = xarray.Dataset({"var": _two_dimensional_coordinate_dataarray()})
+    with pytest.raises(CoordinateDimensionError, match="must be 1-dimensional"):
+        xds.rio.bounds()
+
+
 def test_write_transform__from_read(tmp_path):
     with rioxarray.open_rasterio(
         os.path.join(TEST_COMPARE_DATA_DIR, "small_dem_3m_merged.tif"),