CoordTransform Cadastral & High-Precision Python

Handling Grid Edges and Out-of-Extent Queries in Python

The interesting half of a grid reader is what it does with the queries it cannot serve, and the correct answer — refuse — is the one that takes discipline to implement, because every alternative produces a number that looks like a coordinate. This guide, part of interpolation methods for grid shift surfaces, covers the four boundary cases a production reader meets: a query outside the grid entirely, a query in the outermost ring of cells, a query in the overlap between nested sub-grids, and a query exactly on a node or edge.

The Four Cases, and Why They Are Different

Outside the extent. The surface has no value there. Clamping to the nearest edge node returns a shift computed for a point that may be a hundred kilometres away, and does so silently. The correct response is to raise and let the routing layer decide, exactly as fallback routing strategies for missing grid files describes for a missing grid.

Four boundary cases and their correct handling Four cases. A query outside every extent is refused and routed to the fallback chain. A query in the outermost ring of cells is served bilinearly, with a recorded fallback if the intended kernel was bicubic. A query in a sub-grid overlap is served by the deepest covering grid. A query exactly on a node or a cell edge is served by clamping the index rather than the position, so it lands inside the last cell instead of being rejected. Correct handling The tempting error Outside every extent refuse and route clamp to the edge In the boundary ring bilinear, recorded silent kernel change In a sub-grid overlap deepest grid wins first match in file order Exactly on a node clamp the index index one past the end

Figure — four boundary cases, each with a different correct answer.

In the boundary ring. Bilinear interpolation works right up to the edge, because it only needs the cell it is in. A bicubic kernel does not: it needs a node of margin, so the outermost ring of cells is unservable. Falling back to bilinear there is the standard answer and must be recorded.

In a sub-grid overlap. NTv2 grids nest: a finer sub-grid refines part of its parent. A point inside a child must be served by the child, not the parent, and a reader that stops at the first matching sub-grid in file order can pick the wrong one. The rule is deepest-match-wins.

Exactly on a node or edge. A query on the boundary between two cells must give the same answer whichever cell serves it. Bilinear does this naturally when the weights are formed consistently; it is the floating-point comparison that decides which cell is chosen that breaks, and it breaks asymmetrically — a query at exactly the last node can index one past the end.

Complete Runnable Implementation 

from __future__ import annotations

from dataclasses import dataclass, field

import numpy as np


class OutOfExtent(ValueError):
    """Raised when no sub-grid covers a query position."""


@dataclass(frozen=True)
class SubGrid:
    """One NTv2-style sub-grid, with the parent it refines (None for a root)."""

    name: str
    parent: str | None
    lat0: float
    lon0: float
    dlat: float
    dlon: float
    values: np.ndarray
    depth: int = 0                      # 0 for a root grid, +1 per nesting level

    @property
    def extent(self) -> tuple[float, float, float, float]:
        nlat, nlon = self.values.shape
        return (self.lat0, self.lat0 + (nlat - 1) * self.dlat,
                self.lon0, self.lon0 + (nlon - 1) * self.dlon)

    def covers(self, lat: float, lon: float, eps: float = 1e-12) -> bool:
        """Inclusive of both edges, with a tolerance so a node is never 'outside'."""
        s, n, w, e = self.extent
        return (s - eps) <= lat <= (n + eps) and (w - eps) <= lon <= (e + eps)


@dataclass
class GridSet:
    """A collection of nested sub-grids, resolved deepest-first."""

    grids: list[SubGrid] = field(default_factory=list)

    def select(self, lat: float, lon: float) -> SubGrid:
        """The most refined sub-grid covering the position.

        Deepest match wins: a point inside a child must be served by the child,
        because that is the whole reason the child exists.
        """
        candidates = [g for g in self.grids if g.covers(lat, lon)]
        if not candidates:
            raise OutOfExtent(
                f"({lat:.8f}, {lon:.8f}) is not covered by any sub-grid in this file"
            )
        return max(candidates, key=lambda g: (g.depth, -g.dlat))

    def sample(self, lat: float, lon: float) -> tuple[float, str, str]:
        """Bilinear shift, plus the sub-grid and kernel actually used.

        Returning the provenance alongside the value is what lets a deliverable
        say WHICH grid produced a given coordinate — the question a reviewer asks
        first when two neighbouring parcels disagree.
        """
        grid = self.select(lat, lon)
        nlat, nlon = grid.values.shape
        fi = (lat - grid.lat0) / grid.dlat
        fj = (lon - grid.lon0) / grid.dlon
        # Clamp the INDEX, never the position: a query exactly on the last node is
        # inside the extent and must be served by the last cell, not rejected.
        i = min(max(int(np.floor(fi)), 0), nlat - 2)
        j = min(max(int(np.floor(fj)), 0), nlon - 2)
        u, v = fi - i, fj - j
        block = grid.values[i:i + 2, j:j + 2]
        if np.any(block == -999.0) or not np.all(np.isfinite(block)):
            raise OutOfExtent(
                f"({lat:.8f}, {lon:.8f}) sits on an unmodelled cell of {grid.name}"
            )
        w = np.array([(1 - u) * (1 - v), (1 - u) * v, u * (1 - v), u * v])
        return float(w @ block.reshape(4)), grid.name, "bilinear"

