CoordTransform Cadastral & High-Precision Python

Time-Dependent Transformations and Plate Motion in Python

A coordinate in a modern geocentric frame is a position and an instant, and this topic — part of core transformation fundamentals and standards — is about the instant. The ground moves: a monument in the conterminous United States travels roughly 15–25 mm per year in ITRF, one in south-eastern Australia closer to 70 mm per year, and one on an active margin can move differently from its neighbour ten kilometres away. That motion is invisible while everything stays inside one plate-fixed realization, and it becomes the dominant error the moment a coordinate crosses into a global frame, is compared against an observation made years later, or is submitted without the epoch that makes it interpretable. This page covers what a coordinate epoch is, how to propagate a position between epochs deterministically, and what has to appear in the audit record for the result to mean anything.

The distinction that organises everything below is between a datum transformation and an epoch propagation. A datum transformation moves a coordinate between reference frames — NAD83(2011) to ITRF2014, say — and is the subject of Helmert 7-parameter transformations in Python. An epoch propagation moves a coordinate within one frame from one instant to another, using a velocity. They are different operations with different uncertainties, and folding them into one step is the most common way a time-dependent workflow becomes unreproducible.

Reference Frames, Epochs, and Why Both Are Required

A datum realization such as NAD83(2011) fixes the frame; the coordinate epoch fixes when the coordinates in that frame were valid. NAD83 is plate-fixed: it rotates with the North American plate, so a stable monument in Kansas has coordinates that barely change with time and an epoch that is nearly a formality. ITRF is global: the same monument moves about 15 mm per year in it, so its ITRF coordinates are meaningless without the epoch, and a decade of neglect is 150 mm — three to ten times a typical cadastral tolerance.

A datum transformation and an epoch propagation are different operations A coordinate is identified by a frame and an epoch. A datum transformation changes the frame and leaves the epoch alone, using published parameters and their rates. An epoch propagation changes the epoch and leaves the frame alone, using a site velocity. A complete conversion applies both, in sequence, and records them as two operations with two independent uncertainties rather than one merged figure. Coordinate frame + epoch Datum transformation changes the FRAME Epoch propagation changes the EPOCH Two records, two uncertainties never merged into one figure 14 parameters velocity

Figure — two different operations that are constantly confused: one moves frames, one moves time.

Two consequences follow. First, a transformation between a plate-fixed and a global frame is defined at an epoch, and a published parameter set for such a pair normally carries both a reference epoch and rates of change for the seven parameters. Applying the parameters without the rates, at an epoch far from the reference epoch, produces an error that grows linearly with the interval. Second, a residual computed between a field observation and a published control coordinate is only meaningful once both are at the same epoch. In an active deformation zone, a residual computed across ten years can exceed the tolerance being tested while the survey itself is faultless — and the usual response, widening the tolerance, hides a problem that a one-line epoch propagation would have removed.

The Propagation Model

Propagating a position from epoch t0t_0 to epoch t1t_1 is, to the accuracy anyone needs over a few decades, a linear extrapolation along a velocity vector:

Accumulated motion against elapsed years at four plate rates Four traces of accumulated horizontal motion in metres against elapsed years from zero to twelve. A stable interior at 2 millimetres a year reaches 24 millimetres after twelve years. North America at 18 millimetres a year crosses a 20 millimetre cadastral tolerance in about 1.1 years. Europe at 25 millimetres a year crosses it in under a year, and eastern Australia at 68 millimetres a year crosses it in about four months. 0.00 0.25 0.50 0.75 0 2 4 6 8 10 12 elapsed years m stable interior 2 mm/yr North America 18 mm/yr Europe 25 mm/yr Australia 68 mm/yr

Figure — accumulated ground motion against elapsed time, at four representative plate rates.

X(t1)=X(t0)+v (t1−t0)\mathbf{X}(t_1) = \mathbf{X}(t_0) + \mathbf{v}\,(t_1 - t_0)

