Resampling signals on the 2-sphere#

This notebook demonstrates how to resample spherical signals between different grid resolutions and types using ResampleS2. The module performs separable interpolation along latitudes and longitudes, with proper handling of periodicity and pole values.

Two interpolation modes are available:

bilinear (default): standard linear interpolation in value space, \(f(t) = (1-t)\,f_0 + t\,f_1\)
bilinear-spherical: for fields whose values are angles in radians — wind directions, phases. It interpolates along the shorter arc of the circle, which is identical to bilinear unless two neighbouring samples straddle a \(2\pi\) wrap

Setup#

import math

import torch

from torch_harmonics import ResampleS2
from torch_harmonics.quadrature import precompute_latitudes
from torch_harmonics.utils import load_mola_elevation

import matplotlib.pyplot as plt
from torch_harmonics.plotting import plot_sphere

Test data#

We use the MOLA Mars topography as test data and set up an equiangular grid:

nlat = 257
nlon = 2*(nlat-1)
grid = "equiangular"

xq, wq = precompute_latitudes(nlat, grid=grid)

data = load_mola_elevation(nlat, nlon)

Let us visualize the original data:

plot_sphere(data, projection="orthographic", cmap="turbo", colorbar=True, title="Mars topography MOLA")
plt.show()

../_images/b157aef2608e35cfd077d3382cb66727a2197f8de0996c61985d12a10a899431.png

Grid-to-grid transfer#

ResampleS2 can transfer data between different grid types at the same resolution. Here we resample from equiangular to equiangular as a sanity check:

grid_up = "equiangular"
grid_down = "legendre-gauss"

input_resample = ResampleS2(nlat, nlon, nlat, nlon, grid_in="equiangular", grid_out=grid_up)

data_res = input_resample(data)

plot_sphere(data_res, projection="orthographic", cmap="turbo", colorbar=True, title="interpolated input data")

<cartopy.mpl.geocollection.GeoQuadMesh at 0x7f83789cc080>

../_images/49a0aab25669879c0fc6f490ee61c03981659dd6752debc56ceab14cc127cf27.png

Down- and upsampling#

We define downsampling and upsampling layers with a scale factor of 4. The downsampled grid uses Legendre–Gauss nodes, which is a common choice for spectral methods:

scale_factor = 4

nlat_down = math.ceil(nlat / scale_factor)
nlon_down = math.ceil(nlon / scale_factor)

downsample = ResampleS2(nlat, nlon, nlat_down, nlon_down, grid_in=grid_up)
upsample = ResampleS2(nlat_down, nlon_down, nlat, nlon, grid_in=grid_down)

data_down = downsample(data)
data_up = upsample(data_down)

The four panels show the original, downsampled, upsampled, and residual fields. The residual reveals information lost during the downsampling step:

fig = plt.figure(layout='constrained', figsize=(32, 6))
subfigs = fig.subfigures(1, 4)

# spectral fitting
plot_sphere(data, fig=subfigs[0], projection="orthographic", cmap="turbo", colorbar=True, title="Original")
plot_sphere(data_down, fig=subfigs[1], projection="orthographic", cmap="turbo", colorbar=True, title="down-sampled")
plot_sphere(data_up, fig=subfigs[2], projection="orthographic", cmap="turbo", colorbar=True, title="up-sampled")
plot_sphere((data_up - data), fig=subfigs[3], projection="orthographic", cmap="turbo", colorbar=True, title="residual (after upsampling)")
plt.show()

../_images/4061069d3ebdd8f83bf1c2197767eb16cc626c2083588766958e03068435eb1c.png

print((data_up - data).abs().max())

tensor(0.5244)

Comparison with SHT upsampling#

An alternative to interpolation-based upsampling is to go through spectral space: apply the forward SHT on the coarse grid, then the inverse SHT on the fine grid. This is spectrally exact for band-limited signals but can produce Gibbs-type ringing near sharp features:

from torch_harmonics import RealSHT, InverseRealSHT

sht = RealSHT(nlat_down, nlon_down, grid=grid_down)
isht = InverseRealSHT(nlat, nlon, lmax=sht.lmax, mmax=sht.mmax, grid=grid_up)

data_up = upsample(data_down)
data_sht_up = isht(sht(data_down))

fig = plt.figure(layout='constrained', figsize=(32, 12))
subfigs = fig.subfigures(2, 4)

# spectral fitting
plot_sphere(data, fig=subfigs[0, 0], projection="orthographic", cmap="turbo", colorbar=True, title="Original")
plot_sphere(data_down, fig=subfigs[0, 1], projection="orthographic", cmap="turbo", colorbar=True, title="down-sampled")
plot_sphere(data_up, fig=subfigs[0, 2], projection="orthographic", cmap="turbo", colorbar=True, title="up-sampled")
plot_sphere(data_sht_up, fig=subfigs[0, 3], projection="orthographic", cmap="turbo", colorbar=True, title="SHT up-sampled")
plot_sphere((data_up - data), fig=subfigs[1, 2], projection="orthographic", cmap="turbo", colorbar=True, title="residual (after upsampling)")
plot_sphere((data_sht_up - data), fig=subfigs[1, 3], projection="orthographic", cmap="turbo", colorbar=True, title="residual (after SHT upsampling)")
plt.show()

