How to open a dataset to a wanted resolution

[remi-braun]

Hello,

What I want is to copy the behavior of rasterio.read(out_shape=(count, height, width)) in order to open the dataset to a wanted resolution.
My current workaround is to use open_rasterio then do a interp but it might exists something simpler.
I haven't found anything related to this use case.

Can you help me ? :)

Best Regards

[snowman2]

Currently, decimated reads aren't supported currently. But, I wouldn't turn down a PR with the new feature 😄

Resampling after loading data is currently the only was using rioxarray directly:
https://corteva.github.io/rioxarray/stable/rioxarray.html#rioxarray.raster_array.RasterArray.reproject

rds = rioxarray.open_rasterio(...)
rds = rds.rio.reproject(rds.rio.crs, shape=(height, width))

[remi-braun]

Okay, my bad !
I wasn't looking for a reprojection, but rather a resampling haha

This way seems as fast as rasterio way so it does not bother me
Maybe a PR is too much ? But a line in the documentation would be very nice 😉

[snowman2]

Sounds like a good idea. If you would like to add an example to the documentation, that would be helpful (definitely don't have to).

[remi-braun]

Done 😉

[snowman2]

Thanks 👍

[remi-braun]

Just to add something here, it seems that the rasterio's way and the rioxarray's way have some discrepancy.
You can load any satellite data, the method will give the same result without any reprojection, but will give different values when resampled.

For example, in CI I load L8 data with a very low resolution of 1000 meters and I get different values:

AssertionError: 
Arrays are not almost equal to 5 decimals

Mismatched elements: 33693 / 53820 (62.6%)
Max absolute difference: 0.32614
Max relative difference: 1.6444925
 x: array([[0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.],...
 y: array([[0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.],
       [0., 0., 0., ..., 0., 0., 0.],...

This is also the case for S2, S2-theia and all landsat satellites.
I don't know from where it comes from but this is a bit annoying 😓

[snowman2]

rio.reproject is powered by rasterio.warp.reproject.

[remi-braun]

Yes but their resampling from the read function is made differently.
This was just in order to point out the difference between these two ways of loading resampled data

[cordmaur]

Hi @remi-braun , @snowman2
I don't know if that helps, but I've been using the following function as a workaround to the lack of decimated reading on rioxarray and it greatly improves my performance to create some thumbs.

def decimated_reading(file, out_shape):

    # Open the dataset with rasterio
    ds = rio.open(file)

    # Read it with the desired output shape (decimated reading)
    data=ds.read(out_shape=out_shape).squeeze()
    data[data == ds.nodata] = np.nan

    # create the lats and longs coordinates, considering the bounds and the shape
    longs = np.linspace(ds.bounds.left, ds.bounds.right, out_shape[1])
    lats = np.linspace(ds.bounds.top, ds.bounds.bottom, out_shape[0])

    arr = xr.DataArray(data=data, dims=['y', 'x'], coords=[lats, longs])

    # set CRS and GeoTransform
    arr.rio.set_crs(ds.crs)
    transform = rio.transform.from_bounds(*ds.bounds, height=out_shape[0], width=out_shape[1])
    arr = arr.rio.write_transform(transform)

    return arr

This opens the downscaled version of a 2.5Gb (42k x 36k) DEM in less than 1 sec compared to 35s by using rasterio.reproject.

arr = decimated_reading('/data/anadem_v018/anadem_south_2.tif', out_shape=(420, 360))
arr.plot(vmin=0, aspect=1, figsize=(5, 5))

Wondering if it would be possible to create a PR, but the problem is that the spatial_ref is not set... it would be necessary to reproduce part of the original open_rasterio.
Do you think it is worth working on that for a PR?

[snowman2]

but the problem is that the spatial_ref is not set...

This should help:

arr.rio.write_crs(...)

[snowman2]

If you wanted to start by adding a code example to the documentation, that could be helpful.

[remi-braun]

This looks really nice ! I will try it.
However, I am still wondering how rasterio does its decimated reading with other resamplings than nearest 🤔

[Kirill888]

there is resampling= parameter in the DatasetReader.read, would be nice if those could be exposed via rioxarray.open, I believe currently only .open arguments can be injected.

https://rasterio.readthedocs.io/en/latest/api/rasterio.io.html#rasterio.io.DatasetReader.read

[snowman2]

would be nice if those could be exposed via rioxarray.open

Interesting idea. I think that should be possible.