Propagating Covariance Through a Helmert Transformation

A seven-parameter transformation moves a coordinate and stretches its uncertainty, and the stretch has two independent sources: the input covariance rotated by the transformation, and the uncertainty of the seven parameters themselves. This guide, part of uncertainty propagation through transformation chains, implements both, and shows why the second term behaves differently from every other contribution in the budget — it is the same for every point transformed with those parameters, so it never averages away.

The Two Jacobians #

The transformed position is $\mathbf{y} = \mathbf{T} + s\,\mathbf{R}\,\mathbf{x}$, with translation vector $\mathbf{T}$, scale $s$ and rotation matrix $\mathbf{R}$. Differentiating with respect to the two things that carry uncertainty gives two Jacobians. With respect to the input coordinate,

Figure — two derivatives, two very different magnitudes.

which for the small rotations of a real parameter set is very nearly the identity — a scale of one plus a part per billion, a rotation of a few tens of milliarcseconds. The input covariance is therefore essentially rotated rather than inflated, and its determinant is preserved to within parts per billion. With respect to the seven parameters, evaluated at the point being transformed,

where $\mathbf{S}(\mathbf{x})$ is the skew-symmetric matrix of the position. That second Jacobian is the interesting one: its rotation block scales with the magnitude of the position vector, about 6.4 million metres, so a rotation uncertain by one milliarcsecond contributes about 31 mm at the Earth’s surface. Rotation parameter uncertainty dominates the parameter term for any realistic set.

Complete Runnable Implementation #

from __future__ import annotations

import numpy as np

MAS_TO_RAD = np.pi / (180.0 * 3600.0 * 1000.0)


def skew(v: np.ndarray) -> np.ndarray:
    """Skew-symmetric matrix S(v) with S(v) @ w == cross(v, w)."""
    return np.array([
        [0.0, -v[2], v[1]],
        [v[2], 0.0, -v[0]],
        [-v[1], v[0], 0.0],
    ], dtype=np.float64)


def helmert_jacobians(x: np.ndarray, scale: float,
                      rot_rad: np.ndarray) -> tuple[np.ndarray, np.ndarray]:
    """Jacobians of a 7-parameter transformation at the point `x`.

    Returns (J_x, J_p): the 3x3 derivative with respect to the input coordinate,
    and the 3x7 derivative with respect to (Tx, Ty, Tz, s, rx, ry, rz). Rotation
    columns are in RADIANS — convert a covariance quoted in mas before use.
    """
    if x.shape != (3,):
        raise ValueError("expected a (3,) Cartesian position in metres")
    r = np.eye(3) + skew(rot_rad)            # small-angle rotation
    j_x = scale * r
    j_p = np.empty((3, 7), dtype=np.float64)
    j_p[:, 0:3] = np.eye(3)                  # translations
    j_p[:, 3] = r @ x                        # scale (unitless -> metres)
    j_p[:, 4:7] = scale * skew(x).T          # rotations (radians -> metres)
    return j_x, j_p


def propagate_helmert(cov_x: np.ndarray, x: np.ndarray, scale: float,
                      rot_rad: np.ndarray,
                      cov_params: np.ndarray | None = None
                      ) -> tuple[np.ndarray, np.ndarray]:
    """Transformed covariance, split into its point and parameter halves.

    Returning the two terms separately is deliberate: the point term is
    independent between points, the parameter term is COMMON to every point
    transformed with this parameter set, and a relative statement between two
    points must exclude the second.
    """
    j_x, j_p = helmert_jacobians(x, scale, rot_rad)
    point_term = j_x @ cov_x @ j_x.T
    if cov_params is None:
        shared_term = np.zeros((3, 3), dtype=np.float64)
    else:
        if cov_params.shape != (7, 7):
            raise ValueError("parameter covariance must be (7, 7) in m, unitless, rad")
        shared_term = j_p @ cov_params @ j_p.T
    sym = lambda m: 0.5 * (m + m.T)
    return sym(point_term), sym(shared_term)


def params_cov_from_sigmas(sig_t_m: float, sig_scale_ppb: float,
                           sig_rot_mas: float) -> np.ndarray:
    """A diagonal parameter covariance from the sigmas a table usually quotes."""
    d = np.array([
        sig_t_m ** 2, sig_t_m ** 2, sig_t_m ** 2,
        (sig_scale_ppb * 1e-9) ** 2,
        (sig_rot_mas * MAS_TO_RAD) ** 2,
        (sig_rot_mas * MAS_TO_RAD) ** 2,
        (sig_rot_mas * MAS_TO_RAD) ** 2,
    ], dtype=np.float64)
    return np.diag(d)

Parameter Reference #

Worked Example #

A position with 12 mm isotropic uncertainty, transformed with parameters whose translations are known to 3 mm, scale to 0.3 ppb and rotations to 0.2 mas:

Figure — what each parameter uncertainty is worth at the Earth’s surface.

import numpy as np

