Applying a Helmert 7-Parameter Transform in Python

The Helmert 7-parameter transform (also called the Bursa-Wolf similarity transform) maps geocentric coordinates from one datum to another using three translations, three rotations, and a single scale factor. The core equation is X' = T + (1 + s)·R·X, evaluated on Earth-Centred Earth-Fixed (ECEF) Cartesian coordinates — never on raw latitude and longitude. This page is an implementation reference within the CRS Normalization Workflows topic; read that first for environment setup and version pinning. Below, two independent implementations — a transparent numpy version and a pyproj helmert pipeline — are cross-validated on a control point so you can trust the parameters before pushing them into a production alignment job.

How the 7-Parameter Transform Works

A datum is defined by an ellipsoid and its position and orientation relative to the Earth’s centre of mass. Two datums — say a national realization and WGS84 — differ by a small rigid motion: a shift of origin, a tiny rotation about each axis, and a scale difference of a few parts-per-million. The 7-parameter transform captures exactly this rigid similarity. It cannot model local distortion, which is why survey authorities publish grid shifts for high-accuracy work.

The transform is defined on geocentric ECEF coordinates, so the full workflow has three stages: convert source geographic coordinates to ECEF using the source ellipsoid, apply the seven parameters, then convert the target ECEF back to geographic on the target ellipsoid. The rotations are small (typically under a few arc-seconds), so the exact rotation matrix is replaced by its small-angle linearisation.

$\begin{pmatrix} X' \\ Y' \\ Z' \end{pmatrix} = \begin{pmatrix} t_x \\ t_y \\ t_z \end{pmatrix} + (1 + s)\,\mathbf{R}\begin{pmatrix} X \\ Y \\ Z \end{pmatrix}$

Here (t_x, t_y, t_z) are the translations in metres, s is the scale expressed as a unitless factor (a value given in ppm is multiplied by $10^{-6}$), and $\mathbf{R}$ is the small-angle rotation matrix built from r_x, r_y, r_z in radians:

$\mathbf{R} = \begin{pmatrix} 1 & -r_z & r_y \\ r_z & 1 & -r_x \\ -r_y & r_x & 1 \end{pmatrix}$

This is the position_vector convention (EPSG method 1033). The coordinate_frame convention (EPSG method 1032) flips the sign of every off-diagonal rotation term — equivalently, it negates r_x, r_y, r_z. The magnitude of the rotations is identical; only their sign differs. A parameter set is meaningless without its convention label, and mixing them is the single most common source of a few-metre systematic error.

pyproj (through the PROJ engine) implements the transform as a +proj=helmert operation and handles the arc-second and ppm conversions internally. The numpy version below makes every step explicit so you can audit the sign convention and unit handling yourself, then confirms it matches PROJ to sub-millimetre agreement.

Production-Ready Script

The script converts a geographic control point on a source datum to WGS84 two ways: a hand-written numpy Helmert on ECEF coordinates, and a pyproj helmert pipeline. It then confirms the two agree, which validates both the parameters and the convention.

# numpy>=1.24.0, pyproj>=3.5.0, Python 3.9+
from __future__ import annotations

import numpy as np
from pyproj import Transformer

ARCSEC_TO_RAD = np.pi / 648000.0  # pi / (180 * 3600)


def geodetic_to_ecef(lon, lat, h, a: float, f: float) -> np.ndarray:
    """Geographic (deg, deg, m) -> geocentric ECEF (m) on the given ellipsoid."""
    lon = np.radians(np.asarray(lon, dtype=np.float64))
    lat = np.radians(np.asarray(lat, dtype=np.float64))
    h = np.asarray(h, dtype=np.float64)
    e2 = f * (2.0 - f)                      # first eccentricity squared
    n = a / np.sqrt(1.0 - e2 * np.sin(lat) ** 2)  # prime vertical radius
    x = (n + h) * np.cos(lat) * np.cos(lon)
    y = (n + h) * np.cos(lat) * np.sin(lon)
    z = (n * (1.0 - e2) + h) * np.sin(lat)
    return np.column_stack((x, y, z))