Parameter Reference

Name Type Units Note
lat, lon float degrees Longitude sign convention must match the grid’s
eps float degrees 1e-12 ≈ 0.1 µm; guards float equality at the edge
depth int — Nesting level; deepest match wins
return (shift, grid, kernel) tuple arc-seconds, str, str Provenance travels with the value
OutOfExtent exception — A routing signal, not a bug

Worked Example 

import numpy as np

parent = SubGrid("CANADA", None, 45.0, -124.0, 0.5, 0.5,
                 np.full((5, 5), 0.2000), depth=0)
child = SubGrid("VANCOUVER", "CANADA", 45.5, -123.5, 0.05, 0.05,
                np.full((5, 5), 0.2075), depth=1)
gs = GridSet([parent, child])

print(gs.sample(45.60, -123.40))    # inside the child
print(gs.sample(46.40, -121.00))    # parent only
try:
    gs.sample(50.00, -100.00)
except OutOfExtent as exc:
    print(f"refused: {exc}")
# (0.2075, 'VANCOUVER', 'bilinear')
# (0.2, 'CANADA', 'bilinear')
# refused: (50.00000000, -100.00000000) is not covered by any sub-grid in this file

The middle result is the one to check in a real file: a point that falls inside the parent but outside every child must be served by the parent, and a reader that only ever consults the deepest grid in the file will reject it.

Error from clamping a query to the grid edge Bar chart of the error in metres from clamping an out-of-extent query to the nearest edge node, at four distances beyond the extent: 0.03 metres at one kilometre, 0.15 at five, 0.31 at ten and 1.60 at fifty. The error is the change in the true shift surface over that distance, so it grows with the distance and with the local gradient — and nothing in the returned coordinate indicates any of it. 0.01 0.1 1 10 m 0.03 1 km 0.15 5 km 0.31 10 km 1.6 50 km

Figure — extrapolating past the edge: error against distance beyond the extent.

Validation Check 

def check_edges(gs: GridSet, grid: SubGrid) -> None:
    """A query exactly on each corner node must be served, not rejected."""
    s, n, w, e = grid.extent
    for lat, lon in ((s, w), (s, e), (n, w), (n, e)):
        value, used, _ = gs.sample(lat, lon)
        assert np.isfinite(value), f"corner ({lat}, {lon}) was not served"

Corner queries are the cheapest regression test for an off-by-one in the index clamp, and they fail loudly on the exact input a naive int(floor(...)) gets wrong.

Common Mistakes

Clamping the position instead of the index. Clamping the index keeps a query on the last node inside the last cell, which is correct. Clamping the position moves an out-of-extent query onto the boundary and returns a shift for somewhere else entirely — with no error, no warning, and a plausible-looking result.

Where the extent test belongs in a batch run Four ordered stages. Load the grid set and record its extents once. Test the whole input against those extents in one vectorised pass and report the count of unservable rows immediately. Route those rows to the fallback chain as a deliberate decision. Only then transform the servable rows. Discovering unservable rows at the end of a long run costs the run; discovering them in the first second costs nothing. 1 Load grids + extents 2 Test whole batch, one pass 3 Route unservable rows 4 Transform the rest

Figure — test extents over the whole batch first, not per point halfway through the run.

First-match sub-grid selection. Iterating the sub-grids in file order and taking the first that covers the point picks the parent whenever the parent comes first, which is usually. The nested child exists because the parent is not accurate enough there, so this silently discards the refinement the grid author went to the trouble of publishing.

Extent compared with a bare float equality. A query at exactly the northern limit can fail a strict < test by one unit in the last place after the extent was computed from lat0 + (n-1) * dlat. The tolerance in covers is there for that, and it is one part in 10¹² of a degree — far below any real position and far above float noise.

Frequently Asked Questions

Should an out-of-extent query be an exception or a sentinel return?

An exception, in a per-point API. A sentinel has to be checked by every caller and eventually will not be, whereas an exception propagates until something handles it deliberately. In a batch API the calculus changes: raise for the batch, but report which rows were affected, so the caller can route them rather than losing them.

What if a point is just barely outside the grid?

It is outside. There is no principled distance at which extrapolation becomes acceptable, and “just barely” is where the temptation is strongest. Route it to the fallback chain and record that the fallback was used — that is a defensible decision, and a two-metre extrapolation labelled as a grid shift is not.

How should overlapping sub-grids of the same depth be handled?

They should not overlap, and a file where they do is malformed — worth reporting rather than resolving silently. If one must be chosen, prefer the finer node spacing and record the choice; the tie-break in select does exactly that, and it is a diagnostic, not a feature.

Does rejecting queries slow down a large batch?

Not measurably, because the extent test is a handful of comparisons per point and vectorises trivially. What does cost time is discovering at the end of a two-hour run that 0.3 per cent of the rows were unservable. Test extents up front over the whole batch, before any transformation runs — the same “cheapest rejection first” ordering used in handling CRS mismatches in cadastral datasets.