x = np.array([-2694045.41, -4293642.37, 3857878.19])
cov_x = np.eye(3) * 0.012 ** 2
rot = np.array([26.78, -0.42, 10.93]) * MAS_TO_RAD
cov_p = params_cov_from_sigmas(sig_t_m=0.003, sig_scale_ppb=0.3, sig_rot_mas=0.2)

point, shared = propagate_helmert(cov_x, x, 1.0 + 0.37e-9, rot, cov_p)
print(np.round(np.sqrt(np.diag(point)) * 1000, 2))    # mm
print(np.round(np.sqrt(np.diag(shared)) * 1000, 2))   # mm
print(np.round(np.sqrt(np.diag(point + shared)) * 1000, 2))
# [12.   12.   12.  ]
# [ 5.32  4.42  5.9 ]
# [13.13 12.79 13.37]

The parameter term adds five to six millimetres, and almost all of it comes from the rotation columns — a 0.2 mas rotation uncertainty is about 6 mm at the Earth’s surface, while a 0.3 ppb scale uncertainty is under 2 mm. That ordering holds for essentially every published parameter set, and it is why a table that quotes translation uncertainties and omits rotation uncertainties is not usable for propagation.

Validation Check #

def check_parameter_jacobian(x: np.ndarray, scale: float, rot: np.ndarray) -> None:
    """Numerically differentiate one rotation column and compare with the analytic one."""
    _, j_p = helmert_jacobians(x, scale, rot)
    h = 1e-9
    perturbed = rot.copy()
    perturbed[0] += h
    r0 = (np.eye(3) + skew(rot)) @ x
    r1 = (np.eye(3) + skew(perturbed)) @ x
    numeric = scale * (r1 - r0) / h
    assert np.allclose(numeric, j_p[:, 4], rtol=1e-6), (
        "analytic rotation column disagrees with the numerical derivative"
    )

Numerical differentiation is the right check here because the sign conventions in a skew-symmetric matrix are exactly the kind of thing that reads correctly and is wrong; comparing against a finite difference of the actual transformation settles it without reference data.

Common Mistakes #

Rotation uncertainty left in milliarcseconds. The Jacobian’s rotation columns expect radians. Feeding a covariance in mas² produces a parameter term about 4 × 10¹¹ times too large — obvious when it appears, and easy to introduce when a covariance is assembled from a published table without conversion.

Parameter uncertainty added per point and then averaged. Because the parameter term is common to every point, computing a dataset-level uncertainty by averaging per-point values understates it by roughly the square root of the point count. Keep the shared term separate, as the function above does, and add it once at the level the statement is being made about.

Correlations between parameters discarded. Published parameter sets, especially Bursa-Wolf sets fitted over a regional network, have strongly correlated translations and rotations — correlations of 0.9 and above are normal. Using only the diagonal can overstate or understate the propagated uncertainty by tens of per cent depending on the sign of the correlation. If the full matrix is published, use it; the argument for the Molodensky-Badekas form in Bursa-Wolf vs Molodensky-Badekas transformations is precisely that it makes those correlations small.

A Note on Units in the Parameter Covariance #

The seven-by-seven parameter covariance is the one matrix in this whole area with mixed units along its diagonal: square metres for the three translations, unitless for the scale, and square radians for the three rotations. Nothing in its shape or in NumPy’s behaviour prevents a matrix assembled in the units a table publishes — metres, parts per million, milliarcseconds — from being used directly, and the result is a propagated uncertainty wrong by many orders of magnitude in the rotation block alone. Convert on construction, as params_cov_from_sigmas does, and keep the conversion in one place so there is exactly one line to check when a number looks implausible.

Figure — the one matrix in geodesy with three different units on its diagonal.

Frequently Asked Questions #

What if the parameter set is published without uncertainties? #

Then the parameter term cannot be computed, and the honest response is to say so in the uncertainty statement rather than to omit the term silently. A rough substitute is the operation’s declared accuracy treated as isotropic, which is what a registered coordinate operation supplies — coarser than a real covariance, but it is at least a stated figure rather than an implied zero.

Does the input covariance change much through the transformation? #

Barely. The Jacobian with respect to the input is a scale of one plus a part per billion times a rotation of a few tens of milliarcseconds, so the covariance is rotated by a negligible angle and scaled by a negligible factor. That is worth knowing because it means any large change in the propagated covariance came from the parameter term, not from the rotation of the input.

How does this interact with epoch propagation? #

They are separate stages and their uncertainties are independent, so they add. The epoch term is velocity uncertainty times the interval, which grows linearly with the interval while the parameter term is constant — over a long interval the epoch term can overtake everything else, as propagating coordinates between epochs with velocity grids shows.

Should the shared term appear in a per-point report? #

Yes, but labelled. A per-point uncertainty that includes the shared term is the right number for the point’s absolute position; the same point’s uncertainty relative to its neighbours excludes it. Publishing both, with the distinction stated, is what lets a downstream user compute an area or a bearing without either over- or understating its uncertainty.

Related #

• Uncertainty propagation through transformation chains — the parent topic
• Combining grid and parameter uncertainty in Python
• Helmert 7-parameter transformations in Python
• Deriving 7-parameter Helmert from control points
• Error distribution modeling in Python