CoordTransform Cadastral & High-Precision Python

Deriving 7-Parameter Helmert from Control Points

Estimating the seven Helmert parameters from points whose coordinates are known in both datums is the least-squares step that produces the constants applied in Helmert 7-parameter transformations in Python, the parent topic under Core Transformation Fundamentals & Standards. Given at least three common control points, we linearise the transformation, assemble a design matrix, solve the normal equations for [tx,ty,tz,rx,ry,rz,ds][t_x, t_y, t_z, r_x, r_y, r_z, d_s], and report the residuals that certify the fit.

Least-squares estimation of Helmert parameters Control points known in both the source and target datum feed the linearised design matrix A. The design matrix and the observation vector form the normal equations A-transpose-A. Solving them yields the seven parameters, and back-substitution yields the residual vector used to grade the fit. Common control points (both datums) Design matrix A 3n × 7, linearised Normal equations AᵀA x = Aᵀℓ 7 params + residuals

Figure — control points known in two datums drive the design matrix, the normal equations, and the seven-parameter solution.

The Linearised Observation Equations

Because the rotations and scale are tiny, the position vector transformation linearises cleanly around zero rotation and unit scale. For a single control point with source coordinates (Xs,Ys,Zs)(X_s, Y_s, Z_s) and target (Xt,Yt,Zt)(X_t, Y_t, Z_t), keeping only first-order terms gives three observation equations, one per axis:

Xt−Xs=tx+Zs ry−Ys rz+Xs dsYt−Ys=ty−Zs rx+Xs rz+Ys dsZt−Zs=tz+Ys rx−Xs ry+Zs ds \begin{aligned} X_t - X_s &= t_x + Z_s\,r_y - Y_s\,r_z + X_s\,d_s \\ Y_t - Y_s &= t_y - Z_s\,r_x + X_s\,r_z + Y_s\,d_s \\ Z_t - Z_s &= t_z + Y_s\,r_x - X_s\,r_y + Z_s\,d_s \end{aligned}

Here rx,ry,rzr_x, r_y, r_z are in radians and dsd_s is the dimensionless scale ratio. Each control point contributes three rows to the design matrix A\mathbf{A}, whose columns are the partial derivatives with respect to the seven unknowns x=[tx,ty,tz,rx,ry,rz,ds]⊤\mathbf{x} = [t_x, t_y, t_z, r_x, r_y, r_z, d_s]^\top. With nn points the system is 3n3n equations in 77 unknowns, so three points give the minimum 9>79 > 7 and any surplus is redundancy that the least-squares solution averages. The observation vector ℓ\boldsymbol{\ell} holds the coordinate differences Xt−XsX_t - X_s and so on. Solving the normal equations A⊤A x=A⊤ℓ\mathbf{A}^\top\mathbf{A}\,\mathbf{x} = \mathbf{A}^\top\boldsymbol{\ell} minimises the sum of squared residuals, the same estimator used for a full network in least-squares adjustment for control networks.

Complete Runnable Implementation

The estimator below builds the design matrix row by row, solves the normal equations at float64, and returns both the parameters in published units and the residual vector. Coordinates are centred on their mean before the solve to improve conditioning, and the translation is corrected back to the original origin afterwards.

from __future__ import annotations
from dataclasses import dataclass
import math

import numpy as np

ARCSEC_TO_RAD: float = math.pi / (180.0 * 3600.0)


@dataclass(frozen=True)
class HelmertEstimate:
    """Estimated parameters in published units plus fit diagnostics."""
    tx: float           # metres
    ty: float           # metres
    tz: float           # metres
    rx: float           # arc-seconds
    ry: float           # arc-seconds
    rz: float           # arc-seconds
    ds: float           # ppm
    residuals: np.ndarray   # 3n residual components, metres
    rms_m: float            # root-mean-square residual, metres


def _design_rows(x: float, y: float, z: float) -> np.ndarray:
    """Three linearised rows for one point; columns tx,ty,tz,rx,ry,rz,ds."""
    return np.array([
        [1.0, 0.0, 0.0,  0.0,    z,   -y,   x],
        [0.0, 1.0, 0.0,   -z,  0.0,    x,   y],
        [0.0, 0.0, 1.0,    y,   -x,  0.0,   z],
    ], dtype=np.float64)


