BUG: Fix rio.bounds() for rotation / asymmetric shear in the transform (#847)

docs/history.rst

@@ -4,6 +4,7 @@ History
 Unreleased
 ----------
 - ENH: Add write support for Zarr spatial and proj conventions
+- BUG: Fix `rio.bounds()` for grids with rotation or asymmetric shear in the transform (#847)
 
 0.22.0
 ------

rioxarray/_spatial_utils.py

@@ -58,7 +58,7 @@ def _affine_has_rotation(affine: Affine) -> bool:
     -------
     bool
     """
-    return affine.b == affine.d != 0
+    return affine.b != 0 or affine.d != 0
 
 
 def _resolution(affine: Affine) -> tuple[float, float]:
@@ -89,6 +89,41 @@ def _resolution(affine: Affine) -> tuple[float, float]:
     )
 
 
+def _rotated_bounds(
+    affine: Affine, *, width: int, height: int
+) -> tuple[float, float, float, float]:
+    """
+    Compute the axis-aligned bounds of a grid whose affine has rotation/shear.
+
+    The bounds are the envelope of the four corners of the grid transformed
+    by the affine. This matches how :obj:`rasterio.io.DatasetReader.bounds`
+    handles a non-north-up transform.
+
+    Parameters
+    ----------
+    affine: :obj:`affine.Affine`
+        The affine of the grid.
+    width: int
+        The width of the grid.
+    height: int
+        The height of the grid.
+
+    Returns
+    -------
+    left, bottom, right, top: float
+        Outermost coordinates of the grid.
+    """
+    corners = (
+        affine * (0, 0),
+        affine * (width, 0),
+        affine * (0, height),
+        affine * (width, height),
+    )
+    x_coords = [corner[0] for corner in corners]
+    y_coords = [corner[1] for corner in corners]
+    return min(x_coords), min(y_coords), max(x_coords), max(y_coords)
+
+
 def affine_to_coords(
     affine: Affine, width: int, height: int, *, x_dim: str = "x", y_dim: str = "y"
 ) -> dict[str, numpy.ndarray]:

rioxarray/rioxarray.py

@@ -36,6 +36,7 @@
     _has_spatial_dims,
     _order_bounds,
     _resolution,
+    _rotated_bounds,
     affine_to_coords,
 )
 from rioxarray.crs import crs_from_user_input
@@ -840,6 +841,11 @@ def bounds(self, *, recalc: bool = False) -> tuple[float, float, float, float]:
         left, bottom, right, top: float
             Outermost coordinates of the `xarray.DataArray` | `xarray.Dataset`.
         """
+        transform = self._cached_transform()
+        if transform is not None and _affine_has_rotation(transform):
+            # a rotated/sheared grid is not axis-aligned, so the bounds are
+            # the envelope of the four (transformed) corners of the grid.
+            return _rotated_bounds(transform, width=self.width, height=self.height)
         minx, miny, maxx, maxy = self._unordered_bounds(recalc=recalc)
         resolution_x, resolution_y = self.resolution(recalc=recalc)
         return _order_bounds(

test/integration/test_integration_rioxarray.py

@@ -21,7 +21,11 @@
 from rasterio.windows import Window
 
 import rioxarray
-from rioxarray._spatial_utils import _generate_spatial_coords, _make_coords
+from rioxarray._spatial_utils import (
+    _affine_has_rotation,
+    _generate_spatial_coords,
+    _make_coords,
+)
 from rioxarray.exceptions import (
     DimensionMissingCoordinateError,
     MissingCRS,
@@ -3313,6 +3317,58 @@ def test_bounds__ordered__dataset():
     assert xds.rio.bounds() == (-0.5, -0.5, 4.5, 4.5)
 
 
+@pytest.mark.parametrize(
+    "affine",
+    [
+        Affine(1, 0.3, 0, 0.1, -1, 0),  # asymmetric shear (b != d), see gh #847
+        Affine(1, 0.3, 0, 0.3, -1, 0),  # symmetric shear (b == d)
+        Affine(1, 0.2, 0, 0, -1, 0),  # single off-diagonal
+    ],
+)
+def test_bounds__rotated(affine):
+    # https://github.com/corteva/rioxarray/issues/847
+    # a sheared/rotated affine (b or d nonzero) is detected as having rotation ...
+    assert _affine_has_rotation(affine)
+    width, height = 4, 5
+    image = numpy.zeros((height, width), dtype="uint8")
+    xds = (
+        xarray.DataArray(
+            image,
+            dims=("y", "x"),
+            coords=_generate_spatial_coords(affine=affine, width=width, height=height),
+        )
+        .rio.write_crs("EPSG:4326", inplace=True)
+        .rio.write_transform(affine, inplace=True)
+        .rio.write_coordinate_system(inplace=True)
+    )
+    # ... and the bounds are the axis-aligned envelope of the 4 corners
+    corners = [
+        affine * (0, 0),
+        affine * (width, 0),
+        affine * (0, height),
+        affine * (width, height),
+    ]
+    xs = [corner[0] for corner in corners]
+    ys = [corner[1] for corner in corners]
+    expected = (min(xs), min(ys), max(xs), max(ys))
+    assert_almost_equal(xds.rio.bounds(), expected)
+    # cross-check against rasterio's own bounds for the same transform
+    with rasterio.io.MemoryFile() as memfile:
+        with memfile.open(
+            driver="GTiff",
+            height=height,
+            width=width,
+            count=1,
+            dtype="uint8",
+            crs="EPSG:4326",
+            transform=affine,
+        ) as dataset:
+            dataset.write(image, 1)
+        with memfile.open() as dataset:
+            rasterio_bounds = dataset.bounds
+    assert_almost_equal(xds.rio.bounds(), tuple(rasterio_bounds))
+
+
 @pytest.mark.parametrize(
     "rename",
     [