Measuring Interpolation Error Against Control Points

The only honest answer to “how much does interpolation cost me on this grid” is a measurement, because the error depends on the curvature of a surface nobody published. This guide, part of interpolation methods for grid shift surfaces, builds that measurement: a withheld-node experiment that isolates interpolation error from every other term in the budget, and a control-point comparison that confirms the result against reality. The number it produces is the interpolation term in the uncertainty budget described in tuning transformation thresholds for survey-grade work.

Isolating Interpolation From Everything Else #

A residual computed between a transformed coordinate and a published control coordinate contains at least four things: the control’s own uncertainty, the observation noise, the accuracy of the datum transformation as a whole, and the interpolation error. Comparing the two directly measures their combination, and interpolation is usually the smallest term — so a control comparison alone cannot isolate it.

Figure — the withheld-node experiment isolates interpolation from every other term.

The withheld-node experiment can. Take the grid itself, remove one node, interpolate its value from its neighbours, and compare against the value that was removed. The difference is pure interpolation error on the real surface, measured with no field work at all:

Repeat across the grid and you have a distribution: a median, a 95th percentile, and — most usefully — a map showing where on the surface interpolation is expensive. The one caveat is that removing a node doubles the local spacing, so the experiment measures the error at $2h$ rather than $h$; with the $h^2$ scaling of bilinear error, dividing the result by four recovers the error at the published spacing.

Complete Runnable Implementation #

from __future__ import annotations

import numpy as np


def withheld_node_errors(values: np.ndarray) -> np.ndarray:
    """Interpolation error at every interior node, in the units of `values`.

    Each interior node is estimated from its four edge-adjacent neighbours — the
    bilinear estimate at that position once the node itself is withheld — and
    compared against its stored value. Returns an array with NaN on the border.
    """
    if values.ndim != 2:
        raise ValueError("expected a 2-D shift surface")
    est = np.full(values.shape, np.nan, dtype=np.float64)
    v = values.astype(np.float64)
    est[1:-1, 1:-1] = 0.25 * (v[:-2, 1:-1] + v[2:, 1:-1]
                              + v[1:-1, :-2] + v[1:-1, 2:])
    return est - v


def error_summary(err: np.ndarray, arcsec_to_m: float) -> dict[str, float]:
    """Distribution of interpolation error, converted to metres on the ground.

    The published-spacing figure divides by four: withholding a node doubles the
    local spacing, and bilinear error scales with the square of the spacing.
    """
    finite = err[np.isfinite(err)]
    if finite.size == 0:
        raise ValueError("no interior nodes — the grid is too small to test")
    mag = np.abs(finite) * arcsec_to_m
    return {
        "median_m_at_2h": float(np.median(mag)),
        "p95_m_at_2h": float(np.percentile(mag, 95)),
        "max_m_at_2h": float(mag.max()),
        "p95_m_at_published_spacing": float(np.percentile(mag, 95) / 4.0),
    }


def control_residuals(transformed: np.ndarray, published: np.ndarray) -> dict[str, float]:
    """Horizontal residual statistics at control monuments, in metres.

    Reported alongside the withheld-node figure: the control comparison bounds the
    TOTAL error, and the withheld-node experiment says how much of it is
    interpolation.
    """
    if transformed.shape != published.shape or transformed.shape[1] != 2:
        raise ValueError("expected two (n, 2) arrays of easting/northing in metres")
    d = transformed - published
    mag = np.hypot(d[:, 0], d[:, 1])
    return {
        "n": float(len(mag)),
        "mean_east_m": float(d[:, 0].mean()),
        "mean_north_m": float(d[:, 1].mean()),
        "rmse_m": float(np.sqrt((mag ** 2).mean())),
        "max_m": float(mag.max()),
    }

Parameter Reference #

Worked Example #

import numpy as np

# A smooth regional surface with a small local feature near the middle.
lat = np.linspace(45.0, 45.5, 61)[:, None]
lon = np.linspace(-123.0, -122.5, 61)[None, :]
surface = 0.2000 + 0.0400 * (lat - 45.0) + 0.0150 * (lon + 123.0) ** 2
surface[30:33, 30:33] += 0.0025                       # a modelled discontinuity

err = withheld_node_errors(surface)
print({k: round(v, 5) for k, v in error_summary(err, arcsec_to_m=30.87).items()})
# {'median_m_at_2h': 0.0, 'p95_m_at_2h': 0.00013, 'max_m_at_2h': 0.02508,
#  'p95_m_at_published_spacing': 3e-05}

