Regridding

regrid() maps coefficients in one space to another along an axis of an array; regrid_nd() does several axes, and rebasis.xarray.regrid() works on xarray objects by dimension name. Methods:

• "project" (default): the \(L^2\) projection, conservative.

• "measure": the target elements are averaging kernels.

• "interpolate": point values at the target centers.

• "fit": the source elements are the kernels that measured the data (e.g. bin averages), the target a model space: smooth conservative upsampling.

import matplotlib.pyplot as plt
import numpy as np

import rebasis as rb
from _plot import BLUE, MUTED, ORANGE

Bounded projection

\(L^2\) projection onto a coarse smooth basis rings near discontinuities. bounded=True finds the best approximation that conserves the integral and creates no new extrema (each coefficient within the range of the inputs it overlaps); bounded="positive" and bounded=(lo, hi) apply global bounds.

x = np.linspace(0, 1, 4001)
f = ((x > 0.3) & (x < 0.55)).astype(float)
src = rb.Space.hats(x)
dst = rb.Space.bsplines_clamped(np.linspace(0, 1, 16), 3)
plain = rb.regrid(f, src, dst)
bounded = rb.regrid(f, src, dst, bounded=True)

xx = np.linspace(0, 0.9999, 1000)
E = dst.evaluate(xx)
m = dst.integrals()
fig, ax = plt.subplots(figsize=(7, 3))
ax.plot(x, f, color=MUTED, lw=1, label=f"input (∫ = {f @ src.integrals():.6f})")
ax.plot(xx, E @ plain, color=BLUE, label=f"L2 projection (∫ = {plain @ m:.6f})")
ax.plot(xx, E @ bounded, color=ORANGE, label=f"bounded (∫ = {bounded @ m:.6f})")
ax.legend(loc="upper right");

Smooth conservative upsampling

From bin averages to a spline whose bin averages reproduce them:

edges = np.linspace(0, 2 * np.pi, 9)
avg = np.diff(-np.cos(edges)) / np.diff(edges)  # bin averages of sin
bins = rb.Space.boxes(edges)
model = rb.Space.bsplines_clamped(edges, 3)
coef = rb.regrid(avg, bins, model, method="fit", smoothing=1e-8)

xx = np.linspace(0, 2 * np.pi * 0.9999, 500)
fig, ax = plt.subplots(figsize=(7, 3))
ax.stairs(avg, edges, color=MUTED, label="bin averages")
ax.plot(xx, model.evaluate(xx) @ coef, color=BLUE, label="spline with the same bin averages")
ax.legend(loc="lower left")
print("bin averages reproduced to", np.abs(rb.measure(bins, model) @ coef - avg).max())

bin averages reproduced to 4.289765209719576e-10

Area-weighted latitude–longitude regridding

Grids that are tensor products are regridded one axis at a time, which is exact for separable weights. "cos_lat_deg" weights latitude by area:

lat_f, lon_f = rb.Space.boxes(np.linspace(-90, 90, 181)), rb.Space.boxes(np.linspace(0, 360, 361))
lat_c, lon_c = rb.Space.boxes(np.linspace(-90, 90, 13)), rb.Space.boxes(np.linspace(0, 360, 25))
lat, lon = np.meshgrid(lat_f.centers, lon_f.centers, indexing="ij")
field = np.cos(np.radians(lat)) ** 2 * (1 + 0.3 * np.sin(np.radians(3 * lon)))

coarse = rb.regrid_nd(field, {0: (lat_f, lat_c, "cos_lat_deg"), 1: (lon_f, lon_c)})


def global_mean(v, lat_s, lon_s):
    w = np.outer(lat_s.integrals("cos_lat_deg"), lon_s.integrals())
    return (v * w).sum() / w.sum()


print(global_mean(field, lat_f, lon_f), global_mean(coarse, lat_c, lon_c))

0.6666582044272394 0.6666582044272394

xarray

import xarray as xr

import rebasis.xarray

