diff --git a/docs/history.rst b/docs/history.rst index c381976c..729eb799 100644 --- a/docs/history.rst +++ b/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 ------ diff --git a/rioxarray/_spatial_utils.py b/rioxarray/_spatial_utils.py index c651f012..0e00e914 100644 --- a/rioxarray/_spatial_utils.py +++ b/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]: diff --git a/rioxarray/rioxarray.py b/rioxarray/rioxarray.py index f3509a43..580abe1b 100644 --- a/rioxarray/rioxarray.py +++ b/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( diff --git a/test/integration/test_integration_rioxarray.py b/test/integration/test_integration_rioxarray.py index 7f8ed05b..7817e6de 100644 --- a/test/integration/test_integration_rioxarray.py +++ b/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", [