where X\mathbf{X} is the Cartesian position in metres and v\mathbf{v} the site velocity in metres per year, expressed in the same frame. Everything interesting is in where v\mathbf{v} comes from. There are three sources, in ascending order of fidelity:

  • A rigid plate rotation. The plate is treated as rotating about an Euler pole with angular velocity ω\boldsymbol{\omega}, and the site velocity is the cross product v=ω×X\mathbf{v} = \boldsymbol{\omega} \times \mathbf{X}. This is exact for a rigid plate and captures the bulk motion of a stable interior to within a millimetre or two per year.
  • A velocity grid. A gridded model of horizontal and vertical velocity, interpolated at the site. This captures what a rigid rotation cannot: the strain accumulation near a plate boundary, post-glacial rebound, subsidence over an aquifer.
  • A local determination. The site’s own velocity from a continuous GNSS time series, which is the only option where the deformation is not modelled at all.

The rigid-rotation route and the grid route are covered in detail in computing plate motion model velocities in Python and propagating coordinates between epochs with velocity grids respectively.

Production Implementation: an Epoch-Aware Position

The type below refuses to exist without an epoch when the frame requires one, and refuses to compare two positions at different epochs. That is deliberate: the failure mode this topic exists to prevent is an epoch that was never recorded, and a type that permits it will eventually be handed one.

from __future__ import annotations

from dataclasses import dataclass
from decimal import Decimal, ROUND_HALF_EVEN

import numpy as np

# Frames whose coordinates are only interpretable with a stated epoch. A plate-fixed
# frame still HAS an epoch; a global frame cannot be read without one.
EPOCH_REQUIRED = {"ITRF2014", "ITRF2020", "IGS20", "WGS84(G2139)"}


@dataclass(frozen=True)
class EpochPosition:
    """A Cartesian position with the frame and instant that make it meaningful."""

    xyz: np.ndarray          # (3,) float64, metres
    frame: str               # e.g. "ITRF2014" or "NAD83(2011)"
    epoch: float | None      # decimal year, e.g. 2010.00

    def __post_init__(self) -> None:
        if self.xyz.shape != (3,) or self.xyz.dtype != np.float64:
            raise ValueError("xyz must be a (3,) float64 array of metres")
        if self.frame in EPOCH_REQUIRED and self.epoch is None:
            raise ValueError(f"{self.frame} coordinates require a coordinate epoch")

    def propagate(self, velocity: np.ndarray, to_epoch: float) -> "EpochPosition":
        """Move the position within its own frame to `to_epoch`.

        `velocity` is metres per year in the SAME frame. This is an epoch
        propagation, not a datum transformation: the frame is unchanged, and the
        operation is recorded separately in the audit trail (ISO 19111 treats them
        as distinct coordinate operations).
        """
        if self.epoch is None:
            raise ValueError("cannot propagate a position with no source epoch")
        dt = to_epoch - self.epoch
        return EpochPosition(self.xyz + velocity * dt, self.frame, to_epoch)

    def displacement_from(self, other: "EpochPosition") -> float:
        """Metres between two positions — refuses to compare across epochs."""
        if self.frame != other.frame:
            raise ValueError(f"frame mismatch: {self.frame} vs {other.frame}")
        if self.epoch != other.epoch:
            raise ValueError(
                f"epoch mismatch: {self.epoch} vs {other.epoch}; propagate one "
                f"to the other's epoch before comparing"
            )
        return float(np.linalg.norm(self.xyz - other.xyz))


def quantise(value: float, places: int = 4) -> Decimal:
    """Deterministic output rounding — round-half-to-even, as IEEE 754 specifies."""
    quantum = Decimal(1).scaleb(-places)
    return Decimal(repr(value)).quantize(quantum, rounding=ROUND_HALF_EVEN)

The displacement_from guard is the part worth copying even if nothing else here is useful. Comparing coordinates at different epochs is not an error the arithmetic can detect — the numbers subtract perfectly well — so it has to be a rule enforced by the type.

Motion Rates and What They Cost

Rates vary by an order of magnitude between plates, and the interval that matters follows directly. The table below gives representative horizontal ITRF velocities and the interval over which the accumulated motion reaches 20 mm, a typical cadastral horizontal tolerance.