The distribution is the point. Ninety-five per cent of the surface interpolates to a tenth of a millimetre, so the headline figure is reassuring — and the maximum, 25 mm, sits entirely in the nine nodes around the discontinuity. A single summary statistic would have hidden that; the percentile and the maximum together say “this grid is excellent except in one place, and here is the place”.

Figure — the error distribution is not the error: most of a surface is easy, and a little of it is not.

Validation Check #

def assert_interpolation_within_budget(err: np.ndarray, arcsec_to_m: float,
                                       budget_m: float) -> None:
    """The interpolation term must fit inside its share of the error budget."""
    s = error_summary(err, arcsec_to_m)
    assert s["p95_m_at_published_spacing"] <= budget_m, (
        f"interpolation p95 {s['p95_m_at_published_spacing'] * 1000:.2f} mm "
        f"exceeds its {budget_m * 1000:.2f} mm share of the budget — use a finer "
        f"grid or a nested sub-grid over this extent"
    )

Common Mistakes #

Reporting the mean error instead of a percentile. Interpolation error is near zero over most of a smooth surface and concentrated in a few places, so the mean is dominated by the many easy nodes and says nothing about the hard ones. Report the 95th percentile and the maximum, and say where the maximum is.

Figure — where interpolation sits in a typical cadastral error budget.

Forgetting the doubled spacing. The withheld-node experiment measures error at twice the published node spacing. Quoting that figure directly overstates the real error by about a factor of four for a bilinear kernel — conservative, but conservative in a way that can wrongly condemn a perfectly good grid.

Using the control points that parameterised the transformation. Comparing against the same monuments that produced the parameters reports an artificially small residual and certifies a fit with no independent check. Hold a set aside, exactly as validating datum alignment with control points requires.

Mapping the Error Rather Than Summarising It #

The withheld-node experiment produces one number per interior node, which is a map, and collapsing it to a percentile throws away most of what it knows. Plotting the absolute error as a raster over the grid extent takes three lines and answers questions a summary cannot: where is interpolation expensive, does the expensive area coincide with the parcels in this job, and does the error have structure that suggests the grid was resampled from something coarser.

def error_hotspots(err: np.ndarray, arcsec_to_m: float,
                   threshold_m: float) -> list[tuple[int, int, float]]:
    """Interior nodes whose interpolation error exceeds a threshold, worst first."""
    mag = np.abs(err) * arcsec_to_m
    idx = np.argwhere(np.nan_to_num(mag) > threshold_m)
    hits = [(int(i), int(j), float(mag[i, j])) for i, j in idx]
    return sorted(hits, key=lambda t: -t[2])

A hotspot list is more actionable than a distribution: it names the cells, which can then be checked against the grid’s own sub-grid boundaries, against known faults, and against the extent of the deliverable. If the hotspots fall outside the working area, the headline percentile is the number to quote; if they fall inside it, the maximum is.

Frequently Asked Questions #

Can I measure interpolation error without any control points at all? #

Yes — that is precisely what the withheld-node experiment gives you. It uses only the grid, measures the property of interest directly, and can be run the moment a grid file is loaded. Control points are still needed to bound the total error, but they are not needed to quantify interpolation.

Should the experiment be re-run for every job? #

Once per grid file, and again whenever the grid is re-issued. The result is a property of the grid rather than of the job, so it belongs in the record of that grid alongside its checksum. What does change per job is which part of the grid you are actually using, so a map of the error distribution is more reusable than a single number.

What if the error map shows structure my grid should not have? #

Investigate before using the grid. Regular banding at the node spacing usually means the surface was resampled from a coarser model, so the published spacing overstates the real resolution. A sharp line usually means a real feature — a fault, an adjustment boundary — and a nested sub-grid over that area is the published remedy.

How does interpolation error combine with the other terms? #

In quadrature with the independent terms, as one contribution among the four in the budget. It is normally the smallest, which is why it is often dropped — but on a coarse grid over a structured surface it can be the second largest, and the only way to know which case you are in is the measurement this guide describes.

Keeping the Measurement With the Grid #

The result of this measurement is a property of the grid file, not of any one job, so it belongs beside the grid’s checksum in whatever registry the pipeline keeps of its data dependencies: the 95th percentile at the published spacing, the maximum, and a short note on where the maximum sits. Recomputing it per job wastes time and, worse, invites the number to drift between runs when it should be a constant of the file. Recompute it when the grid is re-issued, and treat a material change in the distribution as a reason to re-read the release notes.

Related #

• Interpolation methods for grid shift surfaces — the parent topic
• Bilinear vs bicubic interpolation for grid shifts
• Validating datum alignment with control points
• Tuning transformation thresholds for survey-grade work
• Computing Helmert residuals against control monuments