Optimizing Transformation Tolerance Thresholds with Propagated Uncertainty

Optimizing a transformation tolerance threshold means replacing a fixed metre or centimetre cap with a per-point gate derived from propagated control-network uncertainty, so that every transformed coordinate is judged against its own statistically defensible bound rather than a single arbitrary constant — the survey-grade acceptance discipline established in tuning transformation thresholds for survey-grade conversion, the parent topic this page extends, and part of the broader Algorithmic Math & Geodetic Workflows reference. A static cutoff such as abs(residual) < 0.005 cannot account for the non-linear amplification of measurement uncertainty at network peripheries, in zones of high meridian convergence, or across chained operations; this page builds a self-contained evaluator that propagates each control point’s covariance through the transformation Jacobian and emits a deterministic PASS / WARN / FALLBACK_REQUIRED verdict suitable for an automated cadastral pipeline.

Figure — probabilistic tolerance gating: PASS, WARN, or fallback.

Concept: a threshold is a propagated 1-sigma bound, not a constant #

Geodetic transformations are lossy when they bridge a curved reference surface to a planar grid, and the measurement uncertainty attached to each control point does not survive that mapping unchanged. The standard first-order treatment carries a point’s variance-covariance matrix $\Sigma_i$ through the local Jacobian $J_i$ of the operation, giving the transformed covariance

The positional 1-sigma magnitude at that point is the square root of the trace of $\Sigma_i^{\text{out}}$, and the acceptance threshold is that magnitude scaled by a confidence coverage factor $k$:

Because $\tau_i$ tracks the geometric distortion encoded in $J_i$ and the survey quality encoded in $\Sigma_i$, it tightens automatically where control is dense and loosens where the network is weak — the spatially adaptive envelope that a scalar cutoff can never reproduce. The covariance matrices consumed here are the same ones produced by a least-squares adjustment for control networks, and the residual distributions the gate inspects feed directly into error distribution modeling in Python. When a residual exceeds the 95% bound, the transformation must be flagged, not forced.

Complete runnable implementation #

The module below propagates per-point covariance through an identity Jacobian (the correct local approximation for an affine transformation on a local grid or a low-order regional shift), computes the dynamic threshold, and routes the verdict deterministically. It is self-contained, uses np.float64 throughout, and raises on malformed input rather than degrading silently.

from __future__ import annotations

import logging
from dataclasses import dataclass
from typing import Any, Callable, Optional

import numpy as np
from numpy.linalg import LinAlgError, norm

logger = logging.getLogger("tolerance")


@dataclass(frozen=True)
class ToleranceVerdict:
    """Structured, audit-ready result of a tolerance evaluation."""
    status: str                 # "PASS" | "WARN" | "FALLBACK_REQUIRED"
    residuals: np.ndarray       # per-point radial residual (metre)
    thresholds: np.ndarray      # per-point dynamic tolerance (metre)
    max_violation_ratio: float  # max(residual / threshold) across the network
    metadata: dict[str, Any]


def propagate_point_std(
    point_covariance: np.ndarray,
    jacobian: np.ndarray,
    base_tolerance_m: float,
) -> float:
    """Propagate one point covariance through J: Sigma_out = J Sigma J^T.

    Returns the positional 1-sigma magnitude sqrt(tr(Sigma_out)), the scalar
    consistent with comparing against a radial residual norm. Falls back to the
    base tolerance on a singular/degenerate covariance (ISO 19111 operation
    accuracy must always resolve to a usable bound).
    """
    try:
        sigma_out = jacobian @ point_covariance @ jacobian.T
        return float(np.sqrt(max(float(np.trace(sigma_out)), 0.0)))
    except LinAlgError:
        logger.error("singular covariance for point; using base tolerance")
        return base_tolerance_m


def evaluate_tolerance(
    source_coords: np.ndarray,
    target_coords: np.ndarray,
    transform_fn: Callable[[np.ndarray], np.ndarray],
    control_covariance: Optional[np.ndarray] = None,
    base_tolerance_m: float = 0.005,
    confidence_k: float = 1.96,
    fallback_multiplier: float = 3.0,
) -> ToleranceVerdict:
    """Gate a coordinate transformation against propagated, per-point tolerance.

    control_covariance, when supplied, has shape (n, dim, dim): one covariance
    matrix per control point. When omitted, every point inherits a flat
    base_tolerance_m. A residual above its threshold yields WARN; a residual
    above fallback_multiplier * threshold yields FALLBACK_REQUIRED, the signal
    to degrade to a lower-order model and quarantine the dataset for review.
    """
    src = np.asarray(source_coords, dtype=np.float64)
    tgt = np.asarray(target_coords, dtype=np.float64)
    if src.shape != tgt.shape or src.ndim != 2:
        raise ValueError("source and target must be matching (n, dim) arrays")
    if src.size == 0:
        raise ValueError("coordinate arrays cannot be empty")

    transformed = np.asarray(transform_fn(src), dtype=np.float64)
    residuals = norm(transformed - tgt, axis=1)
    n_pts, dim = src.shape

    if control_covariance is not None:
        cov = np.asarray(control_covariance, dtype=np.float64)
        if cov.shape != (n_pts, dim, dim):
            raise ValueError(f"control_covariance must have shape {(n_pts, dim, dim)}")
        jac = np.eye(dim, dtype=np.float64)  # identity Jacobian for local affine / low-order shift
        stds = np.array(
            [propagate_point_std(cov[i], jac, base_tolerance_m) for i in range(n_pts)],
            dtype=np.float64,
        )
    else:
        stds = np.full(n_pts, base_tolerance_m, dtype=np.float64)

    thresholds = stds * confidence_k
    if not np.all(thresholds > 0.0):
        raise ValueError("dynamic thresholds must be strictly positive")

    ratios = residuals / thresholds
    max_ratio = float(np.max(ratios))
    if max_ratio <= 1.0:
        status = "PASS"
    elif max_ratio > fallback_multiplier:
        status = "FALLBACK_REQUIRED"
        logger.critical("residuals exceed fallback bound; degrade to affine-only routing")
    else:
        status = "WARN"

    return ToleranceVerdict(
        status=status,
        residuals=residuals,
        thresholds=thresholds,
        max_violation_ratio=max_ratio,
        metadata={
            "n_points": int(n_pts),
            "max_residual_m": round(float(np.max(residuals)), 6),
            "mean_residual_m": round(float(np.mean(residuals)), 6),
            "min_threshold_m": round(float(np.min(thresholds)), 6),
            "confidence_k": confidence_k,
        },
    )

