Interpolation Methods for Grid Shift Surfaces

Every grid-based datum shift ends in the same operation: four (or sixteen) numbers around a query position are combined into one shift value, and the choice of how they are combined sets a floor on the accuracy of everything downstream. This topic, part of core transformation fundamentals and standards, covers that operation in detail — the interpolation kernels, the boundary and no-data cases that break them, and how to measure the error the choice actually costs. The parent standards material assumes a grid has been selected, as described in NADCON vs NTv2: choosing the right datum shift; here we assume the grid is in hand and concentrate on reading it correctly.

Interpolation is worth this much attention because it is the one step in a grid shift where a defensible implementation and a plausible-looking one differ silently. A nearest-node lookup, a bilinear kernel and a bicubic kernel all return a number of the right order of magnitude at any query point. They differ by millimetres to centimetres, they differ most exactly where the shift surface is curving fastest, and nothing in the output distinguishes them.

The Shift Surface and Its Curvature #

A grid shift file stores a sampled version of a continuous distortion surface: the difference between two datums as a function of position. Between nodes the true surface is unknown, and interpolation is a model of it. The error of that model depends on two things — the node spacing $h$ and the curvature of the surface — and for bilinear interpolation the leading term is

Figure — interpolation error against node spacing: bilinear squares, bicubic goes as the fourth power.

where $s$ is the shift component being interpolated. Two consequences fall straight out. The error scales with the square of the spacing, so halving the node spacing quarters it — which is precisely why NTv2 supports nested sub-grids that refine over urban extents. And the error is proportional to curvature, so a smooth regional surface interpolates almost exactly while a surface with local structure — a fault, a historic re-survey boundary, a city adjusted separately — does not.

Bicubic interpolation uses a four-by-four neighbourhood and matches first derivatives at the nodes, reducing the leading error term to order $h^4$. It costs sixteen node reads instead of four and, importantly, it can overshoot: near a step in the surface a bicubic kernel produces values outside the range of its inputs, which a bilinear kernel never does. That trade-off is worked through in bilinear vs bicubic interpolation for grid shifts.

Kernels Compared #

The row that matters operationally is the second: bilinear is what the reference implementations use for horizontal grid shifts, so a manual reader that chooses anything else will disagree with PROJ by more than floating-point noise, and the disagreement will be read as a bug. If you deviate, say so in the audit record — the reasoning in when to use manual grid interpolation over pyproj applies directly.

Production Implementation: a Guarded Bilinear Reader #

from __future__ import annotations

from dataclasses import dataclass

import numpy as np

NULL_SENTINEL = -999.0          # NTv2 / NADCON no-data marker


@dataclass(frozen=True)
class ShiftGrid:
    """A regular shift surface: south-west origin, node spacing, and the values."""

    lat0: float                  # degrees, south-west node
    lon0: float                  # degrees, south-west node
    dlat: float                  # degrees per node, north-positive
    dlon: float                  # degrees per node, east-positive
    values: np.ndarray           # (nlat, nlon) float64, arc-seconds

    def __post_init__(self) -> None:
        if self.values.dtype != np.float64:
            raise ValueError("shift arrays must be float64: float32 leaks ~0.1 mm")
        if self.dlat <= 0 or self.dlon <= 0:
            raise ValueError("node spacing must be positive; flip the array instead")

    def cell(self, lat: float, lon: float) -> tuple[int, int, float, float]:
        """Locate the cell containing (lat, lon) and the fractional position in it."""
        nlat, nlon = self.values.shape
        fi = (lat - self.lat0) / self.dlat
        fj = (lon - self.lon0) / self.dlon
        # The extent test is inclusive of the last node, exclusive beyond it: a
        # query 1e-9 degrees outside is outside, and is a rejection, not a clamp.
        if not (0.0 <= fi <= nlat - 1 and 0.0 <= fj <= nlon - 1):
            raise ValueError(f"({lat:.8f}, {lon:.8f}) is outside the grid extent")
        i = min(int(np.floor(fi)), nlat - 2)
        j = min(int(np.floor(fj)), nlon - 2)
        return i, j, fi - i, fj - j

    def bilinear(self, lat: float, lon: float) -> float:
        """Bilinear shift in arc-seconds. Raises on a no-data corner."""
        i, j, u, v = self.cell(lat, lon)
        corners = np.array([
            self.values[i, j], self.values[i, j + 1],
            self.values[i + 1, j], self.values[i + 1, j + 1],
        ], dtype=np.float64)
        if np.any(corners == NULL_SENTINEL) or not np.all(np.isfinite(corners)):
            raise ValueError(
                f"({lat:.8f}, {lon:.8f}) sits on a cell with an unmodelled corner; "
                f"interpolating across a null node fabricates a shift"
            )
        w = np.array([(1 - u) * (1 - v), (1 - u) * v, u * (1 - v), u * v])
        return float(w @ corners)

