Spatial Pipeline EngineeringMRV Automation

Debugging Silent Datum Shifts in Carbon Pipelines

A datum shift is the most dangerous class of coordinate error in a carbon MRV stack because it is silent: the geometry reprojects without an exception, the map overlays look right at project zoom, and the reported area is off by a fraction of a percent that no eyeball catches. This guide is the failure-mode companion to the Geospatial Coordinate Reference Systems (CRS) Alignment discipline within the MRV Architecture & Carbon Accounting Fundamentals stack, and where its sibling how to align WGS84 to a local CRS in Python for carbon mapping shows the correct construction, this page shows how to prove one is wrong.

The failure is not that coordinates are missing a CRS — that raises loudly and gets fixed on day one. The failure is that a geometry carries a plausible but incorrect datum, or that PROJ quietly substitutes a ballpark transform when the proper NTv2 or NADCON grid is absent, and the parcel lands one to tens of metres from where it belongs. At that magnitude a boundary crosses a stratum edge, a pixel of high-biomass forest is credited to a low-biomass class, and the error propagates untouched into the lineage record that MRV data lineage and provenance tracking will later attest as correct. Detection has to be deliberate, because nothing in the happy path will surface it.

Decision tree for detecting a silent datum shift Geometry that looks aligned enters a four-check diagnostic. First, is the CRS authority and datum ensemble explicit? If not, it is flagged as a datum-shift risk. Next, is the exact NTv2 or NADCON grid available to pyproj? If not, only a ballpark transform exists and the geometry is flagged. Then measured displacement between two datum realizations is compared against a tolerance gate. Displacement under tolerance passes to a verified-datum output; displacement over tolerance is flagged as a datum shift. The verified-datum output is drawn in amber. Geometry looks aligned authority + datum ensemble explicit? NTv2 / NADCON grid available to pyproj? displacement under tolerance? Verified datum realization pinned Flag: silent datum shift quarantine · log realization · block credit yes yes yes no no no

Root Cause Analysis

Silent datum shifts arise from a small number of distinct mistakes, each of which produces a geometrically valid but positionally wrong result. Understanding which one you are looking at determines how you detect it.

First, treating EPSG:4326 as if it were a local datum. WGS84 is a global reference frame; NAD83, ETRS89, GDA2020 and the many national datums are not aligned to it. In North America the WGS84-to-NAD83 offset is on the order of one to two metres and grows over time; in Australia the historical AGD66/AGD84-to-GDA offset exceeds 100 metres. A pipeline that ingests GPS coordinates as EPSG:4326 and then labels them NAD83 without a transform has not moved anything — it has simply mislabelled a one-to-two-metre error as zero. Nothing raises, because both are valid geographic CRSs with the same axis order and units.

Second, an unresolved datum ensemble late-binding to the wrong realization. Modern PROJ models EPSG:4326 and EPSG:4269 as datum ensembles with several member realizations (for WGS84, the G-series frames aligned to successive ITRF epochs). When you ask for a transform without pinning a realization, PROJ late-binds: it defers the choice and may select a null or identity step, assuming the ensemble members are interchangeable. Across ITRF realizations that assumption is worth several centimetres to a few decimetres per decade of plate motion — small per parcel, but systematic across an entire project and directional, so it does not average out over a credited area.

Third, a missing NTv2 or NADCON grid silently downgrading to a ballpark transform. The accurate datum transform for many regions is a gridded correction (.gsb NTv2 files, NADCON grids). If that grid is not installed in the PROJ data directory, PROJ does not fail — it falls back to a lower-accuracy Helmert or null transform and reports a coarser accuracy in metres. The geometry still reprojects. The only evidence is the accuracy figure buried in the transform’s metadata, which a naive to_crs() call never inspects. This is the failure mode that most often reaches production, because it depends on the runtime environment: it passes on a developer’s fully provisioned laptop and fails silently in a stripped container.

Diagnostic Pipeline: Detecting the Shift Before Area Is Computed

The diagnostic inspects candidate transforms rather than trusting the default. It uses pyproj.transformer.TransformerGroup to enumerate every available operation, detects when the best available transform is only a ballpark (grid-missing) fallback, and measures the actual metre displacement between two datum realizations over the project’s area of interest. All three checks run before any area is computed.

from __future__ import annotations

import geopandas as gpd
import numpy as np
import structlog
from pyproj import CRS
from pyproj.aoi import AreaOfInterest
from pyproj.transformer import Transformer, TransformerGroup

log = structlog.get_logger()

# Displacement beyond this (metres) means the datum realizations disagree
# materially over the project extent — a shift we must not credit against.
DISPLACEMENT_TOLERANCE_M = 0.5


