Regridding overview#

Regridding is the process of converting gridded data from one grid to another grid. Xugrid provides tools for 2D and 3D regridding of structured gridded data, represented as xarray objects, as well as (layered) unstructured gridded data, represented as xugrid objects.

A number of regridding methods are provided, based on area or volume overlap, as well as interpolation routines. It currently only supports Cartesian coordinates. See e.g. xESMF instead for regridding with a spherical Earth representation (note: EMSF is not available via conda-forge on Windows).

Here are a number of quick examples of how to get started with regridding.

We’ll start by importing a few essential packages.

import matplotlib.pyplot as plt
import xarray as xr

import xugrid as xu

We will take a look at a sample dataset: a triangular grid with the surface elevation of the Netherlands.

uda = xu.data.elevation_nl()
uda.ugrid.plot(vmin=-20, vmax=90, cmap="terrain")
regridder overview
<matplotlib.collections.PolyCollection object at 0x7f748aeb0e90>

Xugrid provides several “regridder” classes which can convert gridded data from one grid to another grid. Let’s generate a very simple coarse mesh that covers the entire Netherlands.

def create_grid(bounds, nx, ny):
    """Create a simple grid of triangles covering a rectangle."""
    import numpy as np
    from matplotlib.tri import Triangulation

    xmin, ymin, xmax, ymax = bounds
    dx = (xmax - xmin) / nx
    dy = (ymax - ymin) / ny
    x = np.arange(xmin, xmax + dx, dx)
    y = np.arange(ymin, ymax + dy, dy)
    y, x = [a.ravel() for a in np.meshgrid(y, x, indexing="ij")]
    faces = Triangulation(x, y).triangles
    return xu.Ugrid2d(x, y, -1, faces)


grid = create_grid(uda.ugrid.total_bounds, 7, 7)

CentroidLocatorRegridder#

An easy way of regridding is by simply looking in which cell of the original the centroids of the new grid fall.

fig, ax = plt.subplots()
uda.ugrid.plot(vmin=-20, vmax=90, cmap="terrain", ax=ax)
grid.plot(ax=ax, color="red")
ax.scatter(*grid.centroids.T, color="red")
regridder overview
<matplotlib.collections.PathCollection object at 0x7f748b1f90a0>

Xugrid provides the CentroidLocatorRegridder for this:

regridder = xu.CentroidLocatorRegridder(source=uda, target=grid)
result = regridder.regrid(uda)
result.ugrid.plot(vmin=-20, vmax=90, cmap="terrain", edgecolor="red")
regridder overview
<matplotlib.collections.PolyCollection object at 0x7f748b2db830>

OverlapRegridder#

Such a regridding is not appropriate when the new grid cells are so large. Let’s try the OverlapOverregridder instead.

regridder = xu.OverlapRegridder(source=uda, target=grid)
mean = regridder.regrid(uda)
mean.ugrid.plot(vmin=-20, vmax=90, cmap="terrain", edgecolor="red")
regridder overview
<matplotlib.collections.PolyCollection object at 0x7f748c13e000>

By default, the OverlapRegridder computes an area weighted mean. Let’s try again, now with the minimum:

regridder = xu.OverlapRegridder(source=uda, target=grid, method="minimum")
minimum = regridder.regrid(uda)
minimum.ugrid.plot(vmin=-20, vmax=90, cmap="terrain", edgecolor="red")
regridder overview
<matplotlib.collections.PolyCollection object at 0x7f748ae67e60>

Or the maximum:

regridder = xu.OverlapRegridder(source=uda, target=grid, method="maximum")
maximum = regridder.regrid(uda)
maximum.ugrid.plot(vmin=-20, vmax=90, cmap="terrain", edgecolor="red")
regridder overview
<matplotlib.collections.PolyCollection object at 0x7f7495cf0e30>

All regridders also work for multi-dimensional data.

Let’s pretend our elevation dataset contains multiple layers, for example to denote multiple geological strata. We’ll generate five layers, each with a thickness of 10.0 meters.