def estimate_helmert(src_pts: np.ndarray, dst_pts: np.ndarray) -> HelmertEstimate:
    """Least-squares 7-parameter Helmert from >=3 common control points.

    src_pts, dst_pts: (n, 3) geocentric ECEF arrays in metres, same order.
    Returns position-vector-convention parameters and residual diagnostics."""
    src = np.asarray(src_pts, dtype=np.float64)
    dst = np.asarray(dst_pts, dtype=np.float64)
    n = src.shape[0]
    if n < 3:
        raise ValueError("need at least 3 common control points")
    if src.shape != dst.shape or src.shape[1] != 3:
        raise ValueError("src and dst must be matching (n, 3) arrays")

    # Centre on the source mean so the rotation/scale columns are not
    # dominated by the ~6.4e6 m geocentric radius (conditioning).
    centroid = src.mean(axis=0)
    src_c = src - centroid
    dst_c = dst - centroid

    rows = [_design_rows(*p) for p in src_c]
    A = np.vstack(rows)                       # (3n, 7)
    ell = (dst_c - src_c).reshape(-1)          # (3n,)

    normal = A.T @ A                           # AᵀA, 7x7
    x = np.linalg.solve(normal, A.T @ ell)     # solve normal equations
    residuals = A @ x - ell                    # v = Ax - ℓ
    rms = float(np.sqrt(np.mean(residuals ** 2)))

    tx, ty, tz, rx, ry, rz, ds = x
    # Correct translation from the centred origin back to the geocentre.
    scale = 1.0 + ds
    rot = np.array([[1.0, -rz, ry],
                    [rz, 1.0, -rx],
                    [-ry, rx, 1.0]], dtype=np.float64)
    t_geo = np.array([tx, ty, tz]) + (np.eye(3) - scale * rot) @ centroid

    return HelmertEstimate(
        tx=float(t_geo[0]), ty=float(t_geo[1]), tz=float(t_geo[2]),
        rx=rx / ARCSEC_TO_RAD, ry=ry / ARCSEC_TO_RAD, rz=rz / ARCSEC_TO_RAD,
        ds=ds * 1.0e6, residuals=residuals, rms_m=rms,
    )

Parameter Reference

Name Type Units Valid range Notes
src_pts np.ndarray metres (n, 3), n ≥ 3 source-datum ECEF of the control points
dst_pts np.ndarray metres (n, 3), same order target-datum ECEF of the same points
tx, ty, tz float metres solved translations at the geocentric origin
rx, ry, rz float arc-seconds solved rotations, position-vector convention
ds float ppm solved scale difference
residuals np.ndarray metres (3n,) per-component misclosures after the fit
rms_m float metres ≥ 0 root-mean-square residual, the headline quality figure

Worked Example

Generate five control points across Britain, carry them into the target datum with a known parameter set plus a few millimetres of noise, then recover the parameters and inspect the residuals.

def geodetic_to_ecef(lat: float, lon: float, h: float) -> np.ndarray:
    a, f = 6378137.0, 1.0 / 298.257223563
    e2 = f * (2.0 - f)
    p, q = math.radians(lat), math.radians(lon)
    nrad = a / math.sqrt(1.0 - e2 * math.sin(p) ** 2)
    return np.array([(nrad + h) * math.cos(p) * math.cos(q),
                     (nrad + h) * math.cos(p) * math.sin(q),
                     (nrad * (1.0 - e2) + h) * math.sin(p)], dtype=np.float64)


src = np.array([geodetic_to_ecef(50.1, -5.0, 50), geodetic_to_ecef(54.0, -1.0, 120),
                geodetic_to_ecef(57.5, -4.0, 300), geodetic_to_ecef(52.0, 0.5, 30),
                geodetic_to_ecef(53.3, -2.9, 80)])
rot = np.array([[1, 4.08e-7, -1.198e-6], [-4.08e-7, 1, -7.28e-7],
                [1.198e-6, 7.28e-7, 1]])                      # ~ -0.84,-0.15,-0.25″
dst = (1 + 20.4894e-6) * (src @ rot.T) + np.array([-446.448, 125.157, -542.060])
dst = dst + np.random.default_rng(42).normal(0, 0.005, dst.shape)  # 5 mm noise

est = estimate_helmert(src, dst)
print(round(est.tx, 3), round(est.ds, 4), round(est.rms_m, 4))
# -446.419 20.4868 0.0042

The recovered translation of about −446.4-446.4 m and scale of about 20.4920.49 ppm match the values that generated the data, and the residual RMS of roughly 44 mm reflects the injected noise averaged across fifteen observation components and seven estimated unknowns.

Validation Check

Gate the estimate on its residual RMS before the parameters are trusted for production transformations.

assert est.rms_m < 0.050, "fit residual exceeds 50 mm — inspect control quality"
assert abs(est.tx - (-446.448)) < 0.1, "translation drifted from expected"
assert np.max(np.abs(est.residuals)) < 0.020, "a control point carries a gross error"

Common Mistakes

Too few points, or points that are nearly colinear
Three points give the bare minimum of nine equations for seven unknowns, leaving almost no redundancy to detect a blunder. Worse, if the control lies along a line or on a small patch, the rotation and scale columns become nearly linearly dependent and AᵀA is ill-conditioned, so tiny observation noise swings the parameters wildly. Use well-distributed points spanning the working area and keep several redundant ones.
Mixing rotation conventions between estimation and application
This estimator returns position-vector-convention rotations. If you later feed them into a routine that expects the coordinate-frame convention, the sign of rx, ry, rz is wrong and the transformation carries a few metres of rotational error. Record the convention alongside the estimate and match it at application time.
Not centring the coordinates before the solve
The rotation and scale columns of the design matrix are the raw geocentric coordinates, on the order of six million metres, while the translation columns are ones. That six-order-of-magnitude spread makes AᵀA badly conditioned and amplifies rounding. Centring on the control centroid shrinks those columns and stabilises the solution — the same reasoning that motivates the Molodensky-Badekas form covered in the companion comparison page.