Uncertainty Propagation Through Transformation Chains

A coordinate that arrives without an uncertainty is an assertion; a coordinate that arrives with one is evidence. This topic, part of algorithmic math and geodetic workflows, is about carrying the second kind through a pipeline: taking the covariance that came out of an adjustment, transforming it correctly at every stage, adding what each stage contributes, and arriving at a number a surveyor can defend. The individual stages — the Helmert transformation, the grid shift, the projection — are covered in their own topics; what this one adds is the arithmetic that connects them and the discipline that keeps the result honest.

The reason it needs a topic of its own is that uncertainty does not travel the way coordinates do. A coordinate passes through a transformation; a covariance matrix passes through the derivative of that transformation, and the two are only the same thing when the transformation is a translation. Everything below follows from that one fact.

The Propagation Law #

For a function $\mathbf{y} = f(\mathbf{x})$ with Jacobian $\mathbf{J} = \partial f/\partial \mathbf{x}$, the covariance of the output is

Figure — a coordinate passes through the transformation; its covariance passes through the derivative.

where $\boldsymbol{\Sigma}_{\text{stage}}$ is whatever uncertainty the stage itself introduces — the declared accuracy of a coordinate operation, the interpolation error of a grid, the model uncertainty of a geoid. The first term propagates what came in; the second adds what happens here. Omitting either is a common and asymmetric mistake: dropping the first understates the input, dropping the second understates the operation, and both err in the direction that makes a deliverable look better than it is.

Three properties of this law shape the implementation. It is linear, so it is exact for a linear transformation and a first-order approximation otherwise — which is fine for geodetic transformations, whose non-linearity over the range of a coordinate uncertainty is negligible. It requires the Jacobian in the right frame: a covariance expressed in local east-north-up and a Jacobian expressed in Earth-centred Cartesian cannot be multiplied together, and nothing in the shapes prevents you from trying. And it preserves symmetry and positive semi-definiteness exactly in theory and approximately in floating point, which makes both properties useful as assertions.

A Chain, Stage by Stage #

Figure — the conversions add nothing; the physical stages are where uncertainty enters.

The pattern in that table is the useful part: the pure coordinate conversions add nothing and only rotate the existing covariance, while the physical stages — the datum shift, the grid, the epoch — are where uncertainty enters. It is also why the Helmert row is the awkward one: its parameter uncertainty is shared by every point transformed with the same parameters, so it is correlated across the dataset and does not average away over many points. Treating it as independent per-point noise understates the uncertainty of any aggregate, which is the subject of propagating covariance through a Helmert transformation.

Production Implementation #

from __future__ import annotations

from dataclasses import dataclass

import numpy as np


@dataclass(frozen=True)
class Uncertain:
    """A position and its covariance, in one frame, with a stated unit."""

    x: np.ndarray            # (3,) float64
    cov: np.ndarray          # (3, 3) float64, same units squared
    frame: str
    units: str = "m"

    def __post_init__(self) -> None:
        if self.x.shape != (3,) or self.cov.shape != (3, 3):
            raise ValueError("expected a (3,) position and a (3, 3) covariance")
        if not np.allclose(self.cov, self.cov.T, atol=1e-15):
            raise ValueError("covariance is not symmetric")
        if np.min(np.linalg.eigvalsh(self.cov)) < -1e-15:
            raise ValueError("covariance is not positive semi-definite")

    def through(self, jac: np.ndarray, stage_cov: np.ndarray | None,
                new_frame: str | None = None) -> "Uncertain":
        """Propagate through one stage: J S J^T plus what the stage contributes."""
        if jac.shape != (3, 3):
            raise ValueError("expected a (3, 3) Jacobian")
        cov = jac @ self.cov @ jac.T
        if stage_cov is not None:
            if stage_cov.shape != (3, 3):
                raise ValueError("stage covariance must be (3, 3)")
            cov = cov + stage_cov
        cov = 0.5 * (cov + cov.T)          # re-symmetrise against float drift
        return Uncertain(self.x, cov, new_frame or self.frame, self.units)

    def sigma(self) -> np.ndarray:
        """Per-axis one-sigma values — the numbers a report usually quotes."""
        return np.sqrt(np.diag(self.cov))

    def horizontal_95(self) -> float:
        """95% horizontal confidence radius, from the 2-D chi-square scaling.

        2.4477 is the two-dimensional 95% factor. Using the one-dimensional 1.96
        here produces a region that covers about 86%, which is the single most
        common reporting error in this area.
        """
        eig = np.linalg.eigvalsh(self.cov[:2, :2])
        return float(2.4477 * np.sqrt(max(eig)))


def isotropic(sigma: float) -> np.ndarray:
    """A diagonal covariance for a stage with one quoted accuracy figure."""
    return np.eye(3, dtype=np.float64) * (sigma ** 2)

The through method is deliberately narrow: one stage, one Jacobian, one optional stage covariance, and a re-symmetrisation at the end. Chains built from it read as a sequence of stages, which is also how they must appear in the audit record.

Worked Example: a Four-Stage Chain #

