Fallback Routing Strategies for Missing Grid Files

Statutory boundary determination and cadastral surveying require deterministic coordinate operations with documented uncertainty budgets, the same audit-ready discipline established across Core Transformation Fundamentals & Standards. When a transformation pipeline encounters a missing grid shift file, silent degradation to a parametric approximation or an unhandled exception breaks ISO 19111 accuracy reporting and destroys auditability. A deterministic fallback routing strategy intercepts the missing asset, evaluates spatial coverage, enforces tolerance thresholds, and routes to a pre-validated transformation method before any coordinate enters a production dataset. This sub-task sits directly under the standards framework, where CRS definitions, operation methods, and grid availability are treated as explicit pipeline constraints rather than runtime afterthoughts.

Grid shift files (NTv2 .gsb, NADCON .las/.los, vertical .gtx) encode localized deformation surfaces that correct for crustal movement, historical network adjustments, and regional geodetic realignments. When the primary EPSG-registered grid is absent from the PROJ_DATA directory, the engine must traverse a compliance-driven routing table rather than guessing. The precedence is fixed: (1) the highest-resolution national or regional grid, (2) a coarser or legacy grid whose extent still covers the dataset, then (3) a parametric Helmert or Molodensky-Badekas transformation, accepted only while projected positional uncertainty stays within agency-defined tolerances. Cadastral workflows typically enforce ±0.025 m horizontal and ±0.050 m vertical limits; exceeding either triggers a hard pipeline failure that prevents contaminated geometries from propagating into statutory records. Which grid wins each tier depends on jurisdictional mandate and interpolation behaviour, the trade-off analysed in NADCON vs NTv2: choosing the right datum shift.

Figure — compliance-driven grid precedence: primary → secondary → parametric → hard fail.

How the Router Decides: Spec and Tolerance Math #

Before writing code, the routing rule must be expressed numerically so the decision is reproducible across machines. PROJ resolves a coordinate operation by searching PROJ_DATA (and any path in the PROJ_DATA/PROJ_LIB environment variables) for the grid named in the EPSG operation. If the grid is absent, PROJ may still return a ballpark operation — a rigid parametric shift with no local distortion modelling — which is exactly the silent degradation cadastral work cannot tolerate. The router therefore evaluates each candidate operation against a tolerance gate instead of trusting the library default. Building the candidate list with explicit accuracy values is covered in automating datum fallback chains in pyproj.

The total positional uncertainty of any fallback candidate is the root-sum-square of the operation’s own accuracy and the residual uncertainty of the control network it is tied to:

A candidate is admissible only when its total horizontal and vertical uncertainty both stay at or below the statutory tolerances:

where $\tau_h = 0.025\,\text{m}$ and $\tau_v = 0.050\,\text{m}$ for typical cadastral deliverables. Operation accuracy values come from the EPSG registry (the accuracy attribute of each coordinate operation), so the comparison is auditable rather than heuristic. Pinning $u_{net}$ to a real number requires independent check points, the workflow described in validating datum alignment with control points.

Step-by-Step Implementation #

The pipeline is assembled in three steps: validate the grid asset, encode the tolerance gate, then route through the precedence chain. Every snippet is self-contained and runnable on Python 3.10+ with pyproj and numpy installed.

Step 1 — Validate grid availability deterministically #

A missing file and a zero-byte or wrong-format file must both be rejected before PROJ is ever invoked. The check is purely local so it cannot itself introduce positional error.

from pathlib import Path

# ISO 19111-1:2019 §C.5 — a coordinate operation must reference a usable
# transformation resource; a missing or empty grid invalidates the operation.
VALID_GRID_SUFFIXES: frozenset[str] = frozenset(
    {".gsb", ".las", ".los", ".gtx", ".tif"}
)


def validate_grid_path(grid_path: Path) -> bool:
    """Return True only if the grid exists, is non-empty, and is a known format."""
    if not grid_path.exists() or grid_path.stat().st_size == 0:
        return False  # missing or truncated asset — never route to it
    return grid_path.suffix.lower() in VALID_GRID_SUFFIXES

Step 2 — Encode the tolerance gate #

The gate compares a candidate’s total uncertainty against the statutory limits using an explicit epsilon so a value sitting exactly on the boundary is treated deterministically rather than at the mercy of float rounding.

import math

EPSILON: float = 1e-9  # machine-precision slack for boundary equality


