Handling Epoch and Plate Motion in PyProj

A survey positioned by a global correction service and delivered on a national grid comes back eight centimetres out, consistently, in the direction the continent is moving. Nothing is misconfigured in the usual sense: the horizontal datums were transformed, the geoid was applied, the control fits. The missing piece is that both coordinate systems are correct at different moments, and the transformation between them was asked to ignore time.

This page makes the epoch explicit, applies the transformation with it, and checks that it was actually used.

Why a coordinate needs a date

The Earth’s crust moves. Australia travels roughly 7 cm a year, North America 1–2 cm, western Europe about 2.5 cm. Two families of reference frame handle that differently.

A global frame — the ITRF realisations, and the frames that global correction services broadcast in — is fixed to the Earth as a whole. A point on a moving plate therefore has coordinates that change every year, and a position in such a frame is meaningless without the date it refers to.

A plate-fixed frame — most national grids — is attached to its own plate and defined at a fixed reference epoch. Coordinates in it are stable over decades, which is exactly what a cadastre needs.

Transforming between them requires knowing both dates: where the point was when observed, and where the frame says it should be at its reference epoch. Supply only one and the transformation applies a static approximation, which is wrong by the accumulated motion.

A fixed point in two frames over twenty years A single physical monument shown over twenty years. In the plate-fixed national frame its coordinates never change, because the frame moves with the plate. In the global frame its coordinates drift steadily, accumulating about fifty centimetres over the period at a plate velocity of two and a half centimetres per year. A survey observed in the global frame in the present and delivered on the national grid must therefore be moved backwards by the accumulated motion, and the size of that correction is the vertical gap between the two lines at the observation date. 2006 2011 2016 2021 2026 date of observation coordinate value plate-fixed national frame — constant global frame — drifts with the plate ≈ 50 cm accumulated The correction is the gap at the observation date — which is why a transformation without dates cannot compute it.

Figure 1 — The same monument in two frames. Neither coordinate is wrong; they describe the position at different moments, and the difference is a measurable quantity rather than an error.

Minimal reproducible solution

pyproj handles this through 4D coordinates: transform with a time component, and use a TransformerGroup when you need to see which operations are available and how accurate each is.

import numpy as np
from pyproj import CRS, Transformer
from pyproj.transformer import TransformerGroup


def epoch_aware_transform(lon: np.ndarray, lat: np.ndarray, h: np.ndarray,
                          observation_epoch: float,
                          source_crs: str, target_crs: str) -> tuple:
    """Transform positions with their observation epoch supplied.

    observation_epoch is a decimal year — 2026.62 for mid-August 2026. The
    fourth coordinate is what lets PROJ apply a time-dependent operation; omit
    it and PROJ silently uses a static approximation.
    """
    src, tgt = CRS.from_user_input(source_crs), CRS.from_user_input(target_crs)

    group = TransformerGroup(src, tgt, always_xy=True)
    if not group.transformers:
        raise ValueError(f"no operation available from {source_crs} to {target_crs}")
    if group.unavailable_operations:
        names = [op.name for op in group.unavailable_operations[:2]]
        raise RuntimeError(
            f"a more accurate operation exists but its grid is missing: {names} — "
            "set PROJ_NETWORK=ON or install the grid")

    tr = group.transformers[0]
    t = np.full_like(np.asarray(lon, dtype=float), observation_epoch)
    e, n, z, _ = tr.transform(lon, lat, h, t)
    return e, n, z, tr.description, tr.accuracy

Two guards carry the weight. TransformerGroup exposes the unavailable operations as well as the available ones, which is the only way to learn that PROJ chose a coarser path because a grid was missing — the ordinary Transformer.from_crs simply returns the best it can build and says nothing. And passing the epoch as a fourth coordinate is what selects a time-dependent operation at all; a three-argument call cannot, regardless of how the CRS are declared.

The epoch itself should come from the data rather than from datetime.now():

from datetime import datetime, timezone


def decimal_year(when: datetime) -> float:
    """Convert an acquisition timestamp to the decimal year PROJ expects."""
    if when.tzinfo is None:
        raise ValueError("declare the timezone; a naive timestamp is ambiguous")
    utc = when.astimezone(timezone.utc)
    start = datetime(utc.year, 1, 1, tzinfo=timezone.utc)
    end = datetime(utc.year + 1, 1, 1, tzinfo=timezone.utc)
    return utc.year + (utc - start).total_seconds() / (end - start).total_seconds()

Using the flight date rather than the processing date matters when data is reprocessed: a survey flown in 2024 and reprocessed in 2027 must still be transformed at its 2024 epoch, or the reprocessed result will differ from the original delivery by three years of motion.

Edge-case matrix