thickness = xr.DataArray(
    data=[10.0, 10.0, 10.0, 10.0, 10.0],
    coords={"layer": [1, 2, 3, 4, 5]},
    dims=["layer"],
)

We need to make that the face dimension remains last, so we transpose the result.

bottom = (uda - thickness.cumsum("layer")).transpose()
bottom
<xarray.DataArray (layer: 5, mesh2d_nFaces: 5248)> Size: 210kB
array([[ -8.83000004,  -0.18999958,  44.04000092, ...,  -9.72      ,
        -25.82999992, -10.44999999],
       [-18.83000004, -10.18999958,  34.04000092, ..., -19.72      ,
        -35.82999992, -20.44999999],
       [-28.83000004, -20.18999958,  24.04000092, ..., -29.72      ,
        -45.82999992, -30.44999999],
       [-38.83000004, -30.18999958,  14.04000092, ..., -39.72      ,
        -55.82999992, -40.44999999],
       [-48.83000004, -40.18999958,   4.04000092, ..., -49.72      ,
        -65.82999992, -50.44999999]])
Coordinates:
    mesh2d_face_x  (mesh2d_nFaces) float64 42kB ...
    mesh2d_face_y  (mesh2d_nFaces) float64 42kB ...
  * layer          (layer) int64 40B 1 2 3 4 5
  * mesh2d_nFaces  (mesh2d_nFaces) int64 42kB 0 1 2 3 4 ... 5244 5245 5246 5247


We can feed the result to the regridder, which will automatically regrid over all additional dimensions.

mean_bottom = xu.OverlapRegridder(source=bottom, target=grid).regrid(bottom)
mean_bottom
<xarray.DataArray (layer: 5, mesh2d_nFaces: 98)> Size: 4kB
array([[ 98.73378481,  24.75605825,          nan,          nan,
                 nan,          nan,  28.6866454 ,  21.59076039,
                 nan,          nan, -10.30473318, -12.46283808,
                 nan,          nan,   1.98885124,  -0.45315257,
         -7.76052909,  -7.60938601, -13.9746031 , -12.82278477,
        -11.86447095, -13.9699236 ,          nan,          nan,
        -11.32250825, -12.61944108,          nan,          nan,
        -10.39512475,  -8.94184232,   7.71162015,   4.05908219,
         11.63407438,  21.07322886,   8.79753456,  16.4475533 ,
         16.63484974,  15.75826454,  -6.61134759,  -5.75454656,
          7.20785073,  -0.81124398,          nan,          nan,
         69.34102272,          nan,          nan,          nan,
          8.73683325,  17.46905942,          nan, -11.27368816,
                 nan,          nan,          nan,          nan,
                 nan,          nan,  -8.92800775, -11.78767776,
         -4.59914228, -11.1621325 ,          nan,          nan,
                 nan,          nan,  20.30399372,  14.92605801,
          4.65570035,  -0.30620094,  -9.73834196, -10.55580261,
        -11.23052693, -11.40449714, -12.85176061, -11.33053948,
         -7.89456194, -10.43064765, -11.51684987, -10.9563055 ,
...
        -51.86447095, -53.9699236 ,          nan,          nan,
        -51.32250825, -52.61944108,          nan,          nan,
        -50.39512475, -48.94184232, -32.28837985, -35.94091781,
        -28.36592562, -18.92677114, -31.20246544, -23.5524467 ,
        -23.36515026, -24.24173546, -46.61134759, -45.75454656,
        -32.79214927, -40.81124398,          nan,          nan,
         29.34102272,          nan,          nan,          nan,
        -31.26316675, -22.53094058,          nan, -51.27368816,
                 nan,          nan,          nan,          nan,
                 nan,          nan, -48.92800775, -51.78767776,
        -44.59914228, -51.1621325 ,          nan,          nan,
                 nan,          nan, -19.69600628, -25.07394199,
        -35.34429965, -40.30620094, -49.73834196, -50.55580261,
        -51.23052693, -51.40449714, -52.85176061, -51.33053948,
        -47.89456194, -50.43064765, -51.51684987, -50.9563055 ,
        -50.24352139, -49.97530487, -54.90399917, -50.87557562,
        -50.05098298, -50.91804551, -39.44818058, -44.02645019,
        -34.95904013, -31.75848616, -53.71649682, -47.7613762 ,
        -46.45744354, -42.33120932, -51.24098772, -50.25680056,
        -45.92794405, -39.50867478]])