A GNSS-derived position with a 12 mm horizontal and 25 mm vertical one-sigma, propagated five years on a velocity known to ±1 mm/yr, shifted by a grid with a 4 mm interpolation term, and transformed by a Helmert operation with a declared accuracy of 15 mm:

import numpy as np

start = Uncertain(
    x=np.array([-2694045.41, -4293642.37, 3857878.19]),
    cov=np.diag([0.012 ** 2, 0.012 ** 2, 0.025 ** 2]),
    frame="ITRF2014",
)
I = np.eye(3)
after_epoch = start.through(I, isotropic(0.001 * 5.0))       # velocity x interval
after_grid = after_epoch.through(I, isotropic(0.004))        # interpolation
final = after_grid.through(I, isotropic(0.015), new_frame="NAD83(2011)")

print(np.round(final.sigma() * 1000, 2))       # mm
print(f"95% horizontal radius {final.horizontal_95() * 1000:.1f} mm")
# [19.92 19.92 29.77]
# 95% horizontal radius 48.8 mm

Twelve millimetres of observation uncertainty became twenty, and the 95 per cent horizontal radius is 49 mm — over twice the one-sigma figure a report might otherwise have quoted. Neither number is pessimistic; they are what the stated inputs imply.

Correlated and Independent Contributions #

The quadrature sum that combines uncertainty terms assumes they are independent, and two of the terms in a typical chain are not. Helmert parameter uncertainty is shared by every point transformed with that parameter set: if the parameters are biased by five millimetres in a direction, every point moves five millimetres in that direction, and the mean residual over a thousand points is still five millimetres rather than five divided by the square root of a thousand. Grid interpolation error is spatially correlated: neighbouring points sit in the same cell or in adjacent ones, so their interpolation errors are nearly identical, and a parcel boundary can be displaced coherently while every individual corner is inside tolerance.

Figure — averaging shrinks the independent part and leaves the correlated part exactly where it was.

The practical consequence is a distinction between two questions that get conflated. “How uncertain is this point?” is answered by the full quadrature sum, correlated terms included. “How uncertain is the relative geometry between two nearby points?” is answered by the independent terms only, because the correlated terms move both points together and cancel in the difference. A parcel’s area and its bearings depend on the second question; its position in the national frame depends on the first. Reporting one number for both is how a deliverable ends up simultaneously overstating the uncertainty of a boundary dimension and understating the uncertainty of the position.

def relative_uncertainty(point_cov: np.ndarray, shared_cov: np.ndarray) -> np.ndarray:
    """Covariance of the DIFFERENCE between two nearby points.

    Terms common to both points cancel in the difference, so a relative statement
    excludes them. Passing the full point covariance here — including the shared
    parameter term — overstates the uncertainty of a boundary dimension.
    """
    independent = point_cov - shared_cov
    if np.min(np.linalg.eigvalsh(independent)) < -1e-15:
        raise ValueError(
            "shared term exceeds the total: the two covariances are inconsistent"
        )
    return 2.0 * independent          # two independent points, one difference

The factor of two is the part that is easy to lose: the difference of two independent quantities has the sum of their variances, so a relative statement between two points of equal quality is larger by a factor of the square root of two than the independent part of either one — not smaller, as the cancellation of the shared term might suggest at a glance.

Validation: Three Properties Worth Asserting #

Uncertainty code fails quietly, because a wrong covariance is still a plausible covariance. Three assertions catch most of it, and all three are exact consequences of the propagation law rather than tolerances chosen by taste.

The propagated covariance must remain symmetric: any asymmetry beyond floating-point noise means a Jacobian was transposed in one place and not another. It must remain positive semi-definite: a negative eigenvalue means either a subtraction that should have been an addition, or a stage covariance supplied with the wrong sign. And propagating through a transformation and then through its inverse must return the original covariance, which tests the Jacobian and its inverse against each other without needing a reference implementation of either.

def check_chain(start: Uncertain, jac: np.ndarray) -> None:
    """Round-trip a covariance through a Jacobian and its inverse."""
    there = start.through(jac, None)
    back = there.through(np.linalg.inv(jac), None)
    assert np.allclose(back.cov, start.cov, atol=1e-18), (
        "covariance did not survive the round trip — check the Jacobian orientation"
    )
    assert np.allclose(there.cov, there.cov.T, atol=1e-18), "asymmetry introduced"
    assert np.min(np.linalg.eigvalsh(there.cov)) >= -1e-18, "lost definiteness"

The tolerance of 1e-18 looks aggressive until you notice the units: these are variances in square metres, so 1e-18 m² is a nanometre of standard deviation. A round trip that closes to a nanometre has nothing wrong with it; one that closes to a millimetre has a real error in the Jacobian, and the difference between those two outcomes is exactly what a loose tolerance would hide.

Where the Jacobians Come From #

Three of the stages in the table above need a real Jacobian rather than an identity, and all three have closed forms worth writing down once rather than differentiating numerically per point.

