Build utility providing general regridding support between arbitrary projections

[lpilz]

Hi there and hi again Huite!

Huite and I met at EGU24 (which was very nice) and already talked a bit about this usecase, so here is the discussion thread as promised.

I am working on a problem where I have to remap emissions conservatively between grids of two different geographical projections (from lat/lon to LCC projection). I think there are two solutions to this problem. The first one is a monte-carlo approach (which I haven't looked at yet but might in the future for curiosity's sake). And the second one an approach constructing polygons (with additional vertices on the edges in order to represent the curvilinear form) from those projections and intersecting them.

I already implemented the second approach using geopandas and shapely but unfortunately it's awfully slow. So I was wondering whether this could be solved more efficiently by using xugrid.

Here is the Setup: I have two xarray datasets which look similar to these:

import numpy as np
import xarray as xr
import pyproj
ds_wrf = xr.Dataset(data_vars={"E_CO2": (("Time",  "y", "x"), np.zeros((5,3,2)), {'grid_mapping': 'wrf_projection'}), "wrf_projection": xr.Variable(data=pyproj.CRS.from_string('+proj=lcc +x_0=0 +y_0=0 +a=6370000 +b=6370000 +lat_1=48.0 +lat_2=53.0 +lat_0=50.94319152832031 +lon_0=8.5 +type=crs'), dims=[])}, coords={"x": ("x", np.array([-1086097.805118, -1071097.805118])), "y": ("y", np.array([-1267043.06047202, -1252043.06047202, -1237043.06047202])), "XLAT": (("y", "x"), np.array([[38.79788 , 38.81968 ],[38.9285  , 38.950367],[39.05917 , 39.08109 ]])), "XLONG": (("y", "x"), np.array([[-3.8602905, -3.6926575],[-3.8885193, -3.72052  ],[-3.91687  , -3.7485046]]))})
ds_emiss = xr.Dataset(data_vars={"E_CO2":  (("Time", "lat", "lon"), np.random.rand(125).reshape(5,5,5)), "projection": xr.Variable(data=pyproj.CRS.from_epsg(4326), dims=[])}, coords={"lat": ("lat", np.linspace(38.6, 40, 5)), "lon": ("lon", np.linspace(-4, -3.6, 5))})

I would like to do something like this:

ds_wrf = polygonize(ds_wrf, num_additional_vertices=2)
ds_emiss = polygonize(ds_emiss, num_additional_vertices=2)
ds_emiss = change_projection(ds_emiss, target_proj=ds_wrf.wrf_projection.item())
ds_wrf.E_CO2.xugrid.intersect(ds_emiss.E_CO2)

Can somebody maybe tell me if this is a sensible approach at all? Does this work with what functionality xugrid is providing at the moment or would you have to add additional functionality? And if this is easy: Is there a way to save the interpolation matrix mapping between lat,lon to x,y as a sparse matrix somehow?

Thanks a lot!

Lukas

[Huite]

Hi @lpilz,

I was wondering when a question would pop up!
Did you take a look yet at the User Guide for regridding in general, and specifically the overlap regridder?

https://deltares.github.io/xugrid/examples/regridder_overview.html#regridding-overview
https://deltares.github.io/xugrid/examples/overlap_regridder.html

If you already have the original data in a geodataframe, with refined polygons you should be able to use it to make an xugrid dataset directly. You can the relevant variable into a UgridDataArray:

import xugrid

uds = xugrid.UgridDataset.from_geodataframe(gdf)
uda = uds["E_CO2"]

This data can be reprojected in a straightforward manner. This doesn't change the data, it simply reprojects all the vertices of the grid to the new coordinate system:

uda_reprojected = uda.ugrid.to_crs(epsg=...)  # can also be a pyproj crs object

You should also define your target topology in way that xugrid will understand it. I'm guessing this is a structured grid in your preferred coordinate system?
Xugrid will understand "ordinary" structured DataArrays with y and x coordinates. I'll call this da_target below.
Alternatively, you can make your own target topology as well. Valid target type are xugrid.UgridDataArray, xarray.DataArray, or xugrid.Ugrid2d.

To compute all the overlaps, we instantiate an OverlapRegridder:

regridder = xugrid.OverlapRegridder(source=uda_reprojected, target=da_target)

This will immediately run the intersection logic and compute all overlaps. These are stored in a sparse matrix, you can get the values with the .weights property:

overlap_matrix = regridder.weights

This'll return a Dataset with all the relevant data for the regridder, the fields basically match what you need for a scipy CSR matrix.

In general though, the idea is that the regridder can do the remapping for you, so you don't have to deal with these weights matrices.
There's a number of methods available for the regridding, in that case it should be quite straightforward:

regridder = xugrid.OverlapRegridder(source=uda_reprojected, target=da_target, method="...")
remapped_da = regridder.regrid(source=uda_reprojected)

There is machinery to do this in a delayed manner for orthogonal dimensions (e.g. a time dimensions, or a vertical dimension layer), so it'll process large datasets as well (provided they are properly chunked).

I'd be curious to hear what method you need and whether you can use one of the provided methods. See:

https://deltares.github.io/xugrid/api/xugrid.OverlapRegridder.html#xugrid.OverlapRegridder
https://deltares.github.io/xugrid/api/xugrid.RelativeOverlapRegridder.html#xugrid.RelativeOverlapRegridder

Note that it's also possible to provide your own custom reduction method, see the custom reductions part in the docs:
https://deltares.github.io/xugrid/examples/overlap_regridder.html#custom-reductions