Detecting Null Shift Sentinels in Grid Files

A node carrying a no-data sentinel is not a shift of −999 arc-seconds; it is the model saying it has nothing to offer at that position, and a reader that treats it as a number moves the point several kilometres or, worse, a few hundred metres. This guide, part of interpolation methods for grid shift surfaces, covers how sentinels appear in the formats a cadastral pipeline meets, how to detect them robustly across float representations, and why the correct response is a refusal rather than a substitution.

What a Sentinel Means and How It Is Stored #

Grid formats predate any convention for missing data in binary arrays, so they encode it in-band with a magic value. NTv2 and NADCON both use -999.0; some vertical models use -88.8888, others use large positive values, and a few modern products use IEEE NaN, which is the only self-describing choice of the four. The magic value is stored in the same float field as a real shift, so nothing but a comparison distinguishes them.

Figure — what a single unnoticed sentinel corner does, by its interpolation weight.

Two properties make the comparison less trivial than it looks. The value is usually stored as a 32-bit float, so the exact bit pattern read back is -999.0 only if the writer wrote exactly that; a value written as -999.0000001 in a producer’s double precision arrives as something a strict equality test misses. And in an array promoted to float64 for interpolation, a float32 -999.0 converts exactly — but a scaled array, one where the reader has already converted arc-seconds to radians or metres, no longer contains anything resembling −999 at all. Detect sentinels before any scaling, and detect them with a tolerance.

Complete Runnable Implementation #

from __future__ import annotations

from dataclasses import dataclass

import numpy as np

# Sentinels seen in the grid formats a cadastral pipeline actually meets.
KNOWN_SENTINELS = (-999.0, -88.8888, -9999.0)
SENTINEL_TOL = 1e-3          # generous: a real shift is never within 1e-3 of these


@dataclass(frozen=True)
class SentinelPolicy:
    """How a reader recognises and reports unmodelled nodes."""

    values: tuple[float, ...] = KNOWN_SENTINELS
    tol: float = SENTINEL_TOL

    def mask(self, arr: np.ndarray) -> np.ndarray:
        """Boolean mask of unmodelled nodes — True where the node has no value.

        Applied to the RAW array, before any unit conversion: once arc-seconds have
        been scaled to radians the sentinel is no longer recognisable.
        """
        if arr.dtype not in (np.float32, np.float64):
            raise TypeError("apply the sentinel mask to the raw float array")
        out = ~np.isfinite(arr)
        for s in self.values:
            out |= np.isclose(arr, s, rtol=0.0, atol=self.tol)
        return out


class UnmodelledNode(ValueError):
    """Raised when an interpolation cell touches a node with no value."""


def sanitise(arr: np.ndarray, policy: SentinelPolicy = SentinelPolicy()) -> np.ndarray:
    """Return a float64 copy with unmodelled nodes replaced by NaN.

    NaN is self-describing and propagates through arithmetic instead of silently
    contributing a number, so a bug downstream produces NaN rather than a
    plausible coordinate. The mask is taken BEFORE the promotion.
    """
    bad = policy.mask(arr)
    out = arr.astype(np.float64, copy=True)
    out[bad] = np.nan
    return out


def bilinear_strict(values: np.ndarray, fi: float, fj: float) -> float:
    """Bilinear sample that refuses any cell touching an unmodelled node."""
    nlat, nlon = values.shape
    i = min(max(int(np.floor(fi)), 0), nlat - 2)
    j = min(max(int(np.floor(fj)), 0), nlon - 2)
    block = values[i:i + 2, j:j + 2]
    if not np.all(np.isfinite(block)):
        n_bad = int((~np.isfinite(block)).sum())
        raise UnmodelledNode(
            f"cell ({i}, {j}) has {n_bad} unmodelled corner(s); the position is "
            f"outside the modelled area of this grid"
        )
    u, v = fi - i, fj - j
    w = np.array([(1 - u) * (1 - v), (1 - u) * v, u * (1 - v), u * v])
    return float(w @ block.reshape(4))


def coverage_report(values: np.ndarray) -> dict[str, float]:
    """Summary a pipeline can log once per grid, rather than per query."""
    bad = ~np.isfinite(values)
    return {
        "nodes": float(values.size),
        "unmodelled": float(bad.sum()),
        "unmodelled_fraction": float(bad.mean()),
    }

Parameter Reference #

Worked Example #

import numpy as np

raw = np.array([
    [0.2031, 0.2044, -999.0],
    [0.2078, 0.2091, -999.0],
    [0.2115, 0.2128, 0.2141],
], dtype=np.float32)