The two guards are the whole point. An out-of-extent query and a no-data corner both have the same correct answer — refuse — and both have a tempting wrong answer that produces a plausible number. Those cases are covered in their own guides: handling grid edges and out-of-extent queries in Python and detecting null shift sentinels in grid files.

Figure — the path of one query, and the two places a correct reader refuses.

Precision and Tolerance #

The two columns differ by an order of magnitude at every spacing, which is the practical lesson: the published node spacing bounds the error only in combination with the smoothness of the surface, and the only way to know which column you are in is to measure against control, as measuring interpolation error against control points sets out.

Compliance and the Audit Record #

The interpolation is part of the operation, so it belongs in the record: the kernel used, the grid file and its checksum, the node spacing, and — where the deliverable is defended point by point — the indices of the nodes actually read. That last item is unusual and is the main practical reason to write a manual reader at all: it lets a reviewer reproduce a single coordinate from the grid by hand.

Validation #

Two checks catch nearly every interpolation bug. Interpolating at a node must return that node’s value exactly, to the last bit — it tests the indexing, the fractional position and the weight ordering in one line. Interpolating at a cell centre must return the mean of the four corners for a bilinear kernel, which tests the weights independently of the indexing.

def check_kernel(grid: ShiftGrid) -> None:
    at_node = grid.bilinear(grid.lat0 + grid.dlat, grid.lon0 + grid.dlon)
    assert at_node == grid.values[1, 1], "node lookup is not exact — check indexing"

    centre = grid.bilinear(grid.lat0 + 0.5 * grid.dlat, grid.lon0 + 0.5 * grid.dlon)
    expected = float(grid.values[0:2, 0:2].mean())
    assert abs(centre - expected) < 1e-12, "cell-centre value is not the corner mean"

Worked Example: One Cell, by Hand and by Code #

Take a 30-arc-second NTv2 sub-grid whose south-west node sits at 45.000000° N, 123.000000° W, and a query at 45.004000° N, 122.994000° W. The cell spans one node in each direction, so the fractional position inside it is u = 0.48 in latitude and v = 0.72 in longitude. With corner latitude-shift values of 0.2031, 0.2044, 0.2078 and 0.2091 arc-seconds (south-west, south-east, north-west, north-east), the bilinear weights are 0.1456, 0.3744, 0.1344 and 0.3456, and the interpolated shift is 0.20641 arc-seconds.

import numpy as np

grid = ShiftGrid(
    lat0=45.0, lon0=-123.0, dlat=30 / 3600, dlon=30 / 3600,
    values=np.array([[0.2031, 0.2044], [0.2078, 0.2091]], dtype=np.float64),
)
shift = grid.bilinear(45.004, -122.994)
print(f"{shift:.5f} arc-seconds = {shift * 30.87:.4f} m")
# 0.20641 arc-seconds = 6.3719 m

Two things are worth noticing in that arithmetic. The four corner values span 0.0060 arc-seconds — under two tenths of a millimetre of variation across the whole cell — so on this smooth stretch of the surface any reasonable kernel gives the same answer to well within tolerance. And the conversion from arc-seconds to metres is latitude-dependent: 30.87 m per arc-second is the meridian value near 45°, and using a single constant for both components is one of the ways an otherwise correct reader loses several millimetres. The reduction to metres is set out in projection math fundamentals for cadastral surveys.

Interpolating a Batch Without a Python Loop #

The per-point form above is the right shape for explaining the kernel and the wrong shape for a million parcel corners. The vectorised version computes the cell indices and fractional positions for every point as arrays, gathers the four corner values with fancy indexing, and forms the weighted sum in one expression — the same arithmetic, executed once over the whole batch:

def bilinear_batch(grid: ShiftGrid, lat: np.ndarray, lon: np.ndarray) -> np.ndarray:
    """Bilinear shift for arrays of positions. Same kernel, one pass."""
    nlat, nlon = grid.values.shape
    fi = (lat - grid.lat0) / grid.dlat
    fj = (lon - grid.lon0) / grid.dlon
    outside = (fi < 0) | (fi > nlat - 1) | (fj < 0) | (fj > nlon - 1)
    if np.any(outside):
        raise ValueError(f"{int(outside.sum())} point(s) outside the grid extent")

    i = np.clip(np.floor(fi).astype(np.int64), 0, nlat - 2)
    j = np.clip(np.floor(fj).astype(np.int64), 0, nlon - 2)
    u, v = fi - i, fj - j

    c00 = grid.values[i, j]
    c01 = grid.values[i, j + 1]
    c10 = grid.values[i + 1, j]
    c11 = grid.values[i + 1, j + 1]
    stack = np.stack([c00, c01, c10, c11])
    if np.any(stack == NULL_SENTINEL) or not np.all(np.isfinite(stack)):
        raise ValueError("batch touches a cell with an unmodelled corner")

    return ((1 - u) * (1 - v) * c00 + (1 - u) * v * c01
            + u * (1 - v) * c10 + u * v * c11)