def ecef_to_geodetic(xyz: np.ndarray, a: float, f: float):
    """Geocentric ECEF (m) -> geographic (deg, deg, m); Bowring iteration."""
    xyz = np.atleast_2d(np.asarray(xyz, dtype=np.float64))
    x, y, z = xyz[:, 0], xyz[:, 1], xyz[:, 2]
    e2 = f * (2.0 - f)
    b = a * (1.0 - f)
    p = np.hypot(x, y)
    lat = np.arctan2(z, p * (1.0 - e2))     # initial guess
    for _ in range(6):                      # converges in <5 for terrestrial data
        n = a / np.sqrt(1.0 - e2 * np.sin(lat) ** 2)
        h = p / np.cos(lat) - n
        lat = np.arctan2(z, p * (1.0 - e2 * n / (n + h)))
    lon = np.arctan2(y, x)
    n = a / np.sqrt(1.0 - e2 * np.sin(lat) ** 2)
    h = p / np.cos(lat) - n
    return np.degrees(lon), np.degrees(lat), h


def helmert_7param(
    xyz: np.ndarray,
    tx: float, ty: float, tz: float,
    rx: float, ry: float, rz: float,   # arc-seconds
    s: float,                          # ppm
    convention: str = "position_vector",
) -> np.ndarray:
    """Apply the 7-parameter similarity transform to ECEF coordinates (N, 3)."""
    xyz = np.atleast_2d(np.asarray(xyz, dtype=np.float64))
    rx, ry, rz = (v * ARCSEC_TO_RAD for v in (rx, ry, rz))
    scale = 1.0 + s * 1e-6
    if convention == "position_vector":
        r = np.array([[1.0, -rz, ry], [rz, 1.0, -rx], [-ry, rx, 1.0]])
    elif convention == "coordinate_frame":
        r = np.array([[1.0, rz, -ry], [-rz, 1.0, rx], [ry, -rx, 1.0]])
    else:
        raise ValueError(f"Unknown convention: {convention!r}")
    t = np.array([tx, ty, tz])
    return t + scale * (xyz @ r.T)


if __name__ == "__main__":
    # Illustrative parameters: a national datum -> WGS84, position_vector.
    # Always take the values (and convention) from an authoritative registry.
    TX, TY, TZ = -446.448, 125.157, -542.060      # metres
    RX, RY, RZ = -0.1502, -0.2470, -0.8421        # arc-seconds
    S = 20.4894                                    # ppm
    # Source ellipsoid (Airy 1830) and target ellipsoid (WGS84 / GRS80-like).
    SRC_A, SRC_F = 6377563.396, 1.0 / 299.3249646
    WGS_A, WGS_F = 6378137.0, 1.0 / 298.257223563

    # A control point expressed on the source datum.
    lon0, lat0, h0 = -1.54700, 55.00000, 100.0

    # --- (a) numpy path: geographic -> ECEF -> Helmert -> geographic ---
    src_ecef = geodetic_to_ecef(lon0, lat0, h0, SRC_A, SRC_F)
    dst_ecef_np = helmert_7param(src_ecef, TX, TY, TZ, RX, RY, RZ, S,
                                 convention="position_vector")
    lon_np, lat_np, h_np = ecef_to_geodetic(dst_ecef_np, WGS_A, WGS_F)

    # --- (b) pyproj path: a helmert pipeline operating directly on ECEF ---
    # PROJ takes rotations in arc-seconds and scale in ppm as given.
    pipeline = (
        f"+proj=pipeline +step +proj=helmert "
        f"+x={TX} +y={TY} +z={TZ} +rx={RX} +ry={RY} +rz={RZ} +s={S} "
        f"+convention=position_vector"
    )
    ecef_tf = Transformer.from_pipeline(pipeline)
    dst_ecef_pj = np.column_stack(
        ecef_tf.transform(src_ecef[:, 0], src_ecef[:, 1], src_ecef[:, 2])
    )

    residual = np.linalg.norm(dst_ecef_np - dst_ecef_pj)
    print(f"numpy vs pyproj ECEF agreement: {residual * 1e3:.4f} mm")
    print(f"WGS84 control point: lon={lon_np[0]:.8f}  lat={lat_np[0]:.8f}")
    assert residual < 1e-3, "numpy and pyproj disagree — check convention/units"

