aggregate raster data to polygons with area-weighting

[rsignell-usgs]

A very common use case for USGS (and others in hydrology) is to aggregate raster data over some polygon (e.g. total precipitation over a watershed).

xagg is a package to aggregate gridded data in xarray to polygons in geopandas using area-weighting from the relative area overlaps between pixels and polygons. A raster of weights is generated for each polygon, then the weight raster is multiplied by the data raster and summed to get the value for the polygon.

Could xarray-spatial do this kind of calculation?

[ianthomas23]

For this sort of problem I would use datashader to convert the polygons to a raster of the appropriate area and resolution and use this raster as a mask in xarray-spatial. I am not sure that we'd want to accept vector data as inputs to xarray-spatial, but we'd certainly like to expand our documentation significantly to show how it fits into the wider ecosystem and how various other libraries can be used with xarray-spatial. Now that we can output vector data via polygonize, we could certainly show examples of vector to raster and back to vector again (or the other way round) with processing occurring on the most appropriate side (raster vs vector) for the problem in hand.

[rsignell-usgs]

@ianthomas23 , this would be the "point in polygon" approach, right (all weights 0 or 1)?

[ianthomas23]

That is what I meant, yes, but I realise now that doesn't give you the area weighting that you want.

[rsignell-usgs]

Although xagg creates polygons around each pixel (grid square) and then determines the fraction of these contained within the supplied polygon (e.g. watershed), it seems we could also just interpolate the gridded data to finer resolution and then use the point in polygon approach. Do you think that would be faster?

[ianthomas23]

I'll be specific, just in case we are talking at cross purposes. I would:

1. Work out the bounds and resolution (spacing) of the raster data.
2. Create a datashader.Canvas for the same bounds and resolution.
3. Make sure your polygon data is in the same coordinate reference system as your raster data.
4. Call Canvas.polygon to tag pixels in your canvas as 1 if they are inside the polygon(s) and 0 if not.
5. Multiply the returned xr.DataArray by your raster and sum (with some scaling factor if required) to get your "aggregate within polygon" value. So this doesn't even need any of xarray-spatial.
6. So far this is using full pixels, not fractional ones, from your polygon data. To get fractional ones simply scale up the shape of your datashader.Canvas by some integer factor, e.g. 5 in both dimensions, rerun the Canvas.polygon call and do a 5x5 blockwise average on the returned array to make it the same shape as your raster. Now the area of each raster pixel that is within a polygon has a resolution of 1/25=0.04.

This approach does not interpolate your gridded data, it just runs the polygon to raster conversion at a higher resolution. The benefit is that the work occurs in datashader which uses numba and should be really quick.

I should point out that I lean towards a datashader approach because I am used to it, it is not necessarily a better solution that using xagg. Evidently xagg is an excellent idea if you really do need accurate areas per pixel.

[rsignell-usgs]

@ianthomas23, yes, this is what the sort of thing I was imagining, determining area by counting a dense field of pixels, instead of calculations with overlapping geometry in geopandas. But much more clearly defined! I'll see if I can figure this out and report back!

[alaws-USGS]

Hello, I wanted to hop onto this discussion and ask if xarray-spatial had implemented any area weighted methods since this post?