[Question] Saving a cropped tile

[guillemc23]

I've been thinking, for some use cases, saving the whole extension of the raster is not optimal. How do you manage these cases? I believe that, in this line, only some *kwargs get propagated:

eoreader/eoreader/utils.py

Line 524 in a8fac5c

def stack(band_xds: xr.Dataset, **kwargs) -> (xr.DataArray, type):

I think it would be great to allow passing a geometry argument that rakes care of cropping the output and saving only this cropped raster as in:

# Define AOI as a shapely polygon, wkt string or other type of geometry
aoi = shapely.Polygon("POLYGON (1 1, 1 0, 0 1, 0 0)
# Pass the area of interest as geometry argument so that it performs the cropping to this aoi before saving the stack
stack = prod.stack([NDVI, MNDWI, GREEN], pixel_size=20, geometry=aoi)

I guess loading the product would still need to load the whole product but maybe there's no need to save all of it 😄

[remi-braun]

This should already be supported with the window keyword.

You have a notebook on this topic here 😄

Why window you may ask? To mimic rasterio's keyword 😉

[guillemc23]

yes, this definitely makes sense, a bit less usable but will try to implement it soon. Thank you!

[remi-braun]

I didn't want to go into cropping to shapes as i wasn't sure about the real gain compared to the window reading. If you need that you can create a feature issue, but I don't promise anything ;)

[guillemc23]

I've been playing around with this... once you read the window, the output is no longer an S2Product, but instead it's a xarray.Dataset. How would you stack the result this way? I'm using sertit.rasters.write, but I'm not sure this is the way.

Also, it's cropping a bounding box around the GeoDataFrame I'm passing in as window, but my geometry is just the dataframe with the exact footprint I need 🤔

[remi-braun]

• prod.stack and prod.load both output xarray data structures (load = Dataset, stack = DataArray + possibility to write it on disk). The solution is just to use prod.stac(..., window=...). Product objects are like rasterio datasets, they do not have informations about the content of the bands so they don't care on what you read the band or at which resolution.
• window, as the keyword implies, only reads extent of your input geometry. I see you need a crop of a geometry, and this may be cheap to implement, so here is the corresponding issue: ENH: Add a 'geometry' keywords that mimics 'window' but also crops to the input geometry #226.

[guillemc23]

My bad, I couldn't explain it properly 😄 yeah, that's exactly what I was looking for, rasterio windows are by definition rectangular, but adding the logic to crop the output to a geometry would be great. I like the way eoreader handles product metadata and that's my main concern instead of just cropping it using rasterio.mask 😅