def within_tolerance(
    uncertainty_m: float, tolerance_h: float, tolerance_v: float
) -> bool:
    """ISO 19111 accuracy gate: accept a candidate only if it meets both limits."""
    # math.isclose handles the on-the-boundary case; the `<` covers clearly-inside.
    h_ok = math.isclose(uncertainty_m, tolerance_h, abs_tol=EPSILON) or uncertainty_m < tolerance_h
    v_ok = math.isclose(uncertainty_m, tolerance_v, abs_tol=EPSILON) or uncertainty_m < tolerance_v
    return h_ok and v_ok

Step 3 — Route through the precedence chain #

The router walks primary → secondary → parametric, applying validate_grid_path and within_tolerance at each tier and raising on tolerance exceedance so contaminated coordinates can never reach the dataset. Coordinates are rounded with explicit precision control to eliminate floating-point drift between runs.

import logging
from dataclasses import dataclass
from pathlib import Path
from typing import Optional, Sequence, Tuple
import numpy as np
from pyproj import CRS, Transformer
from pyproj.exceptions import ProjError

logging.basicConfig(format="%(asctime)s | %(levelname)s | %(message)s", level=logging.INFO)
logger = logging.getLogger("grid_fallback_router")

COORDINATE_PRECISION: int = 3  # millimetre resolution (1e-3 m) for survey-grade output


@dataclass(frozen=True)
class TransformationResult:
    coordinates: np.ndarray
    method: str
    uncertainty_m: float
    grid_file: Optional[str]
    iso_19111_compliant: bool


class GridFallbackRouter:
    def __init__(self, tolerance_horizontal: float = 0.025, tolerance_vertical: float = 0.050) -> None:
        self.tolerance_h = float(tolerance_horizontal)
        self.tolerance_v = float(tolerance_vertical)

    @staticmethod
    def _estimate_parametric_uncertainty(crs_from: CRS, crs_to: CRS) -> float:
        """Conservative u_total for a 7-parameter Helmert fallback (no local grid).

        Production code should read the EPSG operation `accuracy` attribute and
        combine it with the network uncertainty via root-sum-square.
        """
        return 0.150  # metres — deliberately pessimistic so the gate stays strict

    def _round_precision(self, coords: np.ndarray) -> np.ndarray:
        """Deterministic rounding to survey-grade precision (eliminates FP drift)."""
        return np.around(coords, decimals=COORDINATE_PRECISION)

    def _grid_transform(
        self,
        crs_src: CRS,
        crs_tgt: CRS,
        coords: np.ndarray,
        method: str,
        uncertainty_m: float,
        grid: Path,
    ) -> TransformationResult:
        # ISO 19111-1:2019 §10 — apply the declared grid-based operation.
        transformer = Transformer.from_crs(crs_src, crs_tgt, always_xy=True)
        x, y, z = transformer.transform(coords[:, 0], coords[:, 1], coords[:, 2])
        out = self._round_precision(np.column_stack((x, y, z)))
        logger.info("%s transformation succeeded via %s", method, grid.name)
        return TransformationResult(out, method, uncertainty_m, str(grid), True)

    def route_transformation(
        self,
        source_crs: str,
        target_crs: str,
        coordinates: Sequence[Tuple[float, float, Optional[float]]],
        primary_grid: Optional[Path] = None,
        secondary_grid: Optional[Path] = None,
    ) -> TransformationResult:
        coords = np.array(coordinates, dtype=np.float64)
        crs_src, crs_tgt = CRS.from_user_input(source_crs), CRS.from_user_input(target_crs)

        # Tier 1 — highest-resolution primary grid.
        if primary_grid is not None and validate_grid_path(primary_grid):
            try:
                return self._grid_transform(crs_src, crs_tgt, coords, "grid_shift_primary", 0.01, primary_grid)
            except ProjError:
                logger.warning("Primary grid load failed; evaluating secondary tier.")

        # Tier 2 — coarser/legacy grid with overlapping extent.
        if secondary_grid is not None and validate_grid_path(secondary_grid):
            try:
                return self._grid_transform(crs_src, crs_tgt, coords, "grid_shift_secondary", 0.03, secondary_grid)
            except ProjError:
                logger.warning("Secondary grid load failed; evaluating parametric tier.")

        # Tier 3 — parametric Helmert fallback, gated on tolerance.
        u_param = self._estimate_parametric_uncertainty(crs_src, crs_tgt)
        if within_tolerance(u_param, self.tolerance_h, self.tolerance_v):
            transformer = Transformer.from_crs(crs_src, crs_tgt, always_xy=True)
            x, y, z = transformer.transform(coords[:, 0], coords[:, 1], coords[:, 2])
            out = self._round_precision(np.column_stack((x, y, z)))
            logger.info("Parametric fallback accepted (u=%.4f m).", u_param)
            return TransformationResult(out, "parametric_helmert", u_param, None, True)

        # Tier 4 — hard fail: ISO 19111 forbids committing an out-of-tolerance result.
        raise RuntimeError(
            f"Transformation aborted: projected uncertainty {u_param:.4f} m exceeds "
            f"statutory tolerance (H {self.tolerance_h} m / V {self.tolerance_v} m). "
            "No valid grid shift file available; pipeline halted to preserve dataset integrity."
        )