Geodetic to Cartesian. The derivative of the forward conversion with respect to latitude, longitude and ellipsoidal height involves only the two radii of curvature and trigonometric terms already computed for the conversion itself. Its inverse serves the reverse stage, which is why the round-trip test on the covariance is a genuine check on both.

The Helmert step. With respect to the input coordinate the Jacobian is scale times the rotation matrix, which for a real parameter set is within parts per billion of the identity — so the input covariance is rotated, not inflated. With respect to the seven parameters it is a three-by-seven matrix whose rotation columns scale with the position vector, and that is where the parameter uncertainty enters at millimetre magnitudes.

The projection. For a conformal projection the local Jacobian is a rotation by the grid convergence and a uniform scaling by the point scale factor. That has a pleasant consequence: a conformal projection rotates an error ellipse and changes its size, but cannot change its shape — the ratio of the semi-axes is preserved. An error ellipse that changes eccentricity through a projection step therefore indicates a bug rather than a property of the map.

def projection_jacobian(point_scale: float, convergence_deg: float) -> np.ndarray:
    """Local 2-D Jacobian of a conformal projection: rotate, then scale."""
    g = np.radians(convergence_deg)
    return point_scale * np.array([[np.cos(g), -np.sin(g)],
                                   [np.sin(g), np.cos(g)]], dtype=np.float64)

Building each stage’s Jacobian from the quantities the stage already computes — the radii, the rotation matrix, the scale factor and convergence — keeps the propagation consistent with the transformation it describes. A Jacobian derived independently, from a textbook rather than from the code that ran, will eventually disagree with it, and the disagreement will be invisible because both produce plausible covariances.

Compliance and Reporting #

An uncertainty statement is only reviewable if its inputs are visible. The record should carry each stage’s contribution separately rather than a single total: the input covariance and its source, the stage covariance for each operation and where its figure came from, and the confidence level and dimensionality of any quoted region. A single number labelled “accuracy: 0.05 m” with no dimensionality, no confidence level and no breakdown is the statement this topic exists to replace — and it is the one that cannot be defended when a neighbouring survey disagrees. The reporting conventions are worked through in reporting 95 per cent confidence regions for cadastral points.

Failure Modes #

• Frame mismatch between covariance and Jacobian. A covariance in local east-north-up multiplied by a Cartesian Jacobian produces a matrix of the right shape and the wrong meaning. Carry the frame with the covariance, as the type above does.
• Stage uncertainty omitted. Propagating only the input covariance produces a final uncertainty smaller than the declared accuracy of the operation that was applied — which is self-evidently impossible and routinely published.
• Parameter uncertainty treated as independent per point. Helmert parameter error is common to every point using those parameters; averaging over a million points does not reduce it.
• One-dimensional scale factors used on a two-dimensional region. 1.96 for a 95 per cent ellipse under-covers by roughly nine percentage points.
• Negative eigenvalues ignored. Accumulated floating-point asymmetry can make a covariance slightly indefinite; re-symmetrising each stage keeps it clean, and asserting positive semi-definiteness catches a real bug when it appears.

What a Complete Statement Looks Like #

Pulling the pieces together, a defensible per-point uncertainty statement in a cadastral deliverable carries six things: the horizontal semi-axes and their azimuth, the vertical bound, the confidence level, the dimensionality each figure applies to, the breakdown by contributing term, and whether the figure is absolute or relative to neighbouring points. Five of the six are one number each; the sixth is a short table. Together they let a reviewer recompute the statement rather than trust it, which is the only property that matters when two surveys of the same boundary disagree and both are defended by their own numbers.

Frequently Asked Questions #

Do I need the full covariance, or are per-axis sigmas enough? #

Per-axis sigmas are enough only while nothing rotates. The moment a transformation mixes axes — any rotation, any projection — the correlations matter, and discarding the off-diagonal terms understates or overstates the transformed uncertainty depending on the direction. Carrying a three-by-three matrix costs nine numbers instead of three and removes the question.

Where do I get the stage uncertainty for a coordinate operation? #

From its declared accuracy, which every registered operation carries and which a transformer object will report. That figure is a one-sigma-equivalent horizontal accuracy for the operation as a whole; treating it as isotropic is a simplification, and a reasonable one unless the operation’s documentation says otherwise. An operation with no declared accuracy is a flag in itself.

Should uncertainty be propagated per point or once for the dataset? #

Per point, because the input uncertainties differ per point and the stage contributions may too — a grid’s interpolation error varies across the surface. A dataset-level summary is fine as a headline, but if it is the only number published, every point inherits the average and the weakest points are misrepresented.

Does this apply to a purely two-dimensional workflow? #

Yes, with two-by-two matrices and the same law. The one thing that changes is the confidence scaling, which is dimension-dependent: the 95 per cent factor is 1.96 in one dimension, 2.4477 in two and 2.7955 in three. Using the wrong one is the most common arithmetic error in the whole area.

Related #

• Algorithmic math and geodetic workflows — the parent section
• Propagating covariance through a Helmert transformation
• Combining grid and parameter uncertainty in Python
• Error distribution modeling in Python
• Tuning transformation thresholds for survey-grade work