Feature Request: 3D Zonal Stats

[rbavery]

In some cases you may have a time series of rasters that you would like to calculate zonal statistics for (standard deviation across the time series). Currently zonal statistics supports 2D rasters.

[brendancol]

@rbavery I love it and we should be able to support this without too much trouble...do you see this being a 3D xr.DataArray or a xr.Dataset with several layers?

[rbavery]

@brendancol late reply, apologies! I think in most cases users will want to compute 3D zonal stats on a time series of a single variable (reflectance or backscatter in a single channel). So I think a 3D xr.DataArray would make the most sense as an input, but that may be just because that's my preferred data structure for single variable xarrays.

For my own workflows if I have a multichannel raster timeseries, I'll typically store that as an xr.DataArray with a band dimension (especially if I'm using rioxarray), in addiiton to a time dimension. So it's easy enough to subset this with .sel to get a 3D xr.DataArray.

Here's an example carpentries lesson showing how I'm using this function for a 2D xarray.DataArray, I'd like to make a similar example on an NDVI Landsat timeseries: https://carpentries-incubator.github.io/geospatial-python/12-zonal-statistics-raster/index.html

[giovp]

+1 , would be great to have.

[brendancol]

@thuydotm Let's push this up in priority