Inconsistent behaviour of clip_box about rioxarray HOT 14 CLOSED

The the x-coordinates the same? If not, that is likely the reason.

from rioxarray.

If also the x-coordinate were the same, the two input rasters would overlap completely (as the MRE shows, that's not the case).

A simple diagram that I hope complements the MRE above:

from rioxarray.

from rioxarray.

len(da1_clipped.y),len(da2_clipped.y)

(22018, 22017)

This is true:

xr.testing.assert_equal(da1_clipped.y[1:], da2_clipped.y)

from rioxarray.

When clipping rasters, you cannot be guaranteed that everything will mach up exactly in your example here. Floating point issues and other things can have an impact on the behavior of the clip. If you want to join them together, I recommend using merge:

from rioxarray.merge import merge_arrays

da_merged = merge_arrays([da1_clipped, da2_clipped])

from rioxarray.

I'm very surprised. I was using clip_box as a shorter version of:

da1_clipped = da1.sel(x=slice(bounds[0], bounds[2]), y=slice(bounds[3], bounds[1]))
da2_clipped = da2.sel(x=slice(bounds[0], bounds[2]), y=slice(bounds[3], bounds[1]))
xr.testing.assert_equal(da1_clipped.y, da2_clipped.y) # no exception raised

Out of curiosity, why are the results of clip_box different?

from rioxarray.

The implementation is different. clip_box uses rasterio to determine the window to select:

numeric precision differences when calculating the transform can cause there to be a slight offset.

from rioxarray.

One part I forgot to mention is that the example you showed only considers the centroid of the grid cell whereas the rasterio version considers the entire grid cell.

from rioxarray.

I can see how that's one difference, but wouldn't that impact the "left" and "right" rasters (da1 and da2) in the same way? My main issue with this result from clip_box is its lack of symmetry.

from rioxarray.

wouldn't that impact the "left" and "right" rasters (da1 and da2) in the same way?

Not necessarily. Depending on how close the boundary is to the edge of the grid cells, floating point precision differences could make it include grid cells on one side versus the other. Also, the boundary may be slightly shifted compared to the raster. If the boundary or raster is shifted, there can be differences in the resulting raster.

from rioxarray.

If you want to be sure everything matches up exactly on both sides, one option is to merge the rasters before clipping.

from rioxarray.

I dug into the code and here are the rasterio Windows used for clipping:

da1:

Window(col_off=-20099.16, row_off=12617.999999999884, width=34670.16, height=22016.40000000014)
# this is converted to:
# row_start=12617.999999999884 row_stop=34634.40000000002 col_start=-20099.16 col_stop=14571.000000000004

da2:

Window(col_off=15900.84, row_off=12618.0, width=34670.16, height=22016.40000000014)
# this is converted to:
# row_start=12618.0 row_stop=34634.40000000014 col_start=15900.84 col_stop=50571.0

Note that da1 has row_start=12617.999999999884 . In rio.isel_window it uses math.floor to select the index of the row to ensure any intersecting cells are added. For da1 this converts row_start=12617 and for da2 this is row_start=12618.

And that is why there is an additional row for one and not the other. It all comes down to floating point precision issues.

from rioxarray.

Thanks for looking into this, any reason why math.floor has been chosen instead of round?

from rioxarray.

" it uses math.floor to select the index of the row to ensure any intersecting cells are added."

from rioxarray.