<jupyter_start><jupyter_text>Averaging over a regionAlthough this may not sound like _real_ regridding, averaging a gridded fieldover a region is supported by `ESMF`. This works because the `conservative`regridding method preserves the areal average of the input field. That is, _thevalue at each output grid cell is the average input value over the output gridarea_. Instead of mapping the input field unto rectangular outputs cells, it'smapped unto an irregular mesh defined by an outer polygon. In other words,applying the regridding weights computes the exact areal-average of the inputgrid over each polygon.This process relies on converting `shapely.Polygon` and `shapely.MultiPolygon`objects into `ESMF.Mesh` objects. However, ESMF meshes do not support allfeatures that come with shapely's (Multi)Polyons. Indeed, mesh elements do notsupport interior holes, or multiple non-touching parts, as do `shapely` objects.The `xesmf.SpatialAverager` class works around these issues by computingindependent weights for interior holes and multi-part geometries, beforecombining the weights.Transforming polygons into a `ESMF.Mesh` is a slow process. Users looking forfaster (but approximate) methods may want to explore[regionmask](https://regionmask.readthedocs.io/) or[clisops](https://clisops.readthedocs.io).The following example shows just how simple it is to compute the average overdifferent countries. The notebook used `geopandas`, a simple and efficientcontainer for geometries, and `descartes` for plotting maps. Make sure bothpackages are installed, as they are not `xesmf` dependencies.<jupyter_code>%matplotlib inline
import matplotlib.pyplot as plt
import geopandas as gpd
import pandas as pd
from shapely.geometry import Polygon, MultiPolygon
import numpy as np
import xarray as xr
import xesmf as xe
import warnings

warnings.filterwarnings("ignore")
xr.set_options(display_style='text')<jupyter_output><empty_output><jupyter_text>Simple exampleIn this example we'll create a synthetic global field, then compute its averageover six countries. Download country outlines<jupyter_code># Load some polygons from the internet
regs = gpd.read_file(
    "https://cdn.jsdelivr.net/npm/world-atlas@2/countries-10m.json"
)

# Select a few countries for the sake of the example
regs = regs.iloc[[5, 9, 37, 67, 98, 155]]

# Simplify the geometries to a 0.02 deg tolerance, which is 1/100 of our grid.
# The simpler the polygons, the faster the averaging, but we lose some precision.
regs["geometry"] = regs.simplify(tolerance=0.02, preserve_topology=True)
regs
# Create synthetic global data
ds = xe.util.grid_global(2, 2)
ds = ds.assign(field=xe.data.wave_smooth(ds.lon, ds.lat))
ds
# Display the global field and countries' outline.
fig, ax = plt.subplots()
ds.field.plot(ax=ax, x="lon", y="lat")
regs.plot(ax=ax, edgecolor="k", facecolor="none")<jupyter_output><empty_output><jupyter_text>Compute the field average over each country`xesmf.SpatialAverager` is a class designed to average an `xarray.DataArray`over a list of polygons. It behaves similarly to `xesmf.Regridder`, but hasoptions to deal specifically with polygon outputs. It uses the `conservative`regridding, and can store and reuse weights.<jupyter_code>savg = xe.SpatialAverager(ds, regs.geometry, geom_dim_name="country")
savg<jupyter_output><empty_output><jupyter_text>When called, the `SpatialAverager` instance returns a `DataArray` of averagesover the `geom` dimension, here countries. `lon` and `lat` coordinates are thecentroids each polygon.<jupyter_code>out = savg(ds.field)
out = out.assign_coords(country=xr.DataArray(regs["name"], dims=("country",)))
out<jupyter_output><empty_output><jupyter_text>As the order of the polygons is conserved in the output, we can easily includethe results back into our `geopandas` dataframe.<jupyter_code>regs["field_avg"] = out.values
regs
fig, ax = plt.subplots(figsize=(12, 5))
ds.field.plot(ax=ax, x="lon", y="lat", cmap="Greys_r")
handles = regs.plot(
    column="field_avg", ax=ax, edgecolor="k", vmin=1, vmax=3, cmap="viridis"
)<jupyter_output><empty_output><jupyter_text>Extract the weight mask from the averagerThe weights are stored in a sparse matrix structure in`SpatialAverager.weights`. The sparse matrix can be converted to a fullDataArray, but note that this will increase memory usage proportional to thenumber of polygons.<jupyter_code># Convert sparse matrix to numpy array, it has size : (n_in, n_out)
# So reshape to the same shape as ds + polygons
w = xr.DataArray(
    savg.weights.toarray().reshape(regs.geometry.size, *ds.lon.shape),
    dims=("country", *ds.lon.dims),
    coords=dict(country=out.country, **ds.lon.coords),
)
plt.subplots_adjust(top=0.9)
facets = w.plot(col="country", col_wrap=2, aspect=2, vmin=0, vmax=0.05)
facets.cbar.set_label("Averaging weights")<jupyter_output><empty_output><jupyter_text>This also allows to quickly check that the weights are indeed normalized, thatthe sum of each mask is 1.<jupyter_code>w.sum(dim=["y", "x"]).values<jupyter_output><empty_output>
StarcoderdataJupyterStructuredCleanDedup