Overview of regridding/resampling

[Manangka]

In GitLab by @Huite on Sep 18, 2020, 12:56

@visr @DirkEilander This is a short overview of resampling methods, related to projected and unprojected geospatial grids and meshes.

GDAL (and GDAL via rasterio) are probably the most commonly known. Supports a variety of methods, and is suitable for reprojecting from one coordinate system to another as well. Limitations are primarily:

• Exclusively 2D, no convenience for layers, time. Reprojecting becomes rather costly, since it recomputes weights for every slice.
• Resampling methods are purely based on cell centers, no true evaluation of areal overlap
• Alll C++ in a massive binary, so there’s no easy way to add different resampling methods.

Because I wanted a number of different methods for resampling (e.g. harmonic mean, conservative methods for conductance), I started looking around some time ago. First usable thing I stumbled upon is this library: https://github.com/esa-esdl/gridtools

I added some methods, but it doesn’t take long to realize that almost many methods are a reduction over a collection of weights and values. This can be abstracted, which is what I did for the imod Regridder: https://imod.xyz/api/prepare.html#imod.prepare.Regridder

As well as adding 1D and 3D methods. Since Numba supports JIT, it’s easy to add custom functions (provided they are a reduction over weights and values, linear interpolation doesn’t work this way, and requires a different implementation).

In terms of other numba based interpolation/resampling: https://github.com/EconForge/interpolation.py

(relies on string based code generation -> JIT, which is quirky, but no macros in Python...)

Just found this also while googling: https://github.com/dbstein/fast_interp

Pyresample is 2D, supports a number of methods, mostly focused on satellite imagery, so I think it can do reprojection: https://pyresample.readthedocs.io/en/latest/

Xemsf uses the python bindings to EMSF to do resampling (via regridding weighst netcdf): https://xesmf.readthedocs.io/en/latest/why.html Currently, it only supports rectilinear and curvilinear grids (since that’s what xarray supports). EMSF can do a lot more, but is a very complicated fortran binary (1 million LOC apparently), which isn’t available on Windows.

A bunch of other utilities:

• https://github.com/pysal/tobler: Tobler is a python package for areal interpolation, dasymetric mapping, and change of support.
• https://github.com/corteva/geocube: Tool to convert geopandas vector data into rasterized xarray data.
• https://github.com/corteva/rioxarray: “You can learn how to clip, merge, and reproject rasters in the Usage Examples section of the documentation. Need to export to a raster (GeoTiff)? There is an example for that as well.”
• https://github.com/makepath/xarray-spatial: Xarray-Spatial implements common raster analysis functions using Numba and provides an easy-to-install, easy-to-extend codebase for raster analysis.
• https://github.com/CNES/pangeo-pyinterp: Python library for optimized geo-referenced interpolation.

I’ve made a start with what a Regridder for unstructured 2D meshes looks like: https://github.com/Huite/xugrid/blob/master/xugrid/regrid.py

In contrast to e.g. rectilinear regridding, computing weights is a lot more expensive, so taking the EMSF route and storing the weights seems like a good idea.

I think I’d like to refactor the regridding stuff into something separate, easy-installable, which deals with both structured and unstructured data (since one often wants to convert one to the other). For performance, you’d be able to specialize on structured input; for unstructured meshes, you need some kind of spatial indexing scheme. See: https://github.com/ESM-VFC/xoak/issues/12 (I left a comment there as well).

In terms of unstructured data (2D unstructured, layer for z/depth), the UGRID convention is good and broadly applicable (it supports data on faces, edges, nodes). All mesh definitions I've encountered are very similar: coordinate arrays for the nodes, and index arrays mapping to those nodes.

The closest thing to an "in-memory" UGRID netCDF (like xarray is an "in-memory" netCDF) is the Gridded Dataset: https://github.com/NOAA-ORR-ERD/gridded

I'd like to integrate xarray with it (haven't started yet), but there's also some work coming where xarray will support more flexible indexing for unstructured meshes (see the discussion on the gridded issues).

[Manangka]

In GitLab by @DirkEilander on Sep 21, 2020, 09:04

Thanks for this extensive overview @Huite! I like your idea of making a separate package out of the regridding utilities in imod-python. I'll also have a look at the other references in your description.

[Manangka]

In GitLab by @DirkEilander on Sep 21, 2020, 09:04

mentioned in issue deltares/wflow/hydromt#33

[Manangka]

In GitLab by @DirkEilander on Sep 23, 2020, 17:48

In addition to the methods mentioned above:

• CDO provides several regridding methods including the First-order Conservative Remapping (Jones, P.W. 1999, Monthly Weather Review, 127, 2204-2210) which is the standardized ISMIP6 method. I don't have experience using this package, but there are also python bindings for CDO

• Iris has a regridding method, but requires your data to be presented in an Iris cube format. Xarray has to and from_iris methods to support this.

@Huite Is the CDO First-order Conservative Remapping similar to your conservative method for conductance?

@Huite Are you aware of any library that supports regridding from a geographical (e.g. WGS84) to a projected crs (e.g. UTM) like in gdal_translate/rasterio for 2D data?

[Manangka]

In GitLab by @Huite on Oct 12, 2020, 15:39

@DirkEilander Hmm the CDO first-order might be comparable indeed, although I guess it works for triangular meshes as well?

The issue with the Iris functions is that they use the Scipy interpolation methods, which require a lot of memory and rather slow for large grids in my experience.

If you're talking about 2D data, I take it you mean arbitrary meshes of (convex) cells. If we're talking unstructured meshes, I'd argue reprojection becomes much simpler. Such a mesh typical consists of points (nodes), and the defined by a face-to-node connectivity array. So you simply reproject all the nodes. If you're looking to reproject to a specific mesh, I'd do it in two steps: first reproject the nodes, and then just do cartesian regridding to get to the final destination grid (mind, I still need to write that regridded for 2D meshes...).

For advanced methods, we ought to have a look at: https://github.com/ClimateGlobalChange/tempestremap (There's two papers mentioned in the readme)

And the EMSF people have multiple conservative schemes, also higher order (with publications) if I recall correctly.

[Manangka]

In GitLab by @JoerivanEngelen on May 16, 2023, 18:04

Update: Unstructured regridding is implemented in the xugrid package. Overview unstructured regridding

When this PR for structured regridding is merged, the discussed functionality above is achieved: Regridding for structured and unstructured grids in a separate, more lightweight package. This can also be used in HydroMT.

The plan is to deprecate and eventually replace iMOD Python's regridder module with xugrids regridding: https://gitlab.com/deltares/imod/imod-python/-/issues/328

I'm closing this issue, as the work for this is pretty much done and I'm cleaning up the issue board.

Thanks for the extensive writeup, we can always link colleagues to this issue to show the regridding functionality has been thoroughly investigated.