clean = sanitise(raw)
print(coverage_report(clean))
print(f"served cell: {bilinear_strict(clean, 1.5, 0.5):.5f}")
try:
    bilinear_strict(clean, 0.5, 1.5)
except UnmodelledNode as exc:
    print(f"refused: {exc}")
# {'nodes': 9.0, 'unmodelled': 2.0, 'unmodelled_fraction': 0.2222...}
# served cell: 0.20955
# refused: cell (0, 1) has 2 unmodelled corner(s); the position is outside ...

The refused cell is the one worth studying. Averaging its four corners — two real shifts near 0.204 and two sentinels — gives about −499.4 arc-seconds, which is roughly 15 kilometres of latitude. A reader that skipped the check would produce that number without comment, and a coordinate 15 km out is at least obvious. The dangerous variant is a single sentinel corner weighted at 0.05, which moves the point about 750 metres — enough to be badly wrong and small enough to be mistaken for a datum problem.

Figure — mask before converting units; after the conversion the sentinel is unrecognisable.

Validation Check #

def assert_no_sentinel_leak(values: np.ndarray) -> None:
    """After sanitising, no recognisable sentinel may remain as a finite value."""
    finite = values[np.isfinite(values)]
    for s in KNOWN_SENTINELS:
        assert not np.any(np.isclose(finite, s, rtol=0.0, atol=SENTINEL_TOL)), (
            f"sentinel {s} survived sanitisation — check the raw dtype and the "
            f"order of masking and unit conversion"
        )

Common Mistakes #

Testing equality against -999.0 on a float32 array. The comparison usually works and occasionally does not, depending on how the producer wrote the value and whether the array has been through any arithmetic. A tolerance of 1e-3 costs nothing — no genuine arc-second shift is within a thousandth of −999 — and removes the whole class of near-miss failures.

Figure — the sentinels a cadastral reader actually meets, and how each is detected.

Masking after unit conversion. Convert arc-seconds to radians first and the sentinel becomes −0.004847, which matches nothing in the sentinel list and passes every check. Mask on the raw array, then convert, and let NaN carry the information forward.

Substituting zero for a missing node. A zero shift is a statement: it says the two datums coincide at that node, which is almost never true and is precisely the claim the sentinel exists to avoid making. Substituting zero produces coordinates that look transformed and are not, and the error equals the full datum offset — metres. Refuse instead and let the fallback chain in fallback routing strategies for missing grid files decide.

Reporting Coverage Before the Batch Runs #

Sentinel handling per query is correct and late. A pipeline that transforms a million parcels should test coverage over the whole input against the sanitised grid before the first coordinate is transformed: intersect the input bounding box with the modelled area, count the rows whose cell touches an unmodelled node, and report the number up front. Discovering at the end of a long run that 0.4 per cent of the rows were unservable costs the whole run; discovering it in the first second costs nothing and turns an outage into a routing decision. The same argument, applied to grid extents rather than sentinels, is made in handling grid edges and out-of-extent queries in Python.

Frequently Asked Questions #

Why not interpolate using only the valid corners? #

Because the result would be an average of a smaller neighbourhood with different weights, which is a different estimator with a different and unquantified error — presented as though it were the same bilinear shift. If a reduced-neighbourhood estimate is genuinely wanted, it must be recorded as such, with its own uncertainty. In practice, a cell with a missing corner is at the edge of the modelled area, and the honest answer is that the model does not cover the point.

Should the accuracy nodes be checked for sentinels too? #

Yes. NTv2 records carry latitude and longitude accuracy values alongside the shifts, and those fields have their own sentinels. A pipeline that reports per-node accuracy without checking them will publish −999 as an accuracy figure, which is at least self-evidently wrong — but a mixed record where the shift is valid and the accuracy is not needs a decision, and the defensible one is to report the accuracy as unknown rather than to invent it.

How large a sentinel fraction should raise a flag? #

Any nonzero fraction is worth logging at load time; a fraction above a few per cent is worth a second look at whether the right grid was selected. A national grid whose coverage over your working area is mostly unmodelled is usually the wrong grid, and finding that out when the file is opened is much cheaper than finding it out per query, halfway through a batch.

Do modern formats still use magic numbers? #

Some do, and legacy files certainly do — and a cadastral pipeline reads legacy files for decades. Where a format uses NaN the detection is free and unambiguous, which is why the sanitiser converts everything to NaN internally: one representation to check for downstream, whatever the file used.

Related #

• Interpolation methods for grid shift surfaces — the parent topic
• Handling grid edges and out-of-extent queries in Python
• How to parse NTv2 GSB files with Python
• Python script for NADCON datum transformation
• Fallback routing strategies for missing grid files