The guards are deliberately kept as whole-array tests that raise for the batch rather than per-point exceptions: a batch that touches a null node or leaves the extent is a batch with a routing problem, and quietly dropping the offending rows is how a deliverable ends up with fewer parcels than it started with. The general form of that discipline — every row accounted for — is the subject of vectorizing coordinate transforms with NumPy and Dask.

Nested Sub-Grids: Refinement Where It Is Needed #

The error law says the only way to reduce interpolation error on a fixed surface is to reduce the node spacing, and the only way to reduce it everywhere is to store far more nodes than most of the surface needs. NTv2 resolves that with nesting: a coarse parent covers the whole extent, and finer children refine the parts where the surface has structure — typically urban areas, coastlines and zones with a complicated adjustment history.

Figure — nesting: the deepest covering sub-grid serves the query, not the first one found.

For a reader, nesting adds one rule and one failure mode. The rule is that the deepest covering sub-grid serves the query, because a child exists precisely because its parent was not good enough there. The failure mode is a reader that iterates the sub-grids in file order and takes the first match, which usually picks the parent — silently discarding the refinement and reintroducing the error the child was published to remove. The selection logic is set out in handling grid edges and out-of-extent queries in Python.

Nesting also changes what an audit record has to say. Two adjacent parcels can legitimately be served by different sub-grids at different spacings, which means their interpolation uncertainties differ, and a reviewer comparing them needs to know which grid produced each. Recording the sub-grid name alongside each transformed coordinate costs a string per row and answers the first question anyone asks when neighbouring parcels disagree.

Choosing a Kernel: a Short Decision Procedure #

Three questions settle the choice in nearly every case, in this order.

Is this a horizontal datum shift served by a published grid? If so, use bilinear, because the reference implementations do and agreement with them is worth more than a fraction of a millimetre of smoothness. Deviating is permitted and must be recorded.

Is this a smooth surface at coarse spacing — a geoid model, a velocity grid? Then a cubic kernel is defensible and often what the model’s own documentation specifies, since the leading error term falls from order $h^2$ to order $h^4$ and the surface has no steps for the cubic to ring against.

Does the surface contain steps or discontinuities? Then bilinear, regardless of spacing. A cubic kernel forced through a step overshoots outside the range of its inputs, and the overshoot lands next to the discontinuity — which is exactly where the parcels in dispute are. The comparison is worked through in bilinear vs bicubic interpolation for grid shifts.

Failure Modes #

• Row order assumed. NTv2 stores nodes south-to-north and east-to-west in longitude, with longitude positive west; a reader that assumes the array is north-to-south or west-to-east produces a mirrored surface. The symptom is a shift of the right magnitude in the wrong direction.
• Extent clamped instead of rejected. Clamping a query to the grid edge returns the boundary shift for a point that may be a hundred kilometres away. Reject, and let the fallback logic decide.
• No-data averaged. One -999.0 corner averaged with three valid ones moves the point about 7.5 kilometres — large, but not so large that it fails a range check on the coordinate.
• float32 storage. Grid values stored as float32 carry about 0.1 mm of representation error at arc-second magnitudes, which is small but entirely avoidable by promoting on load.
• Interpolating the wrong quantity. Latitude and longitude shifts must be interpolated independently; interpolating a combined magnitude and re-deriving the components loses the direction information.

Frequently Asked Questions #

Should I use bicubic interpolation to get a better answer than PROJ? #

Not for a horizontal datum shift. The reference implementations use bilinear, so a bicubic reader disagrees with every other tool by an amount larger than floating-point noise, and being different from the reference is a worse position to defend than being marginally less smooth. Bicubic is the right choice for geoid models and other smooth surfaces where the published guidance calls for it.

Does interpolation error accumulate over a large dataset? #

Not as a sum — it is a per-point error, not a drift, so a million points do not accumulate a million times the error. What does happen is that the error is spatially correlated: neighbouring points get similar errors, so a parcel boundary can be displaced coherently even though each point is individually within tolerance. That is why interpolation error belongs in the uncertainty budget rather than being dismissed as noise.

Can I resample a coarse grid to a finer one to reduce the error? #

No. Resampling invents nodes from the same information that was already there; the interpolation error is set by the original sampling and the true curvature of the surface, and no amount of upsampling recovers detail that was never captured. Resampling for performance is fine; resampling as an accuracy improvement is not.

How do I know if my grid is smooth or structured? #

Compute the second differences of the shift arrays across the grid and look at their distribution. A smooth regional surface has small, slowly varying second differences; a structured one has localised spikes at faults, coastlines and old adjustment boundaries. The spikes tell you where interpolation error is concentrated, and they are usually exactly where the interesting parcels are.

Related #

• Core transformation fundamentals and standards — the parent section
• Bilinear vs bicubic interpolation for grid shifts
• Handling grid edges and out-of-extent queries in Python
• Detecting null shift sentinels in grid files
• Understanding NTv2 grid shift files in Python