diff --git a/docs/history.rst b/docs/history.rst index 38747916..f1c9dec1 100644 --- a/docs/history.rst +++ b/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 diff --git a/rioxarray/exceptions.py b/rioxarray/exceptions.py index 2e00b08c..f1c727a5 100644 --- a/rioxarray/exceptions.py +++ b/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.""" diff --git a/rioxarray/rioxarray.py b/rioxarray/rioxarray.py index f3509a43..b4b70b7c 100644 --- a/rioxarray/rioxarray.py +++ b/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,7 +671,7 @@ 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 @@ -678,7 +679,7 @@ 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]) diff --git a/test/integration/test_integration_rioxarray.py b/test/integration/test_integration_rioxarray.py index 245eb5be..07aadd3e 100644 --- a/test/integration/test_integration_rioxarray.py +++ b/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"),