Key implementation notes:

• geodetic_to_ecef and ecef_to_geodetic use each datum’s own ellipsoid — the source ellipsoid going in, the target ellipsoid coming out. Using WGS84 on both ends silently discards the ellipsoid difference and biases height.
• The numpy rotation matrix is written for both conventions in one place. Flip convention and the off-diagonal signs flip — the exact behaviour PROJ applies internally.
• The pyproj pipeline is deliberately ECEF-in, ECEF-out so it isolates the Helmert step for direct numerical comparison. In production you would prepend +step +proj=cart +ellps=<src> and append +step +inv +proj=cart +ellps=WGS84 to accept and return geographic coordinates.
• The assert at the end is the validation: sub-millimetre agreement on ECEF confirms the parameters, the arc-second conversion, and the convention are all consistent. If it fires, the convention is the first thing to check.

Compatibility Matrix

Fallback Strategies

1. Convention sign confusion (position_vector vs coordinate_frame)

A result that is offset by a consistent few metres — same direction for every point — almost always means the rotation signs are inverted. Registries and vendors disagree on which convention they publish, and some omit the label. If you only have one set of rotations, test both: negate r_x, r_y, r_z and re-run against a known control point. The correct convention drives the control-point residual toward zero; the wrong one leaves a stubborn systematic bias.

2. Rotations left in arc-seconds (or degrees) instead of radians

In the numpy path the rotation terms must be radians. An arc-second is about 4.85e-6 rad, so forgetting the ARCSEC_TO_RAD factor scales the rotation by ~206265 and throws points thousands of kilometres. PROJ avoids this by consuming arc-seconds directly — do not “pre-convert” values you pass to +rx/+ry/+rz, or you will double-apply the factor.

3. Applying the transform to geographic coordinates

The parameters are defined in ECEF metres. Passing degrees of latitude and longitude into helmert_7param produces nonsense. Always run the cartesian conversion first. When in doubt, sanity-check that your ECEF magnitudes are on the order of 6.4e6 metres before the Helmert step.

4. Ellipsoid mismatch on the return conversion

The inverse cartesian step must use the target ellipsoid. A common bug is reusing the source ellipsoid for both directions, which leaves a height error of tens of metres and a small horizontal shift. Keep the source and target ellipsoid parameters in separate, clearly named constants — as in the script — and pass the right pair to each conversion.

5. Prefer a datum grid shift for survey-grade accuracy

A single 7-parameter transform is a rigid model; it cannot absorb the local, non-linear distortion between datums that survey networks exhibit. Over a country it is typically accurate to a few decimetres. When you need centimetres, use a published NTv2 or geoid grid (for example OSTN15 for Great Britain, or the appropriate NADCON/NTv2 grid for North America). pyproj will select a grid-based operation automatically when the grids are installed — see the pyproj-based reprojection covered in Reprojecting CAD Coordinates with pyproj Transformer, and confirm grid availability with pyproj.datadir.

For the georeferencing metadata that supplies these parameters in a BIM context, see Reading IFC Georeferencing with ifcopenshell, and for the simpler 2D case where survey control replaces published datum parameters, see Converting CAD Local Coordinates to EPSG:4326.

FAQ

Related Pages

• CRS Normalization Workflows — parent topic: full pipeline from CRS detection through datum transformation and validation
• Reprojecting CAD Coordinates with pyproj Transformer — vectorized reprojection and automatic operation selection once the datum is resolved
• Converting CAD Local Coordinates to EPSG:4326 — the 2D survey-control counterpart when no published datum parameters exist
• Reading IFC Georeferencing with ifcopenshell — extracting the datum and map-conversion metadata that feed these parameters
• Coordinate Transformation & Spatial Alignment — domain overview covering datum alignment, unit conversion, and layer mapping