CoordTransform Cadastral & High-Precision Python

Automating Datum Fallback Chains in pyproj

Automating a datum fallback chain in pyproj is the operation of programmatically selecting the highest-accuracy coordinate operation that is actually available at runtime, then refusing to transform at all when no candidate clears a survey-grade tolerance — typically 0.025 m horizontal for cadastral deliverables under ISO 19111 coordinate-operation accuracy reporting. This page sits beneath fallback routing strategies for missing grid files, which defines the precedence policy; here we build the concrete TransformerGroup driver that walks that precedence deterministically and hard-fails rather than emitting a silently degraded ballpark coordinate.

pyproj datum fallback-chain selection Flowchart: TransformerGroup splits coordinate operations into an accuracy-sorted available .transformers list and an excluded .unavailable_operations list whose grids are missing. A loop walks the available candidates — while candidates remain it reads each candidate's declared accuracy and tests whether it is at or below the tolerance tau; if yes the operation is selected, applied and rounded to Decimal, if no the index advances and the loop repeats; when no candidate is left the driver raises RuntimeError as a deterministic hard fail rather than emit a degraded coordinate. TransformerGroup(src, tgt, always_xy) .transformers available · accuracy-sorted .unavailable_operations grid missing · excluded candidates left? candidate[i].accuracy accuracy ≤ τ ? Select · transform · commit Decimal round (half-even) raise RuntimeError hard fail · no degraded coord yes yes none left no · i++

What pyproj Actually Does When a Grid Is Missing

A common misconception is that pyproj returns the same operation with a “degraded” accuracy number when a grid shift file is absent. It does not. When pyproj.transformer.TransformerGroup is constructed for a datum pair, PROJ enumerates every registered coordinate operation between the two CRSs and sorts them by published accuracy. Any operation whose required grid (an NTv2 .gsb, NADCON .las/.los, or vertical .gtx) cannot be found in PROJ_DATA is moved out of the usable .transformers list and into .unavailable_operations. The .transformers list therefore already represents the fallback ordering — the best available operation is at index 0, and the .best_available flag tells you whether the globally best operation made the cut. Automating the chain means reading that structure explicitly instead of letting a bare Transformer.from_crs(...) quietly hand you whatever survived.

The admissibility test is numerical, not heuristic. Each candidate carries an EPSG-registered operation accuracy uopu_{\text{op}}; combined with the residual uncertainty unetu_{\text{net}} of the control network it ties to, the total standard uncertainty must stay at or below the statutory tolerance τ\tau:

utotal=uop2+unet2  ≤  τu_{\text{total}} = \sqrt{u_{\text{op}}^{2} + u_{\text{net}}^{2}} \;\le\; \tau

Because operation accuracy is sorted ascending, the first candidate that satisfies the gate is also the most accurate admissible one, so a single forward pass over .transformers is sufficient. A grid-free 7-parameter Helmert fallback typically reports 1–5 m, which fails a 0.025 m gate outright — exactly the degradation that must be rejected for legal boundary work rather than rounded and shipped.

Complete Runnable Implementation

The function below instantiates the operation graph, walks the accuracy-sorted available candidates, applies a hard tolerance gate, and converts every PROJ float output to Decimal with round-half-to-even before returning. It raises rather than degrade. It runs as-is on Python 3.10+ with pyproj installed.

import logging
from decimal import Decimal, ROUND_HALF_EVEN, getcontext
from typing import List, Optional, Sequence, Tuple

from pyproj import CRS, Transformer
from pyproj.transformer import TransformerGroup
from pyproj.exceptions import ProjError

getcontext().prec = 18  # high-precision decimal context for cadastral arithmetic
logger = logging.getLogger(__name__)


def execute_datum_fallback_chain(
    source_crs: str,
    target_crs: str,
    max_tolerance_m: Decimal,
    coordinates: Sequence[Tuple[float, float]],
    precision_places: int = 4,
) -> List[Tuple[Decimal, Decimal]]:
    """Select the best AVAILABLE coordinate operation that clears tolerance, then transform.

    Mirrors the ISO 19111-1:2019 accuracy-reporting model: every candidate operation
    carries a declared accuracy, and a result may only be committed when that accuracy
    is at or below the statutory tolerance. Operations whose grid shift file is missing
    are excluded by PROJ into TransformerGroup.unavailable_operations and never selected.

    Raises RuntimeError if no available candidate meets max_tolerance_m (deterministic
    hard fail — no silent degradation to a ballpark Helmert shift).
    """
    if max_tolerance_m <= 0:
        raise ValueError("max_tolerance_m must be a positive Decimal value.")

    # 1. Strict CRS instantiation (EPSG code or urn:ogc:def:crs:EPSG::XXXX).
    try:
        src = CRS.from_user_input(source_crs)
        tgt = CRS.from_user_input(target_crs)
    except ProjError as exc:
        raise ValueError(f"Invalid CRS definition: {exc}") from exc

    # 2. Build the operation graph. Only grid-present / grid-free ops appear in
    #    .transformers; grid-blocked ops fall into .unavailable_operations.
    group = TransformerGroup(src, tgt, always_xy=True)
    if not group.best_available:
        logger.warning(
            "Best operation is unavailable (%d blocked by missing grids); "
            "evaluating degraded fallbacks against tolerance.",
            len(group.unavailable_operations),
        )

    # 3. Walk accuracy-sorted available candidates and take the first that passes.
    selected: Optional[Transformer] = None
    selected_accuracy: Optional[Decimal] = None
    for transformer in group.transformers:
        reported = transformer.accuracy  # metres; None or -1.0 means "unknown"
        if reported is None or reported < 0:
            logger.info("Skipping operation with unknown accuracy metric.")
            continue
        acc = Decimal(str(reported))
        if acc <= max_tolerance_m:
            selected, selected_accuracy = transformer, acc
            break
        logger.info("Rejected op: accuracy %s m exceeds tolerance %s m.", acc, max_tolerance_m)

    # 4. Hard stop — ISO 19111 forbids committing an out-of-tolerance coordinate.
    if selected is None:
        raise RuntimeError(
            f"No available transformation meets {max_tolerance_m} m tolerance "
            f"({len(group.unavailable_operations)} operations blocked by missing grids). "
            "Pipeline halted to prevent silent positional degradation."
        )

    logger.info("Selected operation with accuracy %s m.", selected_accuracy)

    # 5. Transform with explicit Decimal rounding to eliminate float drift.
    quant = Decimal(10) ** -precision_places
    out: List[Tuple[Decimal, Decimal]] = []
    for x, y in coordinates:
        tx, ty = selected.transform(x, y)
        dx = Decimal(str(tx)).quantize(quant, rounding=ROUND_HALF_EVEN)
        dy = Decimal(str(ty)).quantize(quant, rounding=ROUND_HALF_EVEN)
        out.append((dx, dy))
    return out