../_images/cf5054cf9f97a1648dd0cc49ee30870a32270c64f2d1b622b8ed711617e3b697.png

The SHT upsampling produces visible ringing near sharp topographic features (craters, ridges), while the interpolation-based approach stays smooth. The max errors confirm this:

print((data_up - data).abs().max())
print((data_sht_up - data).abs().max())

tensor(0.5244)
tensor(0.6142)

Bilinear vs. spherical interpolation#

ResampleS2 supports a second mode, mode="bilinear-spherical", for fields whose values are angles — wind directions, phases — measured in radians. Rather than blending in value space, it moves a fraction \(t\) along the shorter arc of the circle between the two samples:

\[ f(t) = f_0 + t \cdot \mathrm{wrap}_\pi(f_1 - f_0), \qquad \mathrm{wrap}_\pi(\delta) \in (-\pi, \pi] \]

This is the two-dimensional closed form of spherical linear interpolation (slerp). Since \(\mathrm{wrap}_\pi\) is the identity whenever \(|f_1 - f_0| \le \pi\), the two modes agree exactly on any field whose neighbouring samples differ by less than half a turn. They part company only across a \(2\pi\) wrap: halfway between \(f_0 = 3\) and \(f_1 = -3\) — nearly the same direction, either side of the branch cut — the spherical mode gives \(\pi\), while bilinear gives \(0\), the opposite direction.

The MOLA elevation used here is a smooth scalar field normalised to \([0, 1]\), so no pair of neighbours comes anywhere near \(\pi\). We upsample the same coarse data with both kernels to see what that implies:

upsample_slerp = ResampleS2(nlat_down, nlon_down, nlat, nlon, grid_in=grid_down, mode="bilinear-spherical")

# upsample the SAME coarse data with both kernels
data_up_bilinear = upsample(data_down)
data_up_slerp = upsample_slerp(data_down)

fig = plt.figure(layout='constrained', figsize=(32, 12))
subfigs = fig.subfigures(2, 3)

plot_sphere(data_up_bilinear, fig=subfigs[0, 0], projection="orthographic", cmap="turbo",
            colorbar=True, title="bilinear")
plot_sphere(data_up_slerp, fig=subfigs[0, 1], projection="orthographic", cmap="turbo",
            colorbar=True, title="bilinear-spherical")
plot_sphere((data_up_slerp - data_up_bilinear), fig=subfigs[0, 2], projection="orthographic", cmap="RdBu_r",
            colorbar=True, title="slerp $-$ lerp")

plot_sphere((data_up_bilinear - data).abs(), fig=subfigs[1, 0], projection="orthographic", cmap="turbo",
            colorbar=True, title="bilinear vs original")
plot_sphere((data_up_slerp - data).abs(), fig=subfigs[1, 1], projection="orthographic", cmap="turbo",
            colorbar=True, title="slerp vs original")
plot_sphere((data_up_slerp - data_up_bilinear).abs(), fig=subfigs[1, 2], projection="orthographic", cmap="turbo",
            colorbar=True, title="|slerp $-$ lerp|")
plt.show()

../_images/8b4d908bb9a41c3ee4e7bc092b114af4d0f0f4f13d836dc227f6beea2ec0801e.png

print(f"bilinear vs original:   {(data_up_bilinear - data).abs().max():.4f}")
print(f"slerp vs original:      {(data_up_slerp - data).abs().max():.4f}")
print(f"slerp vs bilinear:      {(data_up_slerp - data_up_bilinear).abs().max():.6f}")

bilinear vs original:   0.5244
slerp vs original:      0.5244
slerp vs bilinear:      0.000004

The two upsampled fields agree to floating-point rounding — around \(4\times10^{-6}\) here, against an aliasing error of \(0.52\) from the downsampling itself. They are not bit-identical only because the two modes evaluate algebraically equal expressions in a different order.

This is by construction rather than luck: elevation is not an angle and nothing in it wraps, so the shorter arc between two samples is just the straight line between them, and the spherical mode reduces to bilinear.

Applying bilinear-spherical to a non-angular field is therefore harmless but pointless — it either matches bilinear, or, if the field happens to span more than \(2\pi\), silently reinterprets its values as angles.

A case where the mode matters: a direction field#

Elevation never wraps, so the two kernels agree. To see where bilinear-spherical earns its keep, the field’s values have to be angles.

Below is a unit vector field \((u, v) = (\cos\varphi, \sin\varphi)\) — think of a wind direction that veers steadily with latitude, as in an Ekman spiral. The direction \(\varphi\) is constant along each latitude ring, so it stays single-valued at the poles, but it winds through \(\pm\pi\) twice between them. Each crossing is a branch cut: two neighbouring samples land either side of it, numerically a full turn apart while pointing in almost the same direction.

