Add Mesh support and a spatial averaging method

[aulemahal]

Averaging some gridded data over a polygon is something that often comes up in our work at Ouranosinc (and elsewhere for sure). As fair as I know, there is no easy way to do this with xarray datasets. There is clisops's subset_shape (that was in xclim.subset before), but it creates a boolean mask with Trues when the grid centroid falls in the polygon. What we want is a weights map, the weights being computed as the fraction of each gridcell that intersects the polygon. That is exactly what ESMF Regridding does with the conservative method.

I suggest to read in a sequence of shapely.geometry.Polygon objects and generate a ESMF.Mesh from them. Using these meshes with a few elements of many nodes, we can call regridding with a conservative method to get the weights. I have working code with ESMpy that does that and would be pleased to make PR if the idea is accepted.

Also, bonus of this issue, we could accept Geopandas Dataframe as input and output to the regridding object. That would add a lot of dependencies though. Shapely alone is more lightweight.

I'm guessing that my use case would require a new method. There are some issues with this:

• "MultiPolygon" objects would be splitted to their components. Weights could be summed afterwards to keep the hierarchy.
• Polygon with holes require a 2-pass treatment where weights computed from the holes are to be subtracted to the main "exterior-ring" weights. Thus, I would suggest a xesmf.intersection_weights() method (or any other name).

[dcherian]

regionmask helps you do this when combined with DataArray.weighted: https://regionmask.readthedocs.io/en/stable/notebooks/mask_2D.html#calculate-weighted-regional-average