Situation Without an epoch With one
Global frame → plate-fixed grid Off by the accumulated motion Correct
Plate-fixed → the same grid No difference No difference
Two global frames, same epoch Small, from frame realisation Same
Reprocessing an old survey Uses today’s motion Uses the flight’s epoch
Slow plate (a few mm/yr) Under a centimetre Same, negligible
Fast plate (7 cm/yr) Decimetres over a decade Correct
Post-seismic deformation Not modelled by a velocity Needs a deformation model
Missing grid for the operation Coarser path chosen silently Detected by the group check

The post-seismic row is the honest limit. A plate velocity model describes steady motion; it does not describe the metre-scale displacements a large earthquake produces, nor the years of decaying afterslip that follow. In affected regions the national agency publishes a deformation model, and the transformation must use it rather than a linear velocity.

What TransformerGroup shows that from_crs hides Two views of the same coordinate transformation request. Using the ordinary constructor, one operation is returned and used, with no indication that it was chosen because a better one could not be built. Using a transformer group, the full candidate list is visible: the most accurate operation requires a grid that is not installed and is listed as unavailable, while the operation actually selected is a lower-accuracy fallback. A note observes that the numeric difference between them is around a decimetre, which is invisible in the returned coordinates and decisive for a survey. Transformer.from_crs returns one transformer no indication of what else existed chosen: 7-parameter fallback accuracy ≈ 0.5 m the call succeeds and says nothing TransformerGroup lists available and unavailable unavailable: time-dependent, grid missing accuracy ≈ 0.02 m available: 7-parameter fallback accuracy ≈ 0.5 m the gap is visible before it is inherited A twenty-five-fold accuracy difference, invisible in the returned coordinates. Checking the group costs one extra object and is the only way to learn the operation was a fallback.

Figure 2 — Why the group matters. The ordinary constructor is a convenience that quietly hides the one fact a survey pipeline needs: whether the operation it got was the operation it wanted.

Verify the fix worked

The check is that the applied displacement matches the plate velocity times the elapsed time, in the direction the plate moves.

import numpy as np


def assert_epoch_correction_plausible(before_en: np.ndarray, after_en: np.ndarray,
                                      years: float, velocity_mm_yr: float,
                                      tol_ratio: float = 0.3) -> None:
    """The shift must match the expected plate motion in size and uniformity."""
    shift = after_en - before_en
    magnitude = np.linalg.norm(shift, axis=1)
    expected_m = abs(velocity_mm_yr) * abs(years) / 1000.0

    median = float(np.median(magnitude))
    assert abs(median - expected_m) <= expected_m * tol_ratio + 0.005, (
        f"applied shift {median * 100:.1f} cm against an expected "
        f"{expected_m * 100:.1f} cm over {years:.1f} years")

    # Plate motion is uniform across a survey-sized area.
    directions = shift / np.maximum(magnitude[:, None], 1e-9)
    spread = float(np.std(np.degrees(np.arctan2(directions[:, 0], directions[:, 1]))))
    assert spread < 5.0, (
        f"shift directions vary by {spread:.1f}° across the survey — plate "
        "motion is uniform at this scale, so this is not plate motion")

The directional-uniformity test is the one that catches a misdiagnosis. Over a few kilometres the plate velocity field is constant to well under a degree, so a correction whose direction varies across the block is doing something else — usually a projection difference being absorbed into what was assumed to be an epoch problem.

When to escalate

  • The correction is applied and a residual offset remains. The frames differ by more than their epochs; check the datum realisations themselves, which is the vertical and horizontal datum question rather than a temporal one.
  • The region has had a significant earthquake since the frame epoch. A linear velocity model is not adequate. The national agency’s published deformation model is the only correct input, and applying plate motion instead will be wrong by the coseismic displacement.
  • Two surveys of the same site years apart disagree by the plate motion. They were transformed at their own epochs and are both correct — in different frames. For change detection, transform both to a common epoch before differencing, and record which epoch that was.

Whatever epoch is used, record it. The epoch belongs in the run manifest beside the CRS strings and the PROJ version, and it belongs in the deliverable’s metadata, because a coordinate without its epoch cannot be transformed correctly by whoever receives it. A client who is told only “EPSG:7844” has been given a frame and not a position in it, and the first thing they will do is transform it at whatever epoch their software defaults to.

← Coordinate Transformation Workflows in PyProj

Change detection between two epochs, done wrong and right Two surveys of the same quarry, flown four years apart and each correctly transformed at its own epoch in a global frame. Differencing them directly attributes the accumulated plate motion — about ten centimetres — to ground movement, producing a uniform apparent displacement across the whole site including areas that did not change. Transforming both to a common epoch first removes that component, leaving only the real change concentrated in the excavated area. A note states that the give-away is the uniformity: real ground movement is never the same everywhere. differenced at different epochs excavated area 10 cm everywhere, including unchanged ground both moved to a common epoch excavated area change only where the ground actually moved Uniformity is the give-away: real ground movement is never identical across an entire site.

Figure 3 — The consequence for repeat surveys. Two correct datasets differenced without a common epoch report the continent’s motion as the client’s.