turns = 2.0


def wrap(x):
    """Fold an angle back into (-pi, pi]."""
    return torch.remainder(x + math.pi, 2 * math.pi) - math.pi


def wind_direction(nlat_, nlon_, grid_):
    """Direction that veers with latitude, constant along each ring."""
    colat, _ = precompute_latitudes(nlat_, grid=grid_)
    lat = math.pi / 2 - colat
    return wrap(2 * turns * lat).unsqueeze(-1).expand(nlat_, nlon_).contiguous()


dir_coarse = wind_direction(nlat_down, nlon_down, grid_down)  # what we upsample from
dir_ref = wind_direction(nlat, nlon, grid_up)  # the exact answer

dir_up_bilinear = upsample(dir_coarse)
dir_up_slerp = upsample_slerp(dir_coarse)

# directions are only defined modulo 2*pi, so the error has to be measured that way
err_bilinear = wrap(dir_up_bilinear - dir_ref).abs()
err_slerp = wrap(dir_up_slerp - dir_ref).abs()

fig = plt.figure(layout='constrained', figsize=(32, 12))
subfigs = fig.subfigures(2, 3)

# twilight is cyclic, so +pi and -pi share a colour, as they should for a direction
plot_sphere(dir_ref, fig=subfigs[0, 0], projection="orthographic", cmap="twilight",
            colorbar=True, title="direction (reference)")
plot_sphere(dir_up_bilinear, fig=subfigs[0, 1], projection="orthographic", cmap="twilight",
            colorbar=True, title="upsampled: bilinear")
plot_sphere(dir_up_slerp, fig=subfigs[0, 2], projection="orthographic", cmap="twilight",
            colorbar=True, title="upsampled: bilinear-spherical")

plot_sphere(wrap(dir_up_slerp - dir_up_bilinear).abs(), fig=subfigs[1, 0], projection="orthographic",
            cmap="magma", colorbar=True, title="|slerp $-$ lerp| (angular)")
plot_sphere(err_bilinear, fig=subfigs[1, 1], projection="orthographic", cmap="magma",
            colorbar=True, title="angular error: bilinear")
plot_sphere(err_slerp, fig=subfigs[1, 2], projection="orthographic", cmap="magma",
            colorbar=True, title="angular error: bilinear-spherical")
plt.show()

../_images/34753d668415fb2689fe5c18c33235b996ef96704a603108b3f2ea05ccb18715.png

for name, err in (("bilinear", err_bilinear), ("bilinear-spherical", err_slerp)):
    print(f"{name:20s} angular error:  max {math.degrees(err.max()):7.2f} deg"
          f"   mean {math.degrees(err.mean()):6.3f} deg"
          f"   points worse than 10 deg: {(err > math.radians(10)).sum().item()}")

bilinear             angular error:  max  135.22 deg   mean  2.867 deg   points worse than 10 deg: 4096
bilinear-spherical   angular error:  max    8.41 deg   mean  0.131 deg   points worse than 10 deg: 0

With \(\varphi = \mathrm{wrap}_\pi(4\,\mathrm{lat})\) the field jumps by \(2\pi\) wherever \(4\,\mathrm{lat} = \pm\pi\), that is at \(\pm 45^\circ\) — the two thin east–west stripes in the bilinear panel. There, bilinear averages \(+\pi\) and \(-\pi\), which denote the same direction, to \(0\), the opposite one, swinging the arrow half a turn the wrong way. Only the single cell straddling the jump is affected, which is why the damage shows up as a thin band rather than a smooth error. On those rings the spherical mode is exact.

The reference panel shows no stripe even though the field jumps there, because twilight is cyclic: \(+\pi\) and \(-\pi\) share a colour, exactly as they should for a direction.

The residual \(8.4^\circ\) reported for the spherical mode is neither interpolation error nor slerp’s fault — it sits on the polar rows and is identical for both modes. The equiangular output grid has rows at \(\pm 90^\circ\) that fall outside the Legendre–Gauss input’s latitude range, so those values come from the pole expansion reducing the adjacent ring, not from interpolation. Away from the poles the spherical mode reproduces this field exactly.

Note also that the error has to be measured modulo \(2\pi\): comparing raw values would flag correct answers as wrong wherever the two sit on opposite sides of the cut.

Conclusion#

Bilinear vs. bilinear-spherical: on scalar fields the two kernels agree exactly, as long as neighbouring samples differ by less than \(\pi\). The spherical mode exists for angle-valued fields, where it follows the shorter arc and so crosses the \(\pm\pi\) branch cut correctly instead of running the long way around the circle.
Interpolation vs. SHT upsampling: SHT-based upsampling is spectrally exact but prone to Gibbs ringing near discontinuities. Interpolation-based resampling trades spectral exactness for smoothness.