# Conservative Region Aggregation with Xarray, Geopandas and Sparse
**URL:** https://discourse.pangeo.io/t/conservative-region-aggregation-with-xarray-geopandas-and-sparse/2715
**Category:** Uncategorized
**Created:** [September 7, 2022, 2:41pm UTC](https://discourse.pangeo.io/t/conservative-region-aggregation-with-xarray-geopandas-and-sparse/2715 "2022-09-07T14:41:26Z")
**Posts on this page:** 20
**Page:** 1
### Author:  [@rabernat](https://discourse.pangeo.io/u/rabernat)
#### Post date: [September 7, 2022, 2:41pm UTC](https://discourse.pangeo.io/t/conservative-region-aggregation-with-xarray-geopandas-and-sparse/2715/1 "2022-09-07T14:41:27Z")
If you just want the full notebook, it’s here: [https://notebooksharing.space/view/c6c1f3a7d0c260724115eaa2bf78f3738b275f7f633c1558639e7bbd75b31456](https://notebooksharing.space/view/c6c1f3a7d0c260724115eaa2bf78f3738b275f7f633c1558639e7bbd75b31456)
## Context
[Conservative regridding](https://climatedataguide.ucar.edu/climate-data-tools-and-analysis/regridding-overview) is an important and expensive computational operation in climate science. As opposed to nearest-neighbor interpolation or reprojecting from one CRS to another, conservative regridding needs to account for the full geometry of the source and target grid. Currently our main tool for conservative regridding is [xESMF](https://github.com/pangeo-data/xESMF/).
Although it is usually discussed in a different context, conservative regridding is funamentally similar to regional aggregation, as performed for example by the [xagg](https://xagg.readthedocs.io/en/latest/index.html) package; both methods require consideration of the cell / region geometry.
So to use GIS terminology, regridding is a many-to-many spatial join operation.
## Goals
We generally rely on packages for regridding. But none of the packages has the performance and features that users really need. (For example, xesmf still cannot run with dask distributed.) Moreover, xesmf is a heavy dependency; it won’t run on windows for example. There is an appetite for alternatives.
Meanwhile. Geopandas has recently gained a lot of great functionality and performance enhancements, meaning that is could be suitable for doing this join. I wanted to explore coding up this sort of regridding from scratch, without using a regridding package, and see how fast it could go.
**Goal:** Regrid a global precipitation dataset into countries _conservatively_, i.e. by exactly partitioning each grid cell into the precise region boundaries.
**Meta Goal:** Demonstrate that we don’t necessarily need a package for this workflow and showcase some of the new capabilities of GeoPandas and Xarray along the way.
My source dataset is the NASA GPCP dataset, via [Pangeo Forge](https://pangeo-forge.org/dashboard/feedstock/42) in Zarr format. It’s a global, daily precipitation dataset with 1x1 degree spatial resolution. The target regions are the [Natural Earth 50m admin boundaries](https://www.naturalearthdata.com/downloads/50m-cultural-vectors/), loaded via a shapefile.
## Approach
I take a three step approach:
- Represent both the original grid and target grid as GeoSeries with Polygon geometry
- Compute their area overlay and turn it into a sparse matrix
- Perform matrix multiplication on the full Xarray dataset (with a time dimension)
## Some Highlights
Check out the [full notebook](https://notebooksharing.space/view/c6c1f3a7d0c260724115eaa2bf78f3738b275f7f633c1558639e7bbd75b31456#displayOptions=) for details. Here are some highlights.
### The original dataset
```python
store = 'https://ncsa.osn.xsede.org/Pangeo/pangeo-forge/gpcp-feedstock/gpcp.zarr'
ds = xr.open_dataset(store, engine='zarr', chunks={})
```
### Creating shapely geometries from grid bounds
```python
points = grid.stack(point=("latitude", "longitude"))
boxes = xr.apply_ufunc(
bounds_to_poly,
points.lon_bounds,
points.lat_bounds,
input_core_dims=[("nv",), ("nv",)],
output_dtypes=[np.dtype('O')],
vectorize=True
)
boxes
```
### Converting that to a GeoDataframe
```python
grid_df= gp.GeoDataFrame(
data={"geometry": boxes.values, "latitude": boxes.latitude, "longitude": boxes.longitude},
index=boxes.indexes["point"],
crs=crs_orig
)
```
### Overlaying this with the region geometries and computing area weights
```auto
overlay = grid_df.overlay(regions_df)
grid_cell_fraction = (
overlay.geometry.area.groupby(overlay.SOVEREIGNT)
.transform(lambda x: x / x.sum())
)
```
### Turning this into a sparse Xarray dataset
```python
multi_index = overlay.set_index(["latitude", "longitude", "SOVEREIGNT"]).index
df_weights = pd.DataFrame({"weights": grid_cell_fraction.values}, index=multi_index)
ds_weights = xr.Dataset(df_weights)
weights_sparse = ds_weights.unstack(sparse=True, fill_value=0.).weights
```
### Applying the matrix multiplication
Note that we can’t just use `xr.dot` because of [einsum implementation · Issue #31 · pydata/sparse · GitHub](https://github.com/pydata/sparse/issues/31).
```python
def apply_weights_matmul_sparse(weights, data):
assert isinstance(weights, sparse.SparseArray)
assert isinstance(data, np.ndarray)
data = sparse.COO.from_numpy(data)
data_shape = data.shape
# k = nlat * nlon
n, k = data_shape[0], data_shape[1] * data_shape[2]
data = data.reshape((n, k))
weights_shape = weights.shape
k_, m = weights_shape[0] * weights_shape[1], weights_shape[2]
assert k == k_
weights_data = weights.reshape((k, m))
regridded = sparse.matmul(data, weights_data)
assert regridded.shape == (n, m)
return regridded.todense()
precip_regridded = xr.apply_ufunc(
apply_weights_matmul_sparse,
weights_sparse,
precip_in_mem,
join="left",
input_core_dims=[["latitude", "longitude", "SOVEREIGNT"], ["latitude", "longitude"]],
output_core_dims=[["SOVEREIGNT"]],
dask="parallelized",
meta=[np.ndarray((0,))]
)
# takes < 10s to regrid over 9000 timesteps
precip_regridded.load()
```
### Plot some data
```python
precip_regridded.sel(SOVEREIGNT="Italy").resample(time="MS").mean().plot()
```
## Thoughts
I think this is a promising way forward for regridding. It removes the xESMF dependency. The challenge will be scaling it up to ultra-high-resolution global grids with millions of points. I will explore that in a follow-up post.
---
### Author:  [@ThomasMGeo](https://discourse.pangeo.io/u/ThomasMGeo)
#### Post date: [September 7, 2022, 3:11pm UTC](https://discourse.pangeo.io/t/conservative-region-aggregation-with-xarray-geopandas-and-sparse/2715/2 "2022-09-07T15:11:24Z")
Hi Ryan, timely topic as I am doing a deep dive on gridding w/r/t xarray.
Any downsides to this compared to xESMF (besides more lines of code written)?
Wrestling with the no package vs package pro’s and con’s. Thanks for sharing this!
---
### Author:  [@rabernat](https://discourse.pangeo.io/u/rabernat)
#### Post date: [September 7, 2022, 3:13pm UTC](https://discourse.pangeo.io/t/conservative-region-aggregation-with-xarray-geopandas-and-sparse/2715/3 "2022-09-07T15:13:15Z")
I’m sure there are many downsides compared to xESMF. I have not attempted to enumerate the pros and cons yet. I’m still at the exploratory stage. If you’re interested in this topic, I’d love to hear your thoughts!
---
### Author:  [@shoyer](https://discourse.pangeo.io/u/shoyer)
#### Post date: [September 7, 2022, 4:12pm UTC](https://discourse.pangeo.io/t/conservative-region-aggregation-with-xarray-geopandas-and-sparse/2715/4 "2022-09-07T16:12:55Z")
For the case of interpolation between rectilinear grids (even on the sphere), you can factorize regridding along each axis. This is less general but makes the entire calculation much simpler, because its feasible to store interpolation weights as dense matrices and to use dense matrix multiplication.
I have some code doing this at scale with Xarray-Beam that we’ve been running on 0.25x0.25 degree inputs. We should be able to source open source it soon if there is interest.
---
### Author:  [@rmcd-mscb](https://discourse.pangeo.io/u/rmcd-mscb)
#### Post date: [September 7, 2022, 6:02pm UTC](https://discourse.pangeo.io/t/conservative-region-aggregation-with-xarray-geopandas-and-sparse/2715/5 "2022-09-07T18:02:56Z")
This work looks really great @rabernat ! We also have been working on a python package for what I’d call grid-to-polygon area weighted intersections or aggregations based on the methods presented here: [R-tree Spatial Indexing with Python – Geoff Boeing](https://geoffboeing.com/2016/10/r-tree-spatial-index-python/). The repo is available here: [Water Mission Area / nhgf / ToolsTeam / gdptools · GitLab (usgs.gov)](https://code.usgs.gov/wma/nhgf/toolsteam/gdptools). Still in development. There are 2 workflows available: 1) using Mike Johnsons openDAP catalog as a source for gridded data (see: [https://mikejohnson51.github.io/opendap.catalog/cat\_params.json](https://mikejohnson51.github.io/opendap.catalog/cat_params.json)) and 2) based on a user-defined gridded source. Hope to have more use-cases available soon in the “preliminary” documentation here: [Simple openDAP Catalog example. — gdptools](https://gdptools.readthedocs.io/en/latest/terraclime_et.html)
---
### Author:  [@Michael\_Sumner](https://discourse.pangeo.io/u/Michael_Sumner)
#### Post date: [September 8, 2022, 12:42am UTC](https://discourse.pangeo.io/t/conservative-region-aggregation-with-xarray-geopandas-and-sparse/2715/6 "2022-09-08T00:42:52Z")
thanks! @rabernat is the zarr url meant to be generally accessible? I can’t tell if my tooling is off , or if I need to be on a system?
---
### Author:  [@rabernat](https://discourse.pangeo.io/u/rabernat)
#### Post date: [September 8, 2022, 12:46am UTC](https://discourse.pangeo.io/t/conservative-region-aggregation-with-xarray-geopandas-and-sparse/2715/7 "2022-09-08T00:46:59Z")
Welcome Michael!
The data should be 100% public. They are stored on [Open Storage Network](https://www.openstoragenetwork.org/) and should be accessible anywhere on the internet. To verify you can access them, take python out of the loop and try running
```auto
$ curl https://ncsa.osn.xsede.org/Pangeo/pangeo-forge/gpcp-feedstock/gpcp.zarr/.zgroup
```
from the command line. You should see
```auto
{
"zarr_format": 2
}
```
---
### Author:  [@Michael\_Sumner](https://discourse.pangeo.io/u/Michael_Sumner)
#### Post date: [September 8, 2022, 1:00am UTC](https://discourse.pangeo.io/t/conservative-region-aggregation-with-xarray-geopandas-and-sparse/2715/8 "2022-09-08T01:00:47Z")
great thanks a lot! I’m totally flailing until I get my tooling practice in - this is a whole-example where I can see my way through to compare some stuff in R - I’m inspired to finally get my python env setup
---
### Author:  [@rabernat](https://discourse.pangeo.io/u/rabernat)
#### Post date: [September 8, 2022, 1:15am UTC](https://discourse.pangeo.io/t/conservative-region-aggregation-with-xarray-geopandas-and-sparse/2715/9 "2022-09-08T01:15:07Z")
Getting a new python environment set up unfortunately can be tricky. If you want something that just works out of the box, check out Pangeo Docker Images (you can use them locally or in the cloud)
> **[GitHub - pangeo-data/pangeo-docker-images: Docker Images For Pangeo...](https://github.com/pangeo-data/pangeo-docker-images)**
>
> Docker Images For Pangeo JupyterHubs and BinderHubs - GitHub - pangeo-data/pangeo-docker-images: Docker Images For Pangeo JupyterHubs and BinderHubs
---
### Author:  [@benbovy](https://discourse.pangeo.io/u/benbovy)
#### Post date: [September 8, 2022, 7:51am UTC](https://discourse.pangeo.io/t/conservative-region-aggregation-with-xarray-geopandas-and-sparse/2715/10 "2022-09-08T07:51:33Z")
Thanks for sharing this @rabernat. That’s the kind of use case I’d like to consider for vectorized Python bindings of [s2geometry](https://github.com/google/s2geometry) and/or [s2geography](https://github.com/paleolimbot/s2geography), to be eventually (hopefully) integrated with Geopandas (and Xarray too, via a geo-extension, why not?!).
This general approach is a nice complement to @shoyer’s special case for rectilinear grids. However, a solution based on shapely (GEOS) planar geometries is still sub-optimal for lat / lon grids where the data projection step may be expensive. Also, naive question: do the polygon shapes distorted by the equal-area projection, combined with their limited number of vertices (grid and/or region geometries) and/or their extent, have any significant effect on the accuracy of the computed area weights?
s2geometry has a lot of features that we could potentially leverage in Python (vectorized) bindings. For example, it provides a [GetOverlapFractions](https://github.com/google/s2geometry/blob/14f5d60875f79ffe275f332eccdf48bcac89a3c4/src/s2/s2polygon.h#L297-L300) function that seems well fitted for this purpose.
---
### Author:  [@benbovy](https://discourse.pangeo.io/u/benbovy)
#### Post date: [September 8, 2022, 8:14am UTC](https://discourse.pangeo.io/t/conservative-region-aggregation-with-xarray-geopandas-and-sparse/2715/11 "2022-09-08T08:14:09Z")
If it is acceptable to perform conservative aggregation of arbitrary polygons at a given finite resolution, I think that something based on s2geometry (or H3) cells may offer great speed-ups. Both libraries support large to very fine spatial resolutions (cf. @darothen’s [comment](https://discourse.pangeo.io/t/discrete-global-grid-systems-dggs-use-with-pangeo/2274/22)).
Like xESMF, a solution based on s2geometry or h3 would also introduce some heavy dependencies, but not heavier than GEOS used by shapely.
---
### Author:  [@rabernat](https://discourse.pangeo.io/u/rabernat)
#### Post date: [September 8, 2022, 12:50pm UTC](https://discourse.pangeo.io/t/conservative-region-aggregation-with-xarray-geopandas-and-sparse/2715/12 "2022-09-08T12:50:31Z")
> [@benbovy](#):
>
> Also, naive question: do the polygon shapes distorted by the equal-area projection, combined with their limited number of vertices (grid and/or region geometries) and/or their extent, have any significant effect on the accuracy of the computed area weights?
This is a great question. I don’t know enough about geography and projections to answer it. In the notebook I did a check that the country areas are preserved. But I suppose the distortions in geometry could affect the weights in other ways.
> [@benbovy](#):
>
> If it is acceptable to perform conservative aggregation of arbitrary polygons at a given finite resolution, I think that something based on s2geometry (or H3) cells may offer great speed-ups.
So for this use case, you are suggesting we first go from lat-lon grid to an arbitrary super high-resolution S3 or H3 grid, then to regions? The intermediate step would introduce some resampling error, and would require a two-stage regridding. Why are you so sure this would be faster? I’d love to see an example of that.
---
### Author:  [@benbovy](https://discourse.pangeo.io/u/benbovy)
#### Post date: [September 8, 2022, 2:24pm UTC](https://discourse.pangeo.io/t/conservative-region-aggregation-with-xarray-geopandas-and-sparse/2715/13 "2022-09-08T14:24:27Z")
> So for this use case, you are suggesting we first go from lat-lon grid to an arbitrary super high-resolution S3 or H3 grid, then to regions? The intermediate step would introduce some resampling error, and would require a two-stage regridding. Why are you so sure this would be faster? I’d love to see an example of that.
Yeah I wrote that a bit too quickly and I’m not sure of anything actually 🙂
With S2 you can approximate any arbitrary region as a union of cells, i.e., `S2CellUnion` via `S2RegionCoverer`, with more or less precision (e.g., [examples here](https://s2geometry.io/devguide/examples/coverings)). [S2CellUnion](https://github.com/google/s2geometry/blob/191fbeef600f39ce799b177a8e56da908e841d27/src/s2/s2cell_union.h#L60) provides the API to compute the area of a cell union and compute the intersection between two cell unions, which I thought would allow some optimization vs. the same computation based on the exact polygon geometries (at the expense of some loss in precision, which can be controlled by parameters). I haven’t looked more into it, though, so I might be completely wrong!
---
### Author:  [@hansmohrmann](https://discourse.pangeo.io/u/hansmohrmann)
#### Post date: [September 9, 2022, 2:53pm UTC](https://discourse.pangeo.io/t/conservative-region-aggregation-with-xarray-geopandas-and-sparse/2715/14 "2022-09-09T14:53:20Z")
Thank you for sharing! the low-dependency calculation of grid weights for arbitrary geometries is very useful.
A comment on the final weight application (which I know isn’t the key part of this demo) - I haven’t been able to see much performance using sparse matrices (on my setup, it takes ~20s to compute the dot product with your custom ufunc). I have found that using dask.tensordot on dense arrays (with a bit more massaging) gives a faster solution (~3x faster in my case).
```auto
# at the top of your "Perform Matrix Multiplication" section
import dask
weight_da = ds_weights.unstack(fill_value=0).weights #no sparse
good_lats = np.isin(precip.latitude.values, weight_da.latitude.values)
precip = precip.isel(latitude=good_lats) # to match lat/lon dims of weights and precip data
precip_in_mem = precip.compute().chunk({"time": "10MB"})
agged_td = dask.array.tensordot(precip_in_mem.data, weight_da, axes=2).compute()
precip_regridded_td = xr.DataArray(data=agged_td, coords=[precip_in_mem['time'], weight_da['SOVEREIGNT']])
```
Might just be my particular setup, or a quirk of dask better optimizing tensordot vs. a generic ufunc.
---
### Author:  [@martinfleis](https://discourse.pangeo.io/u/martinfleis)
#### Post date: [November 5, 2022, 8:48pm UTC](https://discourse.pangeo.io/t/conservative-region-aggregation-with-xarray-geopandas-and-sparse/2715/15 "2022-11-05T20:48:33Z")
> [@benbovy](#):
>
> Also, naive question: do the polygon shapes distorted by the equal-area projection, combined with their limited number of vertices (grid and/or region geometries) and/or their extent, have any significant effect on the accuracy of the computed area weights?
This would be a bigger problem if you’d reproject the grid than country geometries. Larger grids defined by 4 planar points only tend to get messed up in retrojecting and you should densify then before doing that. You can check [pyproj’s solution for reprojecting bounding boxes](https://pyproj4.github.io/pyproj/dev/api/transformer.html#pyproj.transformer.Transformer.transform_bounds).
Also, I should note that a lot of work on interpolation from polygons to polygons has been done in the [Tobler package](https://pysal.org/tobler/api.html) and some of its internals may give a bit of speedup over using `overlay` that comes with a bit of overhead if you are interested only in the weights matrix.
---
### Author:  [@rsignell](https://discourse.pangeo.io/u/rsignell)
#### Post date: [November 8, 2022, 7:43pm UTC](https://discourse.pangeo.io/t/conservative-region-aggregation-with-xarray-geopandas-and-sparse/2715/16 "2022-11-08T19:43:37Z")
I slightly extended the notebook from @rabernat to show a map of mean precip by region (with a little help from @chegint and @ocefpaf): [https://nbviewer.org/gist/b02042191c94b9882af771f13d95693f](https://nbviewer.org/gist/b02042191c94b9882af771f13d95693f)
---
### Author:  [@rabernat](https://discourse.pangeo.io/u/rabernat)
#### Post date: [November 14, 2022, 5:52pm UTC](https://discourse.pangeo.io/t/conservative-region-aggregation-with-xarray-geopandas-and-sparse/2715/17 "2022-11-14T17:52:41Z")
Rich that is awesome! 🤩
---
### Author:  [@allixender](https://discourse.pangeo.io/u/allixender)
#### Post date: [December 6, 2022, 6:29am UTC](https://discourse.pangeo.io/t/conservative-region-aggregation-with-xarray-geopandas-and-sparse/2715/18 "2022-12-06T06:29:08Z")
@benbovy @rabernat Hey, just tagging here. With @annefou and @tinaok we would be keen to explore that more systematically from the DGGS \<-\> unstructured grids side (and myself in general region aggregation over “raster” or “vector” data with DGGS anyway). Just pasting the link to my thread here where we hope to pick up the slack a bit [Discrete Global Grid Systems (DGGS) use with Pangeo - #25 by allixender](https://discourse.pangeo.io/t/discrete-global-grid-systems-dggs-use-with-pangeo/2274/25)
---
### Author:  [@rsignell](https://discourse.pangeo.io/u/rsignell)
#### Post date: [January 18, 2023, 6:35pm UTC](https://discourse.pangeo.io/t/conservative-region-aggregation-with-xarray-geopandas-and-sparse/2715/19 "2023-01-18T18:35:51Z")
@rabernat, in your notebook you mention that it would be nice if this just worked:
```auto
regridded = xr.dot(ds[var], weights_sparse)
```
but pointed out that unfortunately einsum had not been implemented in sparse.
There is a PR here that looks promising: [einsum by jcmgray · Pull Request #564 · pydata/sparse · GitHub](https://github.com/pydata/sparse/pull/564) and the author is looking for testers.
I tried doing:
```python
regridded = sparse.einsum('ij,jk->ik', ds[var], weights_sparse)
```
but got back:
```plaintext
AttributeError: 'DataArray' object has no attribute 'reshape'
```
I don’t really know what I’m doing here, so I stopped and posted. 🙂
---
### Author:  [@TomNicholas](https://discourse.pangeo.io/u/TomNicholas)
#### Post date: [January 18, 2023, 8:17pm UTC](https://discourse.pangeo.io/t/conservative-region-aggregation-with-xarray-geopandas-and-sparse/2715/20 "2023-01-18T20:17:46Z")
@rsignell did you try doing
```python
regridded = xr.dot(ds[var], weights_sparse)
```
with that particular sparse PR installed? Then if your dataset is backed by a sparse array then `xr.dot` _should_ dispatch to the sparse version of einsum…
If it’s still not working at that point we should raise an xarray issue.
[Next page](https://discourse.pangeo.io/t/conservative-region-aggregation-with-xarray-geopandas-and-sparse/2715.md?page=2)