Propagating Coordinates Between Epochs with Velocity Grids

Moving a survey-grade position from the epoch it was observed at to the epoch a deliverable requires is a single linear step — position plus velocity times interval — and the whole difficulty is getting a trustworthy velocity at the point. This guide, part of time-dependent transformations and plate motion within core transformation fundamentals and standards, implements the gridded route: reading a horizontal and vertical velocity model, interpolating it at the site, applying it over the interval, and recording enough that the result can be reproduced. The tolerance it has to satisfy is the same as any other cadastral operation — the propagated uncertainty must stay inside the horizontal and vertical budgets set out in tuning transformation thresholds for survey-grade work.

Why a Grid Rather Than a Plate Rotation #

A rigid plate rotation describes the bulk motion of a plate interior and nothing else. Near a plate boundary the crust deforms: velocities vary smoothly but non-uniformly across tens of kilometres, so two monuments in the same county can move differently by several millimetres a year. A velocity grid stores that field on a regular lattice — typically a national model published on a grid of a few arc-minutes, carrying north, east and up rates in millimetres per year — and interpolating it recovers the local velocity rather than the plate average.

Figure — the propagation path, with the two places it is allowed to refuse.

The arithmetic is unchanged from the rigid case. For a site at $(\varphi, \lambda)$ with interpolated velocity components $v_n, v_e, v_u$ in metres per year, the propagated position at epoch $t_1$ from epoch $t_0$ is

where $M$ and $N$ are the meridian and prime-vertical radii of curvature in metres, and the latitude and longitude increments come out in radians. Working in geodetic coordinates rather than Cartesian ones keeps the north and east components separate, which is what the grid publishes and what a surveyor wants to see in the audit record.

Complete Runnable Implementation #

from __future__ import annotations

from dataclasses import dataclass

import numpy as np

# GRS80 — the ellipsoid behind NAD83, ETRS89, GDA2020 and ITRF alike.
A = 6378137.0
F = 1.0 / 298.257222101
E2 = F * (2.0 - F)


@dataclass(frozen=True)
class VelocityGrid:
    """A regular lattice of site velocities, in metres per year.

    `vn`, `ve`, `vu` are (nlat, nlon) arrays indexed south-to-north and
    west-to-east, matching the row order of the published model. `lat0`/`lon0`
    are the south-west node in degrees and `dlat`/`dlon` the node spacing.
    """

    lat0: float
    lon0: float
    dlat: float
    dlon: float
    vn: np.ndarray
    ve: np.ndarray
    vu: np.ndarray

    def __post_init__(self) -> None:
        if not (self.vn.shape == self.ve.shape == self.vu.shape):
            raise ValueError("velocity components must share a shape")
        if self.vn.dtype != np.float64:
            raise ValueError("velocity grids must be float64: float32 costs ~0.1 mm/yr")

    def interpolate(self, lat: float, lon: float) -> np.ndarray:
        """Bilinear velocity at (lat, lon) in degrees -> (vn, ve, vu) in m/yr.

        Raises outside the grid extent. Extrapolating a deformation model past its
        published extent produces a number with no defensible uncertainty, which is
        the same rule the grid-shift readers apply.
        """
        nlat, nlon = self.vn.shape
        fi = (lat - self.lat0) / self.dlat
        fj = (lon - self.lon0) / self.dlon
        if not (0.0 <= fi <= nlat - 1 and 0.0 <= fj <= nlon - 1):
            raise ValueError(
                f"({lat:.6f}, {lon:.6f}) is outside the velocity grid extent"
            )
        i, j = int(np.floor(fi)), int(np.floor(fj))
        i = min(i, nlat - 2)
        j = min(j, nlon - 2)
        u, v = fi - i, fj - j
        w = np.array([(1 - u) * (1 - v), (1 - u) * v, u * (1 - v), u * v])
        out = np.empty(3, dtype=np.float64)
        for k, comp in enumerate((self.vn, self.ve, self.vu)):
            corners = np.array([comp[i, j], comp[i, j + 1],
                                comp[i + 1, j], comp[i + 1, j + 1]], dtype=np.float64)
            if not np.all(np.isfinite(corners)):
                raise ValueError("velocity grid carries a no-data node at this position")
            out[k] = float(w @ corners)
        return out


def radii_of_curvature(lat_rad: float) -> tuple[float, float]:
    """Meridian (M) and prime-vertical (N) radii in metres, GRS80."""
    s = np.sin(lat_rad)
    w = np.sqrt(1.0 - E2 * s * s)
    n = A / w
    m = A * (1.0 - E2) / (w ** 3)
    return float(m), float(n)


def propagate_epoch(
    lat_deg: float,
    lon_deg: float,
    h_m: float,
    grid: VelocityGrid,
    from_epoch: float,
    to_epoch: float,
) -> tuple[float, float, float, np.ndarray]:
    """Propagate a geodetic position within its frame between two epochs.

    Returns (lat_deg, lon_deg, h_m, velocity_used). The frame is unchanged: this
    is an epoch propagation, recorded as its own coordinate operation.
    """
    v = grid.interpolate(lat_deg, lon_deg)                    # m/yr, north-east-up
    dt = to_epoch - from_epoch                                # decimal years
    lat_rad = np.radians(lat_deg)
    m, n = radii_of_curvature(lat_rad)

    dlat = (v[0] * dt) / m                                    # radians
    dlon = (v[1] * dt) / (n * np.cos(lat_rad))                # radians
    return (
        float(lat_deg + np.degrees(dlat)),
        float(lon_deg + np.degrees(dlon)),
        float(h_m + v[2] * dt),
        v,
    )

