extracting pixel data for individual polygons from a raster using rio.clip()

[JessicaS11]

I've got an Xarray dataset (ds) of geospatial data (DEMs) and a companion Geopandas geodataframe (gdf) with a geometry column of polygons. I'd like to use the polygons to extract pixel data from the DEMs and have it associated with each row in the geodataframe. I can readily use ds['var'].rio.clip(gdf.geometry) to create a new raster clipped to the geometry data, but what is the best way to associate the raster features back to their respective geometries? I've tried doing the clip per-geometry by iterating on the geodataframe:

for datarow in gdf.itertuples(index=True):
    def get_vals(datarow):
            med = np.nanmedian(ds['var'].rio.clip([datarow.geometry]).values)
        return med
    gdf['med'] = gdf.apply(get_vals, axis=1)

but this is quite slow given the size of my raster. In other portions of my code, I've been able to speed up this process by trimming the raster to a buffered bounding box of the geometry prior to clipping. I again need that ability and was about to start turning my existing implementation into a more generalized function. I'm wondering if there's a relevant place/space within rioxarray to house something like this so that I can refactor into a more public space than my own repo.

[snowman2]

I think you are looking for geocube: https://corteva.github.io/geocube/stable/examples/zonal_statistics.html

[JessicaS11]

Thanks @snowman2! Exploring geocube was on my to-do list this week, so this gave me the perfect opportunity to do so. My original implementation ended up being ~3x faster in the context of my workflow (I'm assuming the geocube approach will scale better though), and I'm excited to use geocube more for vector-raster work.