def diagnose_datum_shift(
    gdf: gpd.GeoDataFrame,
    source_crs: str,
    target_crs: str,
    aoi: tuple[float, float, float, float],  # west, south, east, north (deg)
) -> dict:
    """Detect grid-missing (ballpark) transforms and measure datum displacement.

    Returns a report; never mutates the geometry. Raises nothing so callers
    can gate on the flags rather than on an exception.
    """
    area = AreaOfInterest(*aoi)
    src, dst = CRS.from_user_input(source_crs), CRS.from_user_input(target_crs)

    # 1. Enumerate every candidate operation for this pair over the AOI.
    group = TransformerGroup(src, dst, area_of_interest=area, always_xy=True)
    if not group.transformers:
        return {"status": "no_transform", "flag": True}

    best = group.transformers[0]
    # PROJ marks a coarse fallback as "ballpark"; a missing NTv2/NADCON grid
    # is the usual cause. group.unavailable_operations lists what PROJ *wanted*
    # but could not use because the grid file is absent from the data dir.
    grid_missing = len(group.unavailable_operations) > 0
    is_ballpark = bool(best.is_network_enabled is False and best.description
                       and "ballpark" in best.description.lower())

    # 2. Build a grid-aware transformer restricted to non-ballpark ops, and a
    #    strict one that refuses any ballpark substitution, then compare.
    tf_default = Transformer.from_crs(src, dst, always_xy=True,
                                      area_of_interest=area)
    tf_strict = Transformer.from_crs(src, dst, always_xy=True,
                                     area_of_interest=area, only_best=True)

    # 3. Measure displacement at representative vertices of every parcel.
    pts = np.array([(x, y) for geom in gdf.geometry
                    for x, y in geom.exterior.coords[:-1]])
    lon, lat = pts[:, 0], pts[:, 1]

    dx_def, dy_def = tf_default.transform(lon, lat)
    dx_str, dy_str = tf_strict.transform(lon, lat)
    disp = np.hypot(np.asarray(dx_def) - np.asarray(dx_str),
                    np.asarray(dy_def) - np.asarray(dy_str))
    max_disp = float(np.nanmax(disp))

    report = {
        "status": "diagnosed",
        "source_crs": src.to_authority(),
        "target_crs": dst.to_authority(),
        "best_transform": best.description,
        "best_accuracy_m": best.accuracy,
        "grid_missing": grid_missing,
        "is_ballpark": is_ballpark,
        "max_displacement_m": round(max_disp, 4),
        "flag": grid_missing or is_ballpark
        or max_disp > DISPLACEMENT_TOLERANCE_M,
    }
    if report["flag"]:
        log.warning("datum.shift_suspected", **report)
    else:
        log.info("datum.clean", **report)
    return report

Two signals matter most. grid_missing (derived from TransformerGroup.unavailable_operations) tells you PROJ knew a better transform existed but could not load its grid — the container-provisioning failure. max_displacement_m quantifies how far the geometry actually moves between the permissive default and the only_best=True strict transform; when the default has silently fallen back to a ballpark, this gap is the datum shift, measured in metres.

Deterministic Transformation Logic

Once the diagnostic has confirmed the environment can supply an accurate transform, the transformation step enforces it. It resolves the source ensemble to a specific realization, constructs a transformer with only_best=True so PROJ raises rather than substituting a ballpark, re-measures displacement against the pre-transform vertices, and asserts the maximum is under tolerance before returning. There is no path by which a coarse transform reaches the credited-area calculation.

from __future__ import annotations

import geopandas as gpd
import numpy as np
import structlog
from pyproj import CRS
from pyproj.aoi import AreaOfInterest
from pyproj.transformer import Transformer

log = structlog.get_logger()

MAX_SHIFT_M = 0.5  # tolerance gate on residual datum displacement


