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

Figure — a threshold is assembled from its parts, not chosen from a table.

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:

Figure — a threshold set too loose passes bad work; set too tight it rejects good work.

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.

Frequently Asked Questions #

Can I just use the tolerance from the survey specification? #

That is the ceiling, not the threshold. The specification says what the deliverable must achieve; your gate has to sit below it with enough margin to account for the uncertainty you cannot see at gate time. Setting the gate exactly at the specification means half of the work that passes is actually outside it once the unmodelled terms are included.

What if the propagated threshold is tighter than anything achievable? #

Then the uncertainty budget is telling you something true and unwelcome: the control, the transformation or the observations are not good enough for the specification. The answer is to improve one of them or to renegotiate the specification, and both are legitimate. Widening the gate quietly is the one response that is not.

Should thresholds be per-point or per-batch? #

Both, testing different things. A per-point gate catches individual failures; a batch statistic such as RMSE catches a systematic problem that leaves every point just inside its individual limit. A batch that passes on RMSE while several points fail individually is a batch with a real problem, and only the pair of tests reveals it.

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