diff --git a/docs/history.rst b/docs/history.rst index b3718f39..f906f283 100644 --- a/docs/history.rst +++ b/docs/history.rst @@ -5,6 +5,9 @@ Latest ------ - BUG: Fix assigning fill value in `rio.pad_box` (pull #140) - ENH: Add `rio.write_transform` to store cache in GDAL location (issue #129 & #139) +- ENH: Use rasterio windows for `rio.clip_box` (issue #142) +- BUG: Add support for negative indexes in rio.isel_window (pull #145) +- BUG: Write transform based on window in rio.isel_window (pull #145) 0.0.29 ------- diff --git a/rioxarray/rioxarray.py b/rioxarray/rioxarray.py index cdd81518..68319acf 100644 --- a/rioxarray/rioxarray.py +++ b/rioxarray/rioxarray.py @@ -17,12 +17,12 @@ import numpy as np import pyproj import rasterio.warp +import rasterio.windows import xarray from affine import Affine from rasterio.crs import CRS from rasterio.enums import Resampling from rasterio.features import geometry_mask -from rasterio.windows import get_data_window from scipy.interpolate import griddata from rioxarray.crs import crs_to_wkt @@ -642,6 +642,8 @@ def isel_window(self, window): """ Use a rasterio.window.Window to select a subset of the data. + .. warning:: Float indices are converted to integers. + Parameters ---------- window: :class:`rasterio.window.Window` @@ -652,11 +654,29 @@ def isel_window(self, window): :obj:`xarray.Dataset` | :obj:`xarray.DataArray`: The data in the window. """ (row_start, row_stop), (col_start, col_stop) = window.toranges() - row_slice = slice(int(math.floor(row_start)), int(math.ceil(row_stop))) - col_slice = slice(int(math.floor(col_start)), int(math.ceil(col_stop))) - return self._obj.isel( - {self.y_dim: row_slice, self.x_dim: col_slice} - ).rio.set_spatial_dims(x_dim=self.x_dim, y_dim=self.y_dim, inplace=True) + row_start = math.ceil(row_start) if row_start < 0 else math.floor(row_start) + row_stop = math.floor(row_stop) if row_stop < 0 else math.ceil(row_stop) + col_start = math.ceil(col_start) if col_start < 0 else math.floor(col_start) + col_stop = math.floor(col_stop) if col_stop < 0 else math.ceil(col_stop) + row_slice = slice(int(row_start), int(row_stop)) + col_slice = slice(int(col_start), int(col_stop)) + return ( + self._obj.isel({self.y_dim: row_slice, self.x_dim: col_slice}) + .copy() # this is to prevent sharing coordinates with the original dataset + .rio.set_spatial_dims(x_dim=self.x_dim, y_dim=self.y_dim, inplace=True) + .rio.write_transform( + transform=rasterio.windows.transform( + rasterio.windows.Window.from_slices( + rows=row_slice, + cols=col_slice, + width=self.width, + height=self.height, + ), + self.transform(recalc=True), + ), + inplace=True, + ) + ) @xarray.register_dataarray_accessor("rio") @@ -1264,33 +1284,56 @@ def clip_box(self, minx, miny, maxx, maxy, auto_expand=False, auto_expand_limit= f"{_get_data_var_message(self._obj)}" ) + # make sure that if the coordinates are + # in reverse order that it still works resolution_x, resolution_y = self.resolution() + if resolution_y < 0: + top = maxy + bottom = miny + else: + top = miny + bottom = maxy + if resolution_x < 0: + left = maxx + right = minx + else: + left = minx + right = maxx + + # pull the data out + window = rasterio.windows.from_bounds( + left=np.array(left).item(), + bottom=np.array(bottom).item(), + right=np.array(right).item(), + top=np.array(top).item(), + transform=self.transform(recalc=True), + width=self.width, + height=self.height, + ) + cl_array = self.isel_window(window) - clip_minx = minx - abs(resolution_x) / 2.0 - clip_miny = miny - abs(resolution_y) / 2.0 - clip_maxx = maxx + abs(resolution_x) / 2.0 - clip_maxy = maxy + abs(resolution_y) / 2.0 - cl_array = self.slice_xy(clip_minx, clip_miny, clip_maxx, clip_maxy) - if cl_array.rio.width < 1 or cl_array.rio.height < 1: - raise NoDataInBounds( - f"No data found in bounds.{_get_data_var_message(self._obj)}" - ) - - if cl_array.rio.width == 1 or cl_array.rio.height == 1: + # check that the window has data in it + if cl_array.rio.width <= 1 or cl_array.rio.height <= 1: if auto_expand and auto_expand < auto_expand_limit: + resolution_x, resolution_y = self.resolution() return self.clip_box( - clip_minx, - clip_miny, - clip_maxx, - clip_maxy, + minx=minx - abs(resolution_x) / 2.0, + miny=miny - abs(resolution_y) / 2.0, + maxx=maxx + abs(resolution_x) / 2.0, + maxy=maxy + abs(resolution_y) / 2.0, auto_expand=int(auto_expand) + 1, auto_expand_limit=auto_expand_limit, ) - raise OneDimensionalRaster( - "At least one of the clipped raster x,y coordinates" - " has only one point." - f"{_get_data_var_message(self._obj)}" - ) + if cl_array.rio.width < 1 or cl_array.rio.height < 1: + raise NoDataInBounds( + f"No data found in bounds.{_get_data_var_message(self._obj)}" + ) + elif cl_array.rio.width == 1 or cl_array.rio.height == 1: + raise OneDimensionalRaster( + "At least one of the clipped raster x,y coordinates" + " has only one point." + f"{_get_data_var_message(self._obj)}" + ) # make sure correct attributes preserved & projection added _add_attrs_proj(cl_array, self._obj) @@ -1372,7 +1415,9 @@ def clip(self, geometries, crs, all_touched=False, drop=True, invert=False): x_dim=self.x_dim, y_dim=self.y_dim, inplace=True ) cropped_ds = cropped_ds.rio.isel_window( - get_data_window(np.ma.masked_array(clip_mask_arr, ~clip_mask_arr)) + rasterio.windows.get_data_window( + np.ma.masked_array(clip_mask_arr, ~clip_mask_arr) + ) ) if self.nodata is not None and not np.isnan(self.nodata): cropped_ds = cropped_ds.fillna(self.nodata) diff --git a/test/conftest.py b/test/conftest.py index a0d00d89..9251253d 100644 --- a/test/conftest.py +++ b/test/conftest.py @@ -35,7 +35,7 @@ def _assert_attrs_equal(input_xr, compare_xr, decimal_precision): def _assert_xarrays_equal( - input_xarray, compare_xarray, precision=7, skip_xy_check=False + input_xarray, compare_xarray, precision=7, skip_xy_check=False, ): _assert_attrs_equal(input_xarray, compare_xarray, precision) if hasattr(input_xarray, "variables"): diff --git a/test/integration/test_integration_rioxarray.py b/test/integration/test_integration_rioxarray.py index 47badfcb..31ab44dd 100644 --- a/test/integration/test_integration_rioxarray.py +++ b/test/integration/test_integration_rioxarray.py @@ -146,18 +146,29 @@ def _del_attr(input_xr, attr): @pytest.fixture( - params=[xarray.open_dataset, xarray.open_dataarray, rioxarray.open_rasterio] + params=[ + xarray.open_dataset, + xarray.open_dataarray, + rioxarray.open_rasterio, + partial(rioxarray.open_rasterio, parse_coordinates=False), + ] ) def modis_clip(request, tmpdir): return dict( input=os.path.join(TEST_INPUT_DATA_DIR, "MODIS_ARRAY.nc"), compare=os.path.join(TEST_COMPARE_DATA_DIR, "MODIS_ARRAY_CLIP.nc"), + compare_expand=os.path.join( + TEST_COMPARE_DATA_DIR, "MODIS_ARRAY_CLIP_EXPAND.nc" + ), open=request.param, output=str(tmpdir.join("MODIS_CLIP_DUMP.nc")), ) def test_pad_box(modis_clip): + if isinstance(modis_clip["open"], partial): + # SKIP: parse_coodinates=False is not supported + return with modis_clip["open"](modis_clip["input"]) as xdi: # first, clip clipped_ds = xdi.rio.clip_box( @@ -227,30 +238,50 @@ def test_clip_box(modis_clip): modis_clip["compare"] ) as xdc: clipped_ds = xdi.rio.clip_box( - minx=xdi.x[4].values, - miny=xdi.y[6].values, - maxx=xdi.x[6].values, - maxy=xdi.y[4].values, + minx=-7272967.195874103, # xdi.x[4].values, + miny=5048602.8438240355, # xdi.y[6].values, + maxx=-7272503.8831575755, # xdi.x[6].values, + maxy=5049066.156540562, # xdi.y[4].values, ) - _assert_xarrays_equal(clipped_ds, xdc) assert xdi.rio._cached_transform() != clipped_ds.rio._cached_transform() + var = "__xarray_dataarray_variable__" + try: + clipped_ds_values = clipped_ds[var].values + except KeyError: + clipped_ds_values = clipped_ds.values + try: + xdc_values = xdc[var].values + except KeyError: + xdc_values = xdc.values + assert_almost_equal(clipped_ds_values, xdc_values) + assert_almost_equal(clipped_ds.rio.transform(), xdc.rio.transform()) # make sure it safely writes to netcdf clipped_ds.to_netcdf(modis_clip["output"]) def test_clip_box__auto_expand(modis_clip): with modis_clip["open"](modis_clip["input"]) as xdi, modis_clip["open"]( - modis_clip["compare"] + modis_clip["compare_expand"] ) as xdc: clipped_ds = xdi.rio.clip_box( - minx=xdi.x[5].values, - miny=xdi.y[5].values, - maxx=xdi.x[5].values, - maxy=xdi.y[5].values, + minx=-7272735.53951584, # xdi.x[5].values + miny=5048834.500182299, # xdi.y[5].values + maxx=-7272735.53951584, # xdi.x[5].values + maxy=5048834.500182299, # xdi.y[5].values auto_expand=True, ) - - _assert_xarrays_equal(clipped_ds, xdc) + assert xdi.rio._cached_transform() != clipped_ds.rio._cached_transform() + var = "__xarray_dataarray_variable__" + try: + clipped_ds_values = clipped_ds[var].values + except KeyError: + clipped_ds_values = clipped_ds.values + try: + xdc_values = xdc[var].values + except KeyError: + xdc_values = xdc.values + assert_almost_equal(clipped_ds_values, xdc_values) + assert_almost_equal(clipped_ds.rio.transform(), xdc.rio.transform()) # make sure it safely writes to netcdf clipped_ds.to_netcdf(modis_clip["output"]) @@ -261,14 +292,13 @@ def test_clip_box__nodata_error(modis_clip): if hasattr(xdi, "name") and xdi.name: var_match = " Data variable: __xarray_dataarray_variable__" with pytest.raises( - NoDataInBounds, - match=f"Unable to determine bounds from coordinates.{var_match}", + NoDataInBounds, match=var_match, ): xdi.rio.clip_box( - minx=xdi.x[5].values, - miny=xdi.y[7].values, - maxx=xdi.x[4].values, - maxy=xdi.y[5].values, + minx=-8272735.53951584, # xdi.x[5].values + miny=8048371.187465771, # xdi.y[7].values + maxx=-8272967.195874103, # xdi.x[4].values + maxy=8048834.500182299, # xdi.y[5].values ) @@ -286,18 +316,18 @@ def test_clip_box__one_dimension_error(modis_clip): ), ): xdi.rio.clip_box( - minx=xdi.x[5].values, - miny=xdi.y[5].values, - maxx=xdi.x[5].values, - maxy=xdi.y[5].values, + minx=-7272735.53951584, # xdi.x[5].values + miny=5048834.500182299, # xdi.y[5].values + maxx=-7272735.53951584, # xdi.x[5].values + maxy=5048834.500182299, # xdi.y[5].values ) # test exception before raster clipped with pytest.raises(OneDimensionalRaster): xdi.isel(x=slice(5, 6), y=slice(5, 6)).rio.clip_box( - minx=xdi.x[5].values, - miny=xdi.y[7].values, - maxx=xdi.x[7].values, - maxy=xdi.y[5].values, + minx=-7272735.53951584, # xdi.x[5].values + miny=5048371.187465771, # xdi.y[7].values + maxx=-7272272.226799311, # xdi.x[7].values + maxy=5048834.500182299, # xdi.y[5].values ) diff --git a/test/test_data/compare/MODIS_ARRAY_CLIP_EXPAND.nc b/test/test_data/compare/MODIS_ARRAY_CLIP_EXPAND.nc new file mode 100644 index 00000000..dd160447 Binary files /dev/null and b/test/test_data/compare/MODIS_ARRAY_CLIP_EXPAND.nc differ