Parameter Reference #

Worked Example #

A monument at 45.500000° N, 122.600000° W with an ellipsoidal height of 61.230 m, observed at epoch 2015.20, is required at the deliverable’s reference epoch of 2020.00. The interpolated velocity is 7.4 mm/yr north, −12.1 mm/yr east and −1.6 mm/yr up.

Figure — velocity uncertainty at ±1 mm/yr, propagated over four intervals.

import numpy as np

grid = VelocityGrid(
    lat0=45.0, lon0=-123.0, dlat=0.5, dlon=0.5,
    vn=np.full((3, 3), 0.0074), ve=np.full((3, 3), -0.0121), vu=np.full((3, 3), -0.0016),
)
lat, lon, h, v = propagate_epoch(45.5, -122.6, 61.230, grid, 2015.20, 2020.00)
print(f"{lat:.8f} {lon:.8f} {h:.4f}")
print(f"applied {v * 1000} mm/yr over {2020.00 - 2015.20:.2f} yr")
# 45.50000032 -122.60000074 61.2223
# applied [ 7.4 -12.1  -1.6] mm/yr over 4.80 yr

The horizontal position has moved 35 mm north and 58 mm east, and the height has dropped 8 mm. Against a 20 mm horizontal tolerance, ignoring the propagation would have failed the deliverable on its own — with no arithmetic error anywhere in the pipeline.

Validation Check #

def assert_within_budget(v: np.ndarray, dt: float, sigma_v: float, tol_h: float) -> None:
    """A propagation is only admissible while its own uncertainty fits the budget.

    `sigma_v` is the 1-sigma velocity uncertainty in m/yr (from the model's own
    accuracy statement); over a long interval it can exceed the tolerance by itself.
    """
    propagated_sigma = sigma_v * abs(dt)
    assert propagated_sigma <= tol_h, (
        f"velocity uncertainty {propagated_sigma * 1000:.1f} mm over {dt:.2f} yr "
        f"exceeds the {tol_h * 1000:.1f} mm horizontal budget"
    )

assert_within_budget(v, 4.80, sigma_v=0.0010, tol_h=0.020)   # 4.8 mm — passes

Common Mistakes #

Millimetres per year read as metres per year. Published velocity models almost always state rates in mm/yr, and a grid loaded without the conversion propagates a thousand times too far — a monument that should move 35 mm moves 35 m. The tell is obvious once you look, which is exactly why the loader should assert that every rate is within ±0.1 m/yr, as the dataclass above does.

Figure — three unit and convention errors, and what each does to a 35 mm propagation.

The interval computed in days. Subtracting two dates and getting days, then multiplying by an annual rate, overstates the motion by a factor of about 365. It survives a round trip, because the error is symmetric, which is why the time-reversal check in the parent topic tests the displacement against the annual rate directly.

Vertical rate applied to an orthometric height. The velocity model’s up component describes the motion of the ground relative to the ellipsoid. Applying it to an orthometric height silently assumes the geoid moves with the ground, which is not what the model says. Propagate the ellipsoidal height, then re-derive the orthometric height using the approach in applying geoid undulation for orthometric heights.

Frequently Asked Questions #

Is linear propagation good enough over twenty years? #

For steady tectonic motion, yes: the departure from linearity in a stable field is well under a millimetre a decade, far below the velocity model’s own uncertainty. It is not good enough across an earthquake, which displaces a site instantly by an amount no velocity describes, and it is questionable across a period of changing groundwater extraction. Both cases need a co-seismic or time-series correction rather than a longer extrapolation.

What uncertainty should I assign to the velocity? #

The one the model publishes, interpolated the same way the velocity is — most national models supply a per-node uncertainty grid, and using a single national figure discards real spatial variation. Where nothing is published, one to two millimetres a year is a defensible working figure for a stable interior and clearly optimistic near a boundary. Whatever the number, propagate it: over fifteen years, 1 mm/yr becomes 15 mm.

Should I propagate every point, or the dataset as a whole? #

Every point, using its own interpolated velocity, whenever the dataset spans more than a few tens of kilometres in a deforming region. Applying one velocity to a whole county reintroduces exactly the error the grid exists to remove. For a small site in a stable interior a single velocity is fine, and stating that simplification in the audit record makes it reviewable.

Can I propagate and transform the datum in one step? #

Numerically, yes — and it is often how a published pipeline is written. What must not be merged is the record: the epoch propagation and the datum transformation have independent uncertainties, and a reviewer needs to see both, along with the velocity source and the operation identifier. Combining them into one line in the report hides which of the two produced a given millimetre.

Related #

• Time-dependent transformations and plate motion — the parent topic
• Computing plate motion model velocities in Python
• Recording coordinate epochs in cadastral deliverables
• Bilinear vs bicubic interpolation for grid shifts
• Core transformation fundamentals and standards