The always_xy=True argument forces longitude/latitude (easting/northing) input order across every candidate, removing the authority-axis ambiguity that otherwise makes the chain non-portable — the same axis-locking discipline used when setting up high-precision coordinate reference systems.

Parameter and return reference

Name Type Units Valid range / meaning
source_crs / target_crs str EPSG code or URN any registered CRS; URN form locks the datum realization
max_tolerance_m Decimal metres > 0; statutory gate, e.g. 0.025 cadastral, 0.05 engineering control
coordinates Sequence[(x, y)] degrees or metres within the operation’s area of use; (lon, lat) order under always_xy
precision_places int decimal places 0–9; 4 ≈ 0.0001 m, match jurisdictional rounding rules
returns List[(Decimal, Decimal)] same as target CRS rounded round-half-to-even; reproducible across machines
raises RuntimeError — — no available candidate clears the tolerance gate

Minimal Worked Example

Transform a parcel corner from NAD27 to NAD83, demanding 0.05 m accuracy. If the NADCON grids are installed, the chain selects the grid-based operation; if they are absent, every available candidate is a coarse parametric shift and the call raises:

from decimal import Decimal

corner = [(-122.6765, 45.5231)]  # (lon, lat) — Portland, OR area parcel corner

shifted = execute_datum_fallback_chain(
    source_crs="EPSG:4267",      # NAD27 geographic
    target_crs="EPSG:4269",      # NAD83 geographic
    max_tolerance_m=Decimal("0.05"),
    coordinates=corner,
    precision_places=6,
)
print(shifted)
# -> [(Decimal('-122.676497'), Decimal('45.523113'))]
#    selected op accuracy ~= 0.02 m (NADCON5) <= 0.05 m gate

With the grids present the log records the selected operation at roughly 0.02 m, well inside the gate. Delete the NADCON grids from PROJ_DATA and the only survivors in .transformers are metre-level parametric shifts, so the gate rejects them all and the function raises RuntimeError instead of returning a coordinate that would fail a later check against the monument.

Validation Check

Gate the selection before any coordinate is written downstream — assert that an available operation exists and that its declared accuracy clears tolerance:

from decimal import Decimal
from pyproj import CRS
from pyproj.transformer import TransformerGroup

group = TransformerGroup(CRS("EPSG:4267"), CRS("EPSG:4269"), always_xy=True)
best = group.transformers[0]

assert group.transformers, "No available operation: all candidates blocked by missing grids."
assert best.accuracy is not None and Decimal(str(best.accuracy)) <= Decimal("0.05"), (
    f"Best available accuracy {best.accuracy} m exceeds survey-grade tolerance."
)

For multi-monument datasets, follow the selection with an independent residual check against control and aggregate to an RMSE before committing — the rigorous network treatment lives in least-squares adjustment for control networks, and the tolerance values themselves are derived in optimizing transformation tolerance thresholds.

Common Mistakes

Trusting Transformer.from_crs instead of inspecting TransformerGroup
A bare Transformer.from_crs(src, tgt) silently returns the single best available operation and discards the accuracy metadata. When the grid is missing it hands back a ballpark parametric shift with no warning. Always build a TransformerGroup, read best_available and each transformer's .accuracy, and gate explicitly so a metre-level fallback can never masquerade as a grid-shifted result.
Treating an unknown accuracy (None / -1.0) as zero
PROJ reports accuracy is None or -1.0 when an operation declares no uncertainty. Coercing that to 0.0 makes an undocumented operation look perfect and slip through the gate. Skip unknown-accuracy candidates outright — an operation with no declared accuracy cannot satisfy an ISO 19111 accuracy assertion.
Letting float output decide the final digit
PROJ returns IEEE 754 float64, which serialises differently across platforms and accumulates drift over chained operations. Convert each output to Decimal(str(value)) and quantize with ROUND_HALF_EVEN before output so two runs on different machines produce byte-identical, legally defensible coordinates.