From 1af7a06998c5d26caaccf4e8a169af42a680cbec Mon Sep 17 00:00:00 2001 From: snowman2 Date: Wed, 24 Jun 2020 12:20:12 -0500 Subject: [PATCH 1/3] ENH: Use rasterio windows in rio.clip_box --- docs/history.rst | 1 + rioxarray/rioxarray.py | 99 +++++++++++++----- test/conftest.py | 2 +- .../integration/test_integration_rioxarray.py | 84 ++++++++++----- .../compare/MODIS_ARRAY_CLIP_EXPAND.nc | Bin 0 -> 13620 bytes 5 files changed, 131 insertions(+), 55 deletions(-) create mode 100644 test/test_data/compare/MODIS_ARRAY_CLIP_EXPAND.nc diff --git a/docs/history.rst b/docs/history.rst index b3718f39..04c44f1d 100644 --- a/docs/history.rst +++ b/docs/history.rst @@ -5,6 +5,7 @@ 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) 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..a704b512 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.y[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 0000000000000000000000000000000000000000..dd16044725e9fe68bf026392779e0c1c244fda54 GIT binary patch literal 13620 zcmeHOZ)jWB6+ch@Q^%GQrFOhDq?@8`G#^w=e{I<;ZDjqi>?o3u@_K&9YZYR;V2K(&czT2s*OL&{Gz`_@dqjTH#==h0VfL@f>2k}g-Sj2lm zg=BD464NJ>+2V9AD;DFiOcI<3nSh@j#+8koOd>JkqB(J{4RD45X!-Vn6Dh?kgEX8_ zYDyXDTsQ?oVK?5j_~_c1$N3=V01U%H5K~IEdP=IU$<)^wD{AKma}9V+eE^?5?3|<# z1?XpbS2cB6RHd3Jl<9v6={VPhdr8yuvZ6{2*{xOFtGf1pywpJRaNY|zcWzVlhNPAy zz0741dS7dk)8Z31#FohN6FJ2&kf89z4CHMmjz$BTns09)AKX+1vc5|Lp=lni9D8~v zBO;m!m-g0k|8!h@WU?50d{p>V3hVQs^KU)++Wc?y+T$NhVNJMEc_02Re2gFC&i=L$ zuo18kuo18kuo3uwMj#hYpD^6$Idn#;a5ws9a52}vd>3=Q&BaW=I{OpD#YEqb1q|H3 z>0&VidEiDIaU6PF+`8%F$s^E zlM)U9{DpA|Zvc!9FfqCSPxaE(JEJeZpGqeto0u@Thd#XYJ%NQpeb6=zF|Z)y9AeEE zDqqOW#0$lF&zgE#)y}FO?_{i)oJ!`VtgVSyY4&6rRWg&A&KGmjiMF=U&{)9lAM*vG z0k7YiNoJ=?>DlabY2H(l8@lZA2E0U{%onFT-f_P_uzBorJEMgrsG;ZSIyO zgp^dR%Whe!HxzZ5WOhQwGM}h_w37%`(;F-Qxr{MFI9{BOOL{}DE0XFBE_h>e)41SG zW$412DD^<&v5-F+^7(^dKbk-&;tqxazR}2NC=v~hhQg6ZDC9;9aYy{&=y=5M3k9R0(Lgv9 zV6AuGNKM1s5+#+MZu|Yylr|3=6We5IA%OfoWWpEtL90j_TTb#AhhPM}7(d+kS# z^3hFmk90K$zf^wr8l&1nRJ5-m)dUsa;ibt^CWYFSNCs~ERGqjy=Z5F$3KKdsy*I!0=@IyZQj(Z|B7_p( z|M3jmNeyEFd$77RT+hOvi545^?E71HY8LJxw2_5o^v=O|Pr*y1I5v=4*A`>*aE4gM z0TXQi&pbAx!BZr35ir9Aa9`x3XW>=Kxc&furtrkj^`8Pej8F!N`uV3s zdL%M?OE24(HaV*wM!%fbx;{QFz+skRrerX4F+{B*uIfrn7HhJul<`Pr8gxp}Q^(FJ zXW9UUahvD^t4N+IHvB!>SNKWV^s+ePvaBuZ(&~z`BrfT7D)}e) zn+~1;-{^!p&#;SzcemSAKaIGT$*+)OCz`$#YO?!v+}9?L`&GM_N9fq)e~Ezg_fI3o zqBBnl%zpcTep_cxH|+S!j=yfI_{+t1ow{68aDRP3(?w}rssGozSns3Hy!pb(xtM@H zHa$=^eOU%B=x~wkY4qmM)t7Mz) z$IEudQpatjnFynB0R)F}CYb>SbF!&@>wPOxX7sqAr=N~9N|dpd3*Ei!o|rEfna1yW z=my}pdF=7jSWh!L<53!~DHlm+)-VPg__*7AP#$kSey&sZjStMXd7*%4+0ZV!+FI^E z`0_qVgz*;{c)zf5!qm6dPi;as`AAgm6vcH(*QE_|B_!kOjHD~l{Z&~MH{FZ~=0*;k zAh|q2_z#nvLJ}(U@SR>77%=9z<&i$6B00iDy{%t$!5Jb-o|$O(Jcy6{8w#n+skG_kYpmz~)!`}k#PK=8b9SRtkxo~w_FGFOeizR?Na5Ms%9{JN;x#q@f_^b?#9{MH{ zV4*(Jkl#}XIzBV;xap!|vTMJz)J{MUtF)3amBK{PN2QAnneg=99K z6pb9j2qpw52SNFny4?7;n)6&m8J}~h^bap42hSU9C*X#E*)TIuJWs)Eo)hA7Uc7SG P4{x=(@wr}{`K9AuvVU4l literal 0 HcmV?d00001 From 880201516ccfbbe5311156239bc70385001871a7 Mon Sep 17 00:00:00 2001 From: snowman2 Date: Wed, 24 Jun 2020 15:54:50 -0500 Subject: [PATCH 2/3] update history --- docs/history.rst | 2 ++ 1 file changed, 2 insertions(+) diff --git a/docs/history.rst b/docs/history.rst index 04c44f1d..f906f283 100644 --- a/docs/history.rst +++ b/docs/history.rst @@ -6,6 +6,8 @@ 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 ------- From 3e5c3f1be317bc930b0720b3b82466b095a39f47 Mon Sep 17 00:00:00 2001 From: snowman2 Date: Thu, 25 Jun 2020 07:55:07 -0500 Subject: [PATCH 3/3] fix --- test/integration/test_integration_rioxarray.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/integration/test_integration_rioxarray.py b/test/integration/test_integration_rioxarray.py index a704b512..31ab44dd 100644 --- a/test/integration/test_integration_rioxarray.py +++ b/test/integration/test_integration_rioxarray.py @@ -326,7 +326,7 @@ def test_clip_box__one_dimension_error(modis_clip): xdi.isel(x=slice(5, 6), y=slice(5, 6)).rio.clip_box( minx=-7272735.53951584, # xdi.x[5].values miny=5048371.187465771, # xdi.y[7].values - maxx=-7272272.226799311, # xdi.y[7].values + maxx=-7272272.226799311, # xdi.x[7].values maxy=5048834.500182299, # xdi.y[5].values )