ITRF to NAD83(2011) Transformation in Python

Transforming between a global frame and a plate-fixed one is a fourteen-parameter operation, not a seven-parameter one, and treating it as the latter is a systematic error that grows by roughly two centimetres a year. This guide, part of time-dependent transformations and plate motion, implements the time-dependent Helmert transformation between ITRF and NAD83(2011): seven parameters at a reference epoch, seven rates of change, and an explicit coordinate epoch on both sides. It must satisfy the same centimetre-level tolerance as any other cadastral operation, and the whole margin is spent if the rates are dropped.

Why Fourteen Parameters #

The static seven-parameter model described in Helmert 7-parameter transformations in Python relates two frames that are fixed relative to one another. ITRF and NAD83 are not: NAD83 is tied to the North American plate and rotates with it, while ITRF is global, so the relationship between them changes continuously. Published parameter sets therefore come in two halves — values at a reference epoch $t_r$, and rates of change per year:

Figure — seven values, seven rates, and the one term that dominates.

applied to each of the three translations, three rotations and the scale. The dominant terms are the rotation rates: they express the plate rotation itself, and at a mid-latitude North American site they amount to roughly 15–20 mm of horizontal motion per year. A transformation evaluated at the reference epoch and then reused for a decade is out by 150–200 mm.

The second requirement is the coordinate epoch of the input. The fourteen parameters relate the two frames at an instant; they do not move a coordinate through time within a frame. A complete conversion is therefore two operations in sequence — propagate to the epoch at which the transformation is being evaluated, then transform frames — each recorded separately, as the parent topic sets out.

Complete Runnable Implementation #

from __future__ import annotations

from dataclasses import dataclass

import numpy as np

MAS_TO_RAD = np.pi / (180.0 * 3600.0 * 1000.0)     # milliarcseconds -> radians
PPB = 1.0e-9                                        # parts per billion -> unitless


@dataclass(frozen=True)
class TimeDependentHelmert:
    """A 14-parameter frame transformation: 7 values plus 7 rates.

    Translations in metres and metres/year, rotations in mas and mas/year, scale
    in ppb and ppb/year — the units every published table uses. `ref_epoch` is the
    decimal year the values are stated at.
    """

    tx: float; ty: float; tz: float                 # metres
    rx: float; ry: float; rz: float                 # milliarcseconds
    ds: float                                       # ppb
    dtx: float; dty: float; dtz: float              # metres / year
    drx: float; dry: float; drz: float              # mas / year
    dds: float                                      # ppb / year
    ref_epoch: float

    def at(self, epoch: float) -> tuple[np.ndarray, np.ndarray, float]:
        """Parameters evaluated at `epoch` -> (translation_m, rotation_rad, scale)."""
        dt = epoch - self.ref_epoch
        t = np.array([self.tx + self.dtx * dt,
                      self.ty + self.dty * dt,
                      self.tz + self.dtz * dt], dtype=np.float64)
        r = np.array([self.rx + self.drx * dt,
                      self.ry + self.dry * dt,
                      self.rz + self.drz * dt], dtype=np.float64) * MAS_TO_RAD
        s = 1.0 + (self.ds + self.dds * dt) * PPB
        return t, r, s

    def apply(self, xyz: np.ndarray, epoch: float) -> np.ndarray:
        """Transform a Cartesian position at `epoch`, coordinate-frame convention.

        The coordinate-frame (a.k.a. "position vector reversed") convention is the
        one EPSG uses for the ITRF-to-NAD83 parameter sets; flipping the sign of the
        three rotations gives the position-vector convention instead.
        """
        if xyz.shape != (3,) or xyz.dtype != np.float64:
            raise ValueError("xyz must be a (3,) float64 array of metres")
        t, r, s = self.at(epoch)
        rot = np.array([
            [1.0, r[2], -r[1]],
            [-r[2], 1.0, r[0]],
            [r[1], -r[0], 1.0],
        ], dtype=np.float64)                        # small-angle approximation
        return t + s * (rot @ xyz)

Parameter Reference #

Worked Example #

A monument with ITRF2014 Cartesian coordinates at epoch 2020.00, transformed with a representative ITRF2014 → NAD83(2011) parameter set:

Figure — error from evaluating the parameters at their reference epoch and reusing them.

import numpy as np

itrf_to_nad83 = TimeDependentHelmert(
    tx=1.00530, ty=-1.90921, tz=-0.54157,
    rx=26.78138, ry=-0.42027, rz=10.93206, ds=0.36891,
    dtx=0.00079, dty=-0.00060, dtz=-0.00134,
    drx=0.06667, dry=-0.75744, drz=-0.05133, dds=-0.10201,
    ref_epoch=2010.00,
)