The evaluator deliberately separates the numerical decision from the routing decision: when status == "FALLBACK_REQUIRED", the calling pipeline drops from the current model to a simpler one — typically a 3rd-order polynomial down to a 2D affine — to preserve topology while the dataset is escalated for manual surveyor review. That degradation is the same explicit, non-silent contract used for fallback routing when grid files are missing.

Parameter and return-value reference #

Minimal worked example #

Four control points on a local site are transformed by a trial operation that carries a 4 mm easting bias; each point has an isotropic 3 mm standard deviation, so the propagated 95% threshold is $1.96\sqrt{3}\times0.003 \approx 0.0102\,\text{m}$.

import numpy as np

src = np.array([
    [412300.10, 5274500.20, 142.30],
    [412800.45, 5274950.65, 138.90],
    [413100.00, 5275400.10, 151.05],
    [412550.75, 5275100.30, 144.60],
], dtype=np.float64)
target = src.copy()                       # control truth

def trial_transform(p: np.ndarray) -> np.ndarray:
    out = p.copy()
    out[:, 0] += 0.004                     # 4 mm easting bias under test
    return out

sigma = 0.0030                             # 3 mm per-axis control std
cov = np.tile(np.diag([sigma**2, sigma**2, sigma**2]), (src.shape[0], 1, 1))

verdict = evaluate_tolerance(src, target, trial_transform, control_covariance=cov)
print(verdict.status, verdict.metadata["max_residual_m"], round(verdict.max_violation_ratio, 4))
# PASS 0.004 0.3928

Validation check #

In a pipeline, the verdict must be asserted before the parameters are committed to a transformation registry, so a regression run fails loudly on any tolerance escalation rather than shipping a degraded dataset:

assert verdict.status == "PASS", f"tolerance gate rejected the operation: {verdict.metadata}"
assert verdict.max_violation_ratio <= 1.0, "at least one residual exceeded its propagated 95% bound"

Guarding against threshold drift #

A propagated threshold is only as trustworthy as the covariance feeding it, so the inputs themselves need version discipline. Store each $\Sigma_i$ alongside its point coordinates and realization epoch in a tracked repository; any change to control quality, epoch, or the national framework (for example a NAD83(2011) → NATRF2022 transition) then forces re-evaluation instead of silently shifting the envelope. Export max_residual_m and min_threshold_m to monitoring so that a sudden contraction of thresholds — usually a projection-zone boundary artifact or a misconfigured datum shift — surfaces before it quarantines a production batch. For chained operations such as ellipsoid → ECEF → projected grid, propagate covariance at each stage following the ISO 19111 coordinate-operation model rather than re-deriving a single end-to-end bound, and the underlying conversion geometry is covered in geodetic conversion math from ellipsoid to Cartesian.

Common mistakes #

1. Comparing a radial residual against a per-axis sigma. The residual returned here is the Euclidean norm over all axes, so the threshold must be the positional magnitude $\sqrt{\operatorname{tr}(\Sigma^{\text{out}})}$, not a single diagonal entry. Gating a 3D residual against one axis’s $\sigma$ understates the bound by roughly $\sqrt{3}$ and passes datasets that should warn.
2. Feeding standard deviations where covariances are expected. control_covariance holds variances and cross-terms in metre² — passing a (n, dim) array of metres silently mis-shapes the input or, worse, treats millimetre values as variances. Build each block with np.diag([sx**2, sy**2, sz**2]) and only then add off-diagonal correlation.
3. Leaving the identity Jacobian in place for a high-distortion projection. The np.eye(dim) approximation is valid for a local affine or low-order regional shift; for a strongly non-linear projection you must substitute the analytical Jacobian per point, or the propagated threshold will not track the real geometric amplification — the same modelling distinction that separates affine fits from polynomial shift algorithms for regional adjustments.

Related #

• Tuning Transformation Thresholds for Survey-Grade Conversion — parent topic: the statute-anchored acceptance model this evaluator enforces
• Algorithmic Math & Geodetic Workflows — the reference architecture for deterministic geodetic computation
• Least-Squares Adjustment for Control Networks — produces the covariance matrices this gate propagates
• Error Distribution Modeling in Python — turns the gated residuals into a full uncertainty budget
• Validating Datum Alignment with Control Points — independent check-point comparison upstream of threshold tuning