Region Typical ITRF velocity 20 mm accumulates in Practical consequence
Central North America 15–20 mm/yr ~1.1 years Epoch required for any ITRF work
Western Europe 22–27 mm/yr ~0.8 years Epoch required; ETRS89 is the plate-fixed answer
Eastern Australia 65–70 mm/yr ~0.3 years Epoch dominates every other error term
Stable interior, plate-fixed frame < 2 mm/yr > 10 years Epoch still recorded, rarely applied
Active margin (e.g. coastal California) 20–40 mm/yr, non-uniform ~0.6 years Rigid rotation insufficient; grid required

The last row is the one that changes the implementation rather than just the numbers. Where the velocity field is non-uniform, a rigid plate rotation removes the bulk motion and leaves a residual velocity of several millimetres per year that only a grid captures — the same “smooth model versus gridded model” decision that NADCON vs NTv2 makes for horizontal datum shifts.

Compliance: What the Audit Record Must Carry

An epoch propagation is a coordinate operation, and ISO 19111 expects it to be described like one. The record needs the source epoch, the target epoch, the velocity source (model name and version, or grid file and checksum, or the identifier of the local determination), the velocity actually applied at that position, and the uncertainty of that velocity. That last item is routinely omitted and is often the largest term in the propagated uncertainty: a velocity known to ±1 mm/yr, propagated over fifteen years, contributes ±15 mm, which by itself can consume a cadastral budget.

Fields an epoch propagation writes to the audit record Six fields. The source epoch and the target epoch bound the interval. The velocity source names the model and version, or the grid file and its checksum. The velocity applied is the value actually interpolated at that position, not the model average. The velocity uncertainty is required and routinely omitted, and over a long interval it can be the largest term in the budget. The operation identifier records the propagation as its own coordinate operation rather than folding it into the datum shift. source_epoch decimal year where the coordinates were valid target_epoch decimal year where they are valid now velocity_source model + version or grid file + checksum velocity_applied m/yr, N-E-U at this position, not the average velocity_sigma m/yr x interval — often the largest term operation_id string its own operation, not part of the shift

Figure — the six fields an epoch propagation contributes to the audit record.

Two rules keep the record honest. The propagation is recorded separately from any datum transformation applied in the same run, because their uncertainties are independent and a reviewer needs to see both. And an absent epoch is a rejection rather than a default: defaulting to the current date produces a result that changes depending on when the pipeline ran, and leaves no trace of the assumption. The general shape of that record is set out in compliance report generation for agency submission.

Validation: the Round-Trip and the Time-Reversal Test

Two cheap tests catch most implementation errors. The round trip propagates forward and back and requires the original position to within a micrometre; it catches sign errors in the interval and unit errors in the velocity. The time-reversal test propagates a position forward ten years and checks that the displacement equals ten times the annual velocity; it catches an interval computed in days, or in milliseconds since the epoch, which is a mistake that survives the round trip because it is symmetric.

import numpy as np

def check_propagation(pos: EpochPosition, velocity: np.ndarray) -> None:
    """Round-trip and time-reversal checks for an epoch propagation."""
    there = pos.propagate(velocity, pos.epoch + 10.0)
    back = there.propagate(velocity, pos.epoch)
    assert np.allclose(back.xyz, pos.xyz, atol=1e-6), "round trip is not closing"

    moved = np.linalg.norm(there.xyz - pos.xyz)
    expected = 10.0 * np.linalg.norm(velocity)
    assert abs(moved - expected) < 1e-6, (
        f"displacement {moved:.6f} m over 10 years does not match "
        f"10 x the annual velocity ({expected:.6f} m) — check the interval units"
    )

Failure Modes and Their Mitigations

  • Epoch defaulted to “now”. The output changes with the wall-clock time of the run and the assumption is invisible. Reject the input instead, and if a default is genuinely unavoidable, record it as an explicit assumption.
  • Velocity applied in the wrong frame. A velocity expressed in a plate-fixed frame is near zero; the same site’s ITRF velocity is centimetres per year. Applying one where the other belongs produces an error equal to the plate motion itself. Carry the frame with the velocity, as the type above does.
  • Rates of the transformation parameters ignored. A frame-to-frame parameter set with rates, applied at its reference epoch only, drifts by the rate times the interval. This is a silent few millimetres per year.
  • Rigid rotation used on an active margin. The bulk motion is removed and a non-uniform residual of several millimetres per year is left behind, concentrated exactly where the boundary disputes are. Use a velocity grid where one is published.
  • Epoch lost at the format boundary. Most exchange formats have no epoch field. It belongs in the accompanying metadata block, alongside the operation identifier and grid checksum — see exporting ISO 19111 metadata for cadastral deliverables.