Parameter and Return-Value Reference #

Every routing decision is driven by the inputs below. Units and valid ranges are fixed so the behaviour is reproducible in an audit.

Worked Example: NAD83(2011) → NAD83(HARN) in Oregon #

Consider a parcel corner at longitude −123.0294°, latitude 44.0521° (EPSG:6318, NAD83 2011) that must be expressed in NAD83(HARN) (EPSG:4152) for a county recorder who pins legacy plats to the HARN realization. The published NTv2 grid for this pair is us_noaa_nad83_2011_to_nad83_harn_or.gsb. With the grid present, the router applies it directly; with the grid missing, it must drop to the parametric tier and prove the result still clears tolerance.

from pathlib import Path

router = GridFallbackRouter(tolerance_horizontal=0.025, tolerance_vertical=0.050)

result = router.route_transformation(
    source_crs="EPSG:6318",   # NAD83(2011) geographic
    target_crs="EPSG:4152",   # NAD83(HARN) geographic
    coordinates=[(-123.0294, 44.0521, 95.250)],
    primary_grid=Path("/srv/proj_data/us_noaa_nad83_2011_to_nad83_harn_or.gsb"),
    secondary_grid=None,
)

print(result.method, result.uncertainty_m, result.grid_file)
# grid_shift_primary 0.01 /srv/proj_data/us_noaa_nad83_2011_to_nad83_harn_or.gsb

Against the published HARN monument for this corner, the grid-shifted longitude/latitude reproduce the control value to within ≈0.008 m — comfortably inside $\tau_h$. If the .gsb is deleted, validate_grid_path returns False, the primary and secondary tiers are skipped, and the router evaluates the 0.150 m parametric estimate against the 0.025 m gate. Because 0.150 m > 0.025 m, within_tolerance returns False and the router raises rather than emit a coordinate that would fail a later check against the monument. The mathematics behind that grid-shift correction is detailed in projection math fundamentals for cadastral surveys.

Verification and Residual Analysis #

A routed coordinate is only defensible once it has been checked against an independent monument and the residual logged as a structured audit record. The snippet below computes the horizontal residual magnitude, compares it to tolerance, and emits a machine-parseable log line.

import json
import logging
import numpy as np

logger = logging.getLogger("fallback_audit")


def verify_against_control(
    transformed_xy: tuple[float, float],
    control_xy: tuple[float, float],
    method: str,
    tolerance_h: float = 0.025,
) -> bool:
    """Compare a routed point to a control monument and log an audit record."""
    dx = transformed_xy[0] - control_xy[0]
    dy = transformed_xy[1] - control_xy[1]
    residual_m = float(np.hypot(dx, dy))  # RSS horizontal residual, metres
    passed = residual_m <= tolerance_h

    logger.info(
        json.dumps(
            {
                "method": method,
                "residual_m": round(residual_m, 4),
                "tolerance_h_m": tolerance_h,
                "within_tolerance": passed,
            }
        )
    )
    return passed


# Monument check for the worked example (projected metres, EPSG:6339 zone):
assert verify_against_control(
    transformed_xy=(497_882.114, 4_878_233.502),
    control_xy=(497_882.121, 4_878_233.498),
    method="grid_shift_primary",
)  # residual ≈ 0.008 m ≤ 0.025 m

For datasets spanning several monuments, aggregate the residual magnitudes into an RMSE and confirm the maximum residual stays below twice the RMSE before committing — the rigorous network treatment is covered in least-squares adjustment for control networks. Reference the ISO 19111:2019 geospatial information standard for the metadata fields each audit record must carry.

Troubleshooting and Gotchas #

Frequently Asked Questions #

Related References #

• Automating datum fallback chains in pyproj — building the accuracy-sorted candidate list this router walks.
• NADCON vs NTv2: choosing the right datum shift — how each tier’s grid is selected by jurisdiction and interpolation.
• Understanding NTv2 grid shift files in Python — parsing and validating the .gsb binaries the primary tier depends on.
• Validating datum alignment with control points — pinning the network uncertainty term in the tolerance gate.
• Core Transformation Fundamentals & Standards — the parent reference on CRS resolution, grid selection, and ISO 19111 compliance.