Coordinates:
  * layer          (layer) int64 40B 1 2 3 4 5
  * mesh2d_nFaces  (mesh2d_nFaces) int64 784B 0 1 2 3 4 5 ... 92 93 94 95 96 97


Let’s take a slice to briefly inspect our original layer bottom elevation, and the aggregated mean.

section_y = 475_000.0
section = bottom.ugrid.sel(y=section_y)
section_mean = mean_bottom.ugrid.sel(y=section_y)

fig, (ax0, ax1) = plt.subplots(ncols=2, figsize=(10, 5), sharex=True, sharey=True)
section.plot.line(x="mesh2d_s", hue="layer", ax=ax0)
section_mean.plot.line(x="mesh2d_s", hue="layer", ax=ax1)
regridder overview
[<matplotlib.lines.Line2D object at 0x7f748c506bd0>, <matplotlib.lines.Line2D object at 0x7f748c506c60>, <matplotlib.lines.Line2D object at 0x7f748bea8c50>, <matplotlib.lines.Line2D object at 0x7f7494af8680>, <matplotlib.lines.Line2D object at 0x7f748bea8e90>]

BarycentricInterpolator#

All examples above show reductions: from a fine grid to a coarse grid. However, xugrid also provides interpolation to generate smooth fine representations of a coarse grid.

To illustrate, we will zoom in to a part of the Netherlands.

part = uda.ugrid.sel(x=slice(125_000, 225_000), y=slice(440_000, 500_000))
part.ugrid.plot(vmin=-20, vmax=90, cmap="terrain")
regridder overview
<matplotlib.collections.PolyCollection object at 0x7f7494a18350>

We can clearly identify the individual triangles that form the grid. To get a smooth presentation, we can use the BarycentricInterpolator.

We will generate a fine grid.

grid = create_grid(part.ugrid.total_bounds, nx=100, ny=100)

We use the centroids of the fine grid to interpolate between the centroids of the triangles.

regridder = xu.BarycentricInterpolator(part, grid)
interpolated = regridder.regrid(part)
interpolated.ugrid.plot(vmin=-20, vmax=90, cmap="terrain")
regridder overview
<matplotlib.collections.PolyCollection object at 0x7f748c304da0>

Arbitrary grids#

The above examples all feature triangular source and target grids. However, the regridders work for any collection of (convex) faces.

grid = create_grid(part.ugrid.total_bounds, nx=20, ny=15)
voronoi_grid = grid.tesselate_centroidal_voronoi()

regridder = xu.CentroidLocatorRegridder(part, voronoi_grid)
result = regridder.regrid(part)

fig, ax = plt.subplots()
result.ugrid.plot(vmin=-20, vmax=90, cmap="terrain")
voronoi_grid.plot(ax=ax, color="red")
regridder overview
<matplotlib.collections.LineCollection object at 0x7f748b43b260>

Re-use#

The most expensive step of the regridding process is finding and computing overlaps. A regridder can be used repeatedly, provided the source topology is kept the same.

part_other = part - 50.0
result = regridder.regrid(part_other)
result.ugrid.plot(vmin=-20, vmax=90, cmap="terrain")
regridder overview
<matplotlib.collections.PolyCollection object at 0x7f748b0f8e90>

Total running time of the script: (0 minutes 9.365 seconds)

Gallery generated by Sphinx-Gallery