Choosing a Reference Epoch for a Project

Every project that touches a time-dependent frame has to pick one epoch and hold everything to it, and the choice is a project decision rather than a technical one. Three conventions are in common use, and each is defensible in its own context.

The national reference epoch is whatever the agency’s published control is expressed at — 2010.00 for NAD83(2011), 2020.00 for GDA2020, and so on. Adopting it means the published control needs no propagation, which removes an operation and its uncertainty from every comparison. This is the right default for cadastral work, because the deliverable is going to be compared against that control eventually whatever you do.

The epoch of observation means coordinates are published at the instant they were measured. It is the most faithful record of what was observed and the least convenient for anyone downstream, because every comparison then requires a propagation. It suits a monitoring network, where the whole point is how positions change, and suits a boundary survey badly.

A project epoch — often the midpoint of a long campaign — is the compromise: one propagation per observation, applied once at the start, and then a self-consistent dataset. It is worth choosing deliberately rather than inheriting whatever the first day of field work happened to be.

Whichever is chosen, two rules keep it usable. The epoch must appear in the deliverable, not only in the project file, because the two are separated the moment data is shared. And a single project must not mix epochs silently: if a late addition arrives at a different epoch, propagate it and record the propagation, rather than appending coordinates that look like all the others and are not.

Datum Transformations That Carry Rates

The seven-parameter transformations described elsewhere in this section are static: one set of numbers, valid at all times. Transformations between a plate-fixed frame and a global one are not, and the published parameter sets for those pairs come with a reference epoch and seven rates of change, one per parameter. Evaluating the parameters is a linear step before the transformation itself:

P(t)=P(tr)+P˙ (t−tr)P(t) = P(t_r) + \dot{P}\,(t - t_r)

The rotation rates dominate, because they express the plate rotation and are multiplied by the Earth’s radius when they reach a coordinate — roughly 31 mm per milliarcsecond at the surface. A parameter set evaluated at its reference epoch and reused for a decade is therefore out by 150 to 200 mm at mid-latitudes, which is an order of magnitude beyond cadastral tolerance and is produced by omitting a single argument. The full implementation is in ITRF to NAD83(2011) transformation in Python, and the sign-convention trap it shares with every Helmert transformation is covered in position vector vs coordinate frame rotation conventions.

Frequently Asked Questions

Does a plate-fixed frame like NAD83 need an epoch at all?

Yes, though it is applied less often. NAD83(2011) coordinates in the stable interior change by well under a millimetre a year, so the epoch rarely alters the numbers — but it is still part of the coordinate’s identity, and it becomes essential the moment the position is transformed into a global frame or compared against one. On an active margin inside the same frame, it matters immediately: the plate-fixed model does not describe local deformation, and positions there move measurably in NAD83 too.

Is the epoch the same thing as the date of survey?

Usually but not necessarily, and the difference is worth stating explicitly. The coordinate epoch is the instant the coordinates are valid for; the survey date is when the observations were made. They coincide when coordinates are published at the epoch of observation, and they differ whenever coordinates have been propagated to a common reference epoch for a project — which is normal practice. Record both.

Can I ignore vertical velocity?

Only if you have checked. Vertical rates are smaller than horizontal ones in most stable regions — a millimetre or two a year — but post-glacial rebound reaches a centimetre a year in parts of Scandinavia and northern Canada, and subsidence over a pumped aquifer can exceed several centimetres a year. Against a vertical tolerance that is typically tighter than the horizontal one, both are significant over a decade.

How do I choose between a plate motion model and a velocity grid?

By whether the deformation at your site is described by rigid-plate motion. In a stable interior, a plate rotation is accurate to a millimetre or two per year and needs no data files. Near a boundary, in a subsidence basin, or anywhere a published velocity grid exists for the region, use the grid — its existence is usually evidence that the rigid model was not good enough there.