def reproject_with_datum_gate(
    gdf: gpd.GeoDataFrame,
    source_realization: str,   # e.g. "EPSG:9057" (WGS84 G1762), NOT the ensemble
    target_crs: str,           # projected, area-honest target e.g. "EPSG:6933"
    aoi: tuple[float, float, float, float],
) -> gpd.GeoDataFrame:
    """Apply an explicit, grid-based transform with a displacement gate.

    Fails closed: a missing grid or an over-tolerance shift raises before any
    geometry is returned for area accounting.
    """
    area = AreaOfInterest(*aoi)
    src = CRS.from_user_input(source_realization)
    if src.is_deprecated:
        raise ValueError(f"deprecated CRS {source_realization}; pin a live realization.")

    # only_best=True: PROJ refuses to silently downgrade to a ballpark transform.
    # If the NTv2/NADCON grid is absent, this raises instead of guessing.
    try:
        transformer = Transformer.from_crs(
            src, CRS.from_user_input(target_crs),
            always_xy=True, area_of_interest=area, only_best=True,
        )
    except Exception as exc:  # pyproj raises when no non-ballpark op is available
        log.error("datum.grid_unavailable", source=source_realization, error=str(exc))
        raise RuntimeError(
            "required datum-transform grid is not installed; provision PROJ data "
            "before running credited-area accounting.") from exc

    # Sanity-measure the shift on a geographic round-trip before committing.
    probe = gdf.to_crs(4326)
    lon = np.array([geom.centroid.x for geom in probe.geometry])
    lat = np.array([geom.centroid.y for geom in probe.geometry])
    fx, fy = transformer.transform(lon, lat)
    # Compare against a strict no-op ballpark to expose residual disagreement.
    ballpark = Transformer.from_crs(src, CRS.from_user_input(target_crs),
                                    always_xy=True, area_of_interest=area)
    bx, by = ballpark.transform(lon, lat)
    residual = float(np.nanmax(np.hypot(np.asarray(fx) - np.asarray(bx),
                                        np.asarray(fy) - np.asarray(by))))
    if residual > MAX_SHIFT_M:
        log.warning("datum.gate_breached", residual_m=round(residual, 4),
                    tolerance_m=MAX_SHIFT_M)
        raise ValueError(
            f"datum displacement {residual:.3f} m exceeds gate {MAX_SHIFT_M} m; "
            "the ballpark and grid transforms disagree — do not credit this parcel.")

    out = gdf.to_crs(target_crs)
    out.attrs["datum_transform"] = transformer.description
    out.attrs["source_realization"] = src.to_authority()
    out.attrs["residual_shift_m"] = round(residual, 4)
    log.info("datum.reprojected", target=target_crs,
             transform=transformer.description, residual_m=round(residual, 4))
    return out

The load-bearing choices are only_best=True, which converts a silent grid-missing fallback into a raised exception, and the pinned source_realization (a specific EPSG code such as EPSG:9057 for WGS84 G1762) instead of the ensemble EPSG:4326. Together they make the transform reproducible and force any environment that cannot supply the accurate grid to fail closed rather than credit displaced area.

Compliance Gating & Audit Trail Generation

A datum-shift gate is only defensible if its verdict is recorded. Every parcel that passes carries the pinned realization, the exact transform description, and the measured residual shift into the lineage record; every parcel that fails is quarantined with the same evidence, so a verifier can distinguish a geometry that was checked and cleared from one that was never tested. Those records are what MRV data lineage and provenance tracking exposes to an auditor, and what ground-truth alignment for carbon models relies on when it reconciles field-plot coordinates against modelled geometry — a plot-to-pixel comparison is meaningless if the two sit in different realizations.

The gates a submission must show:

  1. Realization pinned, not ensemble. The exported metadata names a specific datum realization and transform operation code, so the reprojection is reproducible byte-for-byte. An ensemble label alone is treated as a finding, because it admits late binding.
  2. Grid provenance recorded. The NTv2/NADCON grid file and its version are logged. A run whose grid_missing flag was ever true is blocked from crediting, satisfying ISO 14064-3’s requirement that geolocation be traceable and reproducible.
  3. Displacement under tolerance. The residual shift in metres is stored per parcel. Verra VM-series methodologies size credited area from geometry; a documented sub-tolerance residual demonstrates the area was not inflated or deflated by a coordinate artefact.
  4. Quarantine transparency. Flagged parcels are held with their diagnostic report attached rather than dropped, so a CSRD ESRS E1 reviewer sees the exclusion and its reason instead of an unexplained gap in the reported extent.

Production Integration

Wire the diagnostic and the gated transform into the pipeline in a fixed ingest → diagnose → transform → validate → export → submit sequence, mirroring the HowTo steps above:

  1. Ingest. Read parcels with geopandas, and resolve each source CRS to a specific realization; refuse any layer whose CRS is an unpinned ensemble, following the construction rules in how to align WGS84 to a local CRS in Python for carbon mapping.
  2. Diagnose. Run diagnose_datum_shift over the project area of interest to enumerate candidate transforms, catch a grid-missing ballpark fallback, and measure inter-realization displacement.
  3. Transform. Call reproject_with_datum_gate with only_best=True so a missing grid raises rather than substituting a coarse transform.
  4. Validate. Assert every parcel’s residual shift is under the tolerance gate; route breaches to quarantine with their diagnostic report intact.
  5. Export. Serialize the pinned realization, transform description, grid version, and residual shift into the parcel attributes and the lineage store.
  6. Submit. Forward only verified geometry to credited-area accounting, and hand the lineage record to registry submission.

Provision the PROJ data directory as a first-class dependency of the container image, not an afterthought — pin the NTv2/NADCON grids you rely on and assert their presence at start-up, because the whole class of silent shift is ultimately an environment-provisioning failure that masquerades as a data-quality one. With the diagnostic measuring displacement before area is ever computed and the transform failing closed when the accurate grid is absent, a datum shift stops being a silent metres-scale corruption of credited area and becomes a logged, quarantined, auditable event.