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.

Anatomy of a TransformerGroup result A TransformerGroup exposes four members. The transformers list holds usable operations sorted by declared accuracy. The unavailable_operations list holds operations whose grid files are missing from the PROJ data directory. The best_available flag is true only when the first entry of transformers is also the most accurate operation known. Each transformer carries an accuracy in metres, which is the value the tolerance gate compares against. .transformers list runnable, accuracy-sorted .unavailable_operations list grid file not on disk .best_available bool False means a better op exists .transformers[i].accuracy float | None metres; None means unstated

Figure — what a TransformerGroup hands back, and which half of it can actually run.

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:

Declared accuracy of each candidate operation Bar chart of the declared accuracy in metres for five candidate operations returned for one NAD83 realisation pair: NTv2 grid 0.015, NADCON pair 0.05, published Helmert 0.3, geocentric translation 1.0, and a null operation at 2.0 metres. A tolerance of 0.05 metres admits only the first two. 0.01 0.1 1 10 m 0.015 NTv2 grid 0.05 NADCON 0.3 Helmert 1 translation 2 null op tolerance tau = 0.05 m

Figure — a typical candidate list: only the operations left of the gate may run.

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.

Frequently Asked Questions

Why not just let the library pick the best available operation?

Because ‘best available’ silently becomes ‘best of what is left’ when a grid file is missing, and nothing in the return value tells the caller that happened. The result is a coordinate that is a decimetre or a metre out, in a file that looks exactly like a correct one. Inspecting the operation group makes the choice visible and lets you refuse when the best remaining option is not good enough for the job.

What should happen when no candidate meets the tolerance?

Raise. It is tempting to return the best available result with a warning, but a warning in a log is not a property of the data, and by the time the coordinates reach a deliverable nobody remembers it. A hard failure forces the real fix, which is either installing the missing grid or accepting a documented lower accuracy as a deliberate, recorded decision.

How do I make grid availability reproducible across machines?

Pin it. Record which grid files were present and their checksums in the audit record, and provision them explicitly rather than relying on whatever a developer happened to download. A pipeline whose accuracy depends on the contents of a cache directory will produce different numbers on the build server than on a workstation, and the difference will not show up until someone compares two deliverables.