healpix_resample.clough_tocher (Delaunay / Clough-Tocher C1 cubic regridding)#

healpix_resample.clough_tocher.CloughTocherResampler maps unstructured (lon_deg, lat_deg, val) samples onto a HEALPix grid using a Delaunay triangulation of the input samples, with a Clough-Tocher C1 cubic Bezier macro-element fitted per triangle — the same construction behind scipy.interpolate.griddata(method='cubic'). Unlike BicubicResampler, this is a genuine bivariate interpolant (exact at the input sample points), not a radial kernel sum.

Note on geodesy: the package manages the HEALPix authalic definition and the Earth ellipsoid with WGS84 through its geometry helper (healpix_geo).

Why this exists#

Comparing BicubicResampler against scipy.interpolate.griddata(method='cubic') on the same smooth, curved synthetic field shows griddata produces visibly fewer small-scale artifacts. The likely cause: BicubicResampler’s Npt=16 nearest-neighbour search is discrete — which cells belong to a given sample’s neighbourhood can flip between one output cell and its immediate neighbour, causing small jumps even though the weight formula itself is continuous for a fixed neighbour set. A Delaunay/Clough-Tocher construction doesn’t have this failure mode: adjacent triangles share two vertices, and the macro-element is built, by construction, to be C1 continuous (value and gradient) across shared edges.

What “Clough-Tocher” means here#

  1. Local gnomonic (central) projection. Delaunay triangulation and Clough-Tocher patches are inherently planar constructions, so samples are projected from (lon_deg, lat_deg) to a local tangent plane about the centroid of the sample set. Gnomonic projection is used specifically because it maps great-circle geodesics to straight lines, so triangle edges in the projected plane correspond to true geodesics on the sphere. This distorts for large angular extents — this resampler is intended for regional/local input extents (a single patch, not a global dataset in one construction call); __init__ raises if any sample falls outside the projection’s well-behaved hemisphere around the centroid. See docs/user-guide/regrid_to_healpix_parent_cell_subsetting.md for processing a large/global dataset one local patch at a time.

  2. Delaunay triangulation (scipy.spatial.Delaunay, CPU/NumPy — the one unavoidable non-torch step; no mature GPU-native Delaunay library exists) of the projected samples.

  3. Gradient estimation: a local least-squares plane fit at each triangulation vertex, using that vertex’s direct Delaunay 1-ring neighbours. Precomputed as two sparse (N, N) matrices Gx, Gy (geometry only). This is not claimed to be bit-identical to scipy.interpolate.CloughTocher2DInterpolator, which uses a different (Nielson-style) gradient estimator — it’s a correct, standard Clough-Tocher construction, not a scipy clone.

  4. Output validity = inside the convex hull. Like scipy.interpolate.griddata, this resampler does not extrapolate: a HEALPix cell is only kept in cell_ids if its projected center falls inside the Delaunay triangulation.

  5. Per-triangle Clough-Tocher macro-element: each triangle is split into 3 micro-triangles by connecting its vertices to the centroid; each micro-triangle carries a cubic Bezier patch built from corner values/gradients and a construction that guarantees C1 continuity across the internal split edges and across each triangle’s shared edges with its neighbours.

Everything except the Delaunay triangulation itself is precompiled into one sparse (N, K) matrix M, so resample(val) is a single batched sparse matmul (hval = y @ M), matching every other resampler in this package.

Class hierarchy and construction#

CloughTocherResampler does not subclass KNeighborsResampler — that base class’s KNN/Gaussian- threshold geometry model (Npt, sigma_m, threshold) doesn’t apply to a triangulation-based construction.

from healpix_resample import CloughTocherResampler

op = CloughTocherResampler(
    lon_deg=lon,       # (N,) sample longitudes, degrees -- regional/local extent
    lat_deg=lat,       # (N,) sample latitudes, degrees
    level=level,       # HEALPix level, nside = 2**level
    device="cuda",
    dtype=torch.float64,
    nest=True,
)

Key parameters#

  • lon_deg, lat_deg: unstructured sample coordinates in degrees, shape (N,). Must span a regional/local extent (see above).

  • level: HEALPix level (nside = 2**level).

  • nest, ellipsoid: same meaning as every other resampler.

  • device, dtype: PyTorch placement and numerical type for Gx, Gy, M.

  • out_cell_ids: optional caller-supplied subset of output cell ids, intersected with the convex-hull criterion.

  • candidate_ring_expand: how many fine-level HEALPix rings around the sample footprint to consider as candidate output cells before filtering to “inside the convex hull” — doesn’t affect correctness (over-generous candidates are simply dropped by the hull test), only whether the hull is fully covered.

Stored attributes#

  • N, K, cell_ids: as in every other resampler.

  • points2d: (N, 2) NumPy array — samples projected to the local gnomonic tangent plane.

  • tri: the underlying scipy.spatial.Delaunay triangulation.

  • Gx, Gy: sparse (N, N) torch tensors — grad_x = Gx @ f, grad_y = Gy @ f for any sample-space field f.

  • M: sparse CSR (N, K) torch tensor — hval = y @ M.

Methods#

resample(val)#

Sample space → HEALPix cell space: hval = val @ self.M. Accepts val of shape (N,) or (B, N), np.ndarray or torch.Tensor; returns the same type/batch shape, wrapped in ResampleResults(cell_data, cell_ids).

invert(hval)#

Not implemented. Unlike the KNN-based resamplers (which get a natural MT “for free” from a symmetric neighbour search), Delaunay/Clough-Tocher only defines a mapping from scattered samples to arbitrary query points, not the reverse. Calling invert() raises NotImplementedError with a pointer to planning/05_clough_tocher_resampler.md for the reasoning.

Notes and practical tips#

  • No benefit over BicubicResampler/BilinearResampler on a perfectly affine field — all interpolate a plane exactly; the advantage shows up on fields with genuine curvature.

  • Cells outside the sample convex hull are silently excluded from cell_ids, not given a garbage extrapolated value.

  • For a large/global dataset, use subset_for_parent_cell (see docs/user-guide/regrid_to_healpix_parent_cell_subsetting.md) to process one regional patch at a time — this resampler is not designed to be constructed once over a global dataset.