Skip to content

Commit

Permalink
Merge 3e5c3f1 into 8605ffe
Browse files Browse the repository at this point in the history
  • Loading branch information
@snowman2
snowman2 committed Jun 25, 2020
2 parents 8605ffe + 3e5c3f1 commit b0200f7
Show file tree
Hide file tree
Showing 5 changed files with 133 additions and 55 deletions.
3 changes: 3 additions & 0 deletions docs/history.rst
View file
Expand Up @@ -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
-------
Expand Down
99 changes: 72 additions & 27 deletions rioxarray/rioxarray.py
View file
Expand Up @@ -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
Expand Down Expand Up @@ -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`
Expand All @@ -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")
Expand Down Expand Up @@ -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)
Expand Down Expand Up @@ -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)
Expand Down
2 changes: 1 addition & 1 deletion test/conftest.py
View file
Expand Up @@ -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"):
Expand Down
84 changes: 57 additions & 27 deletions test/integration/test_integration_rioxarray.py
View file
Expand Up @@ -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(
Expand Down Expand Up @@ -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"])

Expand All @@ -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
)


Expand All @@ -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
)


Expand Down
Binary file added test/test_data/compare/MODIS_ARRAY_CLIP_EXPAND.nc
View file
Binary file not shown.

0 comments on commit b0200f7

Please sign in to comment.