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?

[brendancol]

@rsignell-usgs @alaws-USGS just wanted to chime in here, hope all are well. Like @ianthomas23 mentioned, xarray-spatial doesn't handle the rasterization step (yet...), but it does have zonal.apply() which takes an arbitrary callable on irregular zonals defined by contiguous values / labels.

Assuming you rasterized using datashader or another method...the xarray-spatial steps would look like:

  import numpy as np
  import xarray as xr
  from xrspatial import zonal

  # -- Assume rasterization is already done --
  # zones: 2D raster where each pixel holds a polygon/watershed ID
  # precip: 2D raster of precipitation values (same grid)
  # weights: 2D raster of fractional pixel overlap (0.0–1.0) from
  #          sub-pixel rasterization (e.g. datashader at 5x oversampling)

  # Area-weighted aggregation: multiply values by pixel overlap fractions
  weighted_precip = precip * weights

  # Compute per-zone stats (sum = total weighted precip, count, mean, etc.)
  results = zonal.stats(zones=zones, values=weighted_precip,
                        stats_funcs=['sum', 'mean', 'count'])

  # results is a DataFrame with one row per watershed:
  #   zone  |  sum   |  mean  |  count
  #   ------+--------+--------+-------
  #   1     |  42.3  |  0.85  |  50
  #   2     |  67.1  |  1.12  |  60

  For dask-backed rasters (large datasets that don't fit in memory):

  import dask.array as da

  # Chunk everything consistently
  zones_dask = zones.chunk({'y': 512, 'x': 512})
  weighted_dask = weighted_precip.chunk({'y': 512, 'x': 512})

  # Same API — returns a dask DataFrame, computed lazily
  results = zonal.stats(zones=zones_dask, values=weighted_dask,
                        stats_funcs=['sum', 'mean', 'count'])