da = xr.DataArray(field, dims=("lat", "lon"), coords={"lat": lat_f.centers, "lon": lon_f.centers})
rebasis.xarray.regrid(
    da,
    {"lat": np.arange(-87.5, 90, 5.0), "lon": np.arange(2.5, 360, 5.0)},
    weight={"lat": "cos_lat_deg"},
)

<xarray.DataArray (lat: 36, lon: 72)> Size: 21kB
array([[0.00387015, 0.00415114, 0.00440307, ..., 0.00304633, 0.00329826,
        0.00357925],
       [0.0195573 , 0.02097725, 0.02225034, ..., 0.01539421, 0.0166673 ,
        0.01808725],
       [0.05045494, 0.0541182 , 0.0574026 , ..., 0.03971479, 0.04299918,
        0.04666245],
       ...,
       [0.05045494, 0.0541182 , 0.0574026 , ..., 0.03971479, 0.04299918,
        0.04666245],
       [0.0195573 , 0.02097725, 0.02225034, ..., 0.01539421, 0.0166673 ,
        0.01808725],
       [0.00387015, 0.00415114, 0.00440307, ..., 0.00304633, 0.00329826,
        0.00357925]], shape=(36, 72))
Coordinates:
  * lat      (lat) float64 288B -87.5 -82.5 -77.5 -72.5 ... 72.5 77.5 82.5 87.5
  * lon      (lon) float64 576B 2.5 7.5 12.5 17.5 ... 342.5 347.5 352.5 357.5

xarray.DataArray

• lat: 36
• lon: 72

• 0.00387 0.004151 0.004403 0.004609 ... 0.003046 0.003298 0.003579

  array([[0.00387015, 0.00415114, 0.00440307, ..., 0.00304633, 0.00329826,
          0.00357925],
         [0.0195573 , 0.02097725, 0.02225034, ..., 0.01539421, 0.0166673 ,
          0.01808725],
         [0.05045494, 0.0541182 , 0.0574026 , ..., 0.03971479, 0.04299918,
          0.04666245],
         ...,
         [0.05045494, 0.0541182 , 0.0574026 , ..., 0.03971479, 0.04299918,
          0.04666245],
         [0.0195573 , 0.02097725, 0.02225034, ..., 0.01539421, 0.0166673 ,
          0.01808725],
         [0.00387015, 0.00415114, 0.00440307, ..., 0.00304633, 0.00329826,
          0.00357925]], shape=(36, 72))

• Coordinates: (2)
  • lat
    (lat)
    float64
    -87.5 -82.5 -77.5 ... 82.5 87.5

    array([-87.5, -82.5, -77.5, -72.5, -67.5, -62.5, -57.5, -52.5, -47.5, -42.5,
           -37.5, -32.5, -27.5, -22.5, -17.5, -12.5,  -7.5,  -2.5,   2.5,   7.5,
            12.5,  17.5,  22.5,  27.5,  32.5,  37.5,  42.5,  47.5,  52.5,  57.5,
            62.5,  67.5,  72.5,  77.5,  82.5,  87.5])

  • lon
    (lon)
    float64
    2.5 7.5 12.5 ... 347.5 352.5 357.5

    array([  2.5,   7.5,  12.5,  17.5,  22.5,  27.5,  32.5,  37.5,  42.5,  47.5,
            52.5,  57.5,  62.5,  67.5,  72.5,  77.5,  82.5,  87.5,  92.5,  97.5,
           102.5, 107.5, 112.5, 117.5, 122.5, 127.5, 132.5, 137.5, 142.5, 147.5,
           152.5, 157.5, 162.5, 167.5, 172.5, 177.5, 182.5, 187.5, 192.5, 197.5,
           202.5, 207.5, 212.5, 217.5, 222.5, 227.5, 232.5, 237.5, 242.5, 247.5,
           252.5, 257.5, 262.5, 267.5, 272.5, 277.5, 282.5, 287.5, 292.5, 297.5,
           302.5, 307.5, 312.5, 317.5, 322.5, 327.5, 332.5, 337.5, 342.5, 347.5,
           352.5, 357.5])

Missing data: skipna=True regrids with NaNs set to zero and divides by the regridded validity mask, so outputs are averages of the valid inputs (na_thres sets the minimum valid fraction).