xyz_itrf = np.array([-2_694_045.4131, -4_293_642.3672, 3_857_878.1948])
xyz_nad83 = itrf_to_nad83.apply(xyz_itrf, epoch=2020.00)
print(np.round(xyz_nad83, 4))
# [-2694044.6272 -4293643.5806  3857877.6165]

# The same transformation evaluated at the reference epoch instead:
xyz_wrong = itrf_to_nad83.apply(xyz_itrf, epoch=2010.00)
print(f"{np.linalg.norm(xyz_nad83 - xyz_wrong):.4f} m of drift over 10 years")
# 0.2038 m of drift over 10 years

Two hundred millimetres over ten years, from a single omitted argument, on a job whose tolerance is twenty. That is what the rates are for.

Validation Check #

def check_against_control(computed: np.ndarray, published: np.ndarray,
                          tol_m: float = 0.020) -> None:
    """Compare a transformed position against a published NAD83(2011) coordinate."""
    residual = float(np.linalg.norm(computed - published))
    assert residual <= tol_m, (
        f"residual {residual:.4f} m exceeds the {tol_m:.3f} m tolerance — check "
        f"the coordinate epoch and the rotation sign convention before the parameters"
    )

The message names the two likeliest causes in order. A residual of a few decimetres growing with the interval is an epoch problem; a residual of roughly twice the rotation effect that does not grow is a sign convention problem.

Common Mistakes #

The rotation sign convention reversed. The position-vector and coordinate-frame conventions differ by the sign of all three rotations. At ITRF-to-NAD83 magnitudes — tens of milliarcseconds — the difference is about a metre at the Earth’s surface, so it is obvious once tested against control and invisible if never tested. Published tables state which convention they use; record it in the audit block alongside the parameters.

The epoch of the input confused with the epoch of the transformation. They are usually the same and occasionally not. The parameters are evaluated at the instant the two frames are being related, which must be the instant the input coordinates are valid at. Feeding coordinates propagated to 2010.00 into parameters evaluated at 2020.00 is a decimetre-scale error with no diagnostic signature other than a failed comparison.

A static seven-parameter set used for a plate-fixed pair. Some registries publish an ITRF-to-NAD83 operation with no rates for convenience. It is correct only at its reference epoch, and its accuracy statement usually says so. Check the operation’s declared accuracy and its epoch validity before using it, in the same way the fallback logic in automating datum fallback chains in pyproj checks accuracy before selecting an operation.

Sequencing the Two Operations #

A complete conversion from an observed ITRF position to a published NAD83(2011) coordinate is two operations, and their order is fixed by what each one means. First propagate within ITRF from the observation epoch to the epoch at which the transformation is to be evaluated; then apply the fourteen-parameter transformation at that epoch. Doing it the other way — transforming first and then propagating in NAD83 — requires a NAD83 velocity, which is a different quantity from the ITRF velocity and is near zero in the stable interior, so the two routes disagree by very nearly the full plate motion over the interval.

Figure — propagate inside ITRF first, then change frames; the reverse needs a different velocity.

Both operations belong in the audit record separately, with their own accuracies: the propagation carries the velocity uncertainty times the interval, and the transformation carries the parameter uncertainty. Merging them into a single reported figure hides which of the two produced a given millimetre, and it is the propagation term that grows without bound as the interval lengthens.

Frequently Asked Questions #

Which ITRF realisation should I transform from? #

The one your positioning service actually delivered — realisations differ from one another by centimetres, and treating ITRF2014 output as ITRF2020 input introduces exactly that. Post-processing services state the realisation and the epoch in their reports; carry both through the pipeline rather than reconstructing them later from the processing date.

Is WGS84 close enough to ITRF for cadastral work? #

The recent WGS84 realisations agree with the corresponding ITRF realisations at the centimetre level, so for many purposes yes — but “WGS84” without a realisation suffix is an ensemble with about two metres of internal ambiguity, and treating that as ITRF is not defensible. If the source metadata says only WGS84, the honest record says the realisation is unknown, and the uncertainty statement has to reflect it.

Do I need the rates if my epoch is close to the reference epoch? #

Within a year or so of the reference epoch the rate contribution is a couple of centimetres, which is already at the tolerance for cadastral work; beyond that it dominates. Including the rates costs nothing at runtime and removes an entire class of error, so there is no operational reason to omit them.

How do I check my implementation without published control? #

Transform a coordinate forward and back through the same parameter set at the same epoch and require closure to a micrometre — that catches sign, transposition and unit errors in the rotation matrix. Then transform a position at two epochs ten years apart and confirm that the difference matches ten times the rate terms. Neither test needs external data, and together they catch most implementation mistakes short of a wrong parameter table.

Related #

• Time-dependent transformations and plate motion — the parent topic
• Helmert 7-parameter transformations in Python
• Position vector vs coordinate frame rotation conventions
• Computing plate motion model velocities in Python
• Core transformation fundamentals and standards