Computing Plate Motion Model Velocities in Python

A plate motion model turns a position into a velocity with three numbers and a cross product, and for a site in a stable plate interior it is accurate to a millimetre or two per year — enough for almost any cadastral epoch propagation. This guide, part of time-dependent transformations and plate motion, implements the Euler-pole computation, converts the resulting Cartesian velocity into the north-east-up components a surveyor reads, and states precisely where the rigid-plate assumption stops being defensible.

The Euler Pole and the Cross Product #

A rigid plate on a sphere moves by rotating about an axis through the Earth’s centre. That rotation is described by an angular velocity vector $\boldsymbol{\omega}$ whose direction is the rotation axis — the Euler pole — and whose magnitude is the rotation rate. The velocity of any site on the plate is then

Figure — the cross product: velocity is perpendicular to both the axis and the position.

with $\mathbf{X}$ the site’s Cartesian position in metres and $\mathbf{v}$ its velocity in metres per year. Published models give $\boldsymbol{\omega}$ either as Cartesian components in radians per year (often quoted as milliarcseconds per year) or as a pole latitude, pole longitude and rotation rate in degrees per million years. Both parameterisations describe the same vector, and converting between them is where most of the unit errors live.

Two properties of the cross product explain the behaviour of the result. The velocity is always perpendicular to both the rotation axis and the position vector, so it is horizontal only to the extent the position is perpendicular to the axis — a rigid rotation therefore produces a small but non-zero vertical component, which is real and should not be discarded. And the magnitude grows with the perpendicular distance from the rotation axis, which is why sites far from the Euler pole move fastest.

Complete Runnable Implementation #

from __future__ import annotations

import numpy as np

A = 6378137.0
F = 1.0 / 298.257222101
E2 = F * (2.0 - F)

# Milliarcseconds per year -> radians per year: the unit most tables publish.
MAS_PER_YEAR_TO_RAD = np.pi / (180.0 * 3600.0 * 1000.0)


def omega_from_pole(pole_lat_deg: float, pole_lon_deg: float,
                    rate_deg_per_myr: float) -> np.ndarray:
    """Cartesian angular velocity (rad/yr) from a published Euler pole.

    `rate_deg_per_myr` is degrees per million years, the convention used by the
    geological plate-motion literature.
    """
    rate = np.radians(rate_deg_per_myr) / 1.0e6            # rad/yr
    plat, plon = np.radians(pole_lat_deg), np.radians(pole_lon_deg)
    return rate * np.array([
        np.cos(plat) * np.cos(plon),
        np.cos(plat) * np.sin(plon),
        np.sin(plat),
    ], dtype=np.float64)


def omega_from_mas(wx: float, wy: float, wz: float) -> np.ndarray:
    """Cartesian angular velocity (rad/yr) from mas/yr components."""
    return np.array([wx, wy, wz], dtype=np.float64) * MAS_PER_YEAR_TO_RAD


def geodetic_to_cartesian(lat_deg: float, lon_deg: float, h_m: float) -> np.ndarray:
    """GRS80 geodetic -> Earth-centred Cartesian, metres."""
    lat, lon = np.radians(lat_deg), np.radians(lon_deg)
    s = np.sin(lat)
    n = A / np.sqrt(1.0 - E2 * s * s)
    return np.array([
        (n + h_m) * np.cos(lat) * np.cos(lon),
        (n + h_m) * np.cos(lat) * np.sin(lon),
        (n * (1.0 - E2) + h_m) * s,
    ], dtype=np.float64)


def enu_rotation(lat_deg: float, lon_deg: float) -> np.ndarray:
    """Rotation from Earth-centred Cartesian into local east-north-up."""
    lat, lon = np.radians(lat_deg), np.radians(lon_deg)
    sl, cl = np.sin(lat), np.cos(lat)
    so, co = np.sin(lon), np.cos(lon)
    return np.array([
        [-so, co, 0.0],
        [-sl * co, -sl * so, cl],
        [cl * co, cl * so, sl],
    ], dtype=np.float64)


def plate_velocity(lat_deg: float, lon_deg: float, h_m: float,
                   omega: np.ndarray) -> tuple[np.ndarray, np.ndarray]:
    """Site velocity from a rigid plate rotation.

    Returns (cartesian_velocity, enu_velocity), both in metres per year. The ENU
    triple is what belongs in a survey record; the Cartesian one is what an epoch
    propagation in a geocentric frame consumes.
    """
    xyz = geodetic_to_cartesian(lat_deg, lon_deg, h_m)
    v_xyz = np.cross(omega, xyz)                    # rad/yr x m = m/yr
    v_enu = enu_rotation(lat_deg, lon_deg) @ v_xyz
    return v_xyz, v_enu

Parameter Reference #

Worked Example #

Using an approximate ITRF2014 rotation for the North American plate — Euler pole at −4.85° N, −85.98° E rotating at 0.194 °/Myr — the velocity at a monument near Portland, Oregon at 45.500° N, 122.600° W, 61 m:

Figure — horizontal ITRF speeds for four plates, at a representative mid-plate site.

import numpy as np

omega = omega_from_pole(-4.85, -85.98, 0.194)
v_xyz, v_enu = plate_velocity(45.500, -122.600, 61.0, omega)
print(np.round(v_enu * 1000, 2))    # east, north, up in mm/yr
# [-13.02   8.05  -0.06]

Roughly 13 mm/yr west and 8 mm/yr north — the familiar bulk motion of the North American plate in a global frame, and the reason ITRF coordinates in the conterminous United States are unusable without an epoch. The vertical component is a fraction of a millimetre, which is the rigid-rotation answer and is not a claim about actual uplift.

Validation Check #

def check_rigid_rotation(omega: np.ndarray, lat: float, lon: float, h: float) -> None:
    """A rigid rotation has two exact properties worth asserting."""
    xyz = geodetic_to_cartesian(lat, lon, h)
    v = np.cross(omega, xyz)
    # 1. The velocity is perpendicular to the rotation axis.
    assert abs(float(v @ omega)) < 1e-12, "velocity is not perpendicular to omega"
    # 2. A site ON the Euler pole axis does not move.
    on_axis = omega / np.linalg.norm(omega) * np.linalg.norm(xyz)
    assert np.linalg.norm(np.cross(omega, on_axis)) < 1e-12, "pole site is moving"

Both assertions are exact consequences of the cross product rather than tolerances, so a failure means a genuine implementation error — usually a transposed rotation matrix or a pole converted with the wrong unit.

Common Mistakes #

Degrees per million years used directly as radians per year. The two differ by about 5.7 × 10⁷, so the resulting velocity is absurd rather than subtly wrong — which is fortunate. The subtler version is mixing the two conventions across a parameter set, where one component was converted and the others were not.

Figure — where a rigid rotation is enough, and where only a grid will do.

East-north-up order assumed to be north-east-up. Velocity grids commonly publish north first, while the rotation matrix above returns east first. Swapping them puts the motion at ninety degrees to the truth, and the magnitude looks perfectly plausible. Name the components in the return value and in the audit record rather than relying on position.

A rigid model used across a plate boundary. The whole premise is that the plate is rigid. Within a few hundred kilometres of an active boundary it is not, and the residual velocity after removing the rigid rotation can reach several millimetres a year — with exactly the spatial structure that matters for adjacent parcels. Use a published velocity grid there, as described in propagating coordinates between epochs with velocity grids.

Comparing a Rigid Model Against a Published Grid #

The decision between a plate rotation and a velocity grid is easy to make empirically. Evaluate both at a scatter of positions across the working area and difference them: agreement within a millimetre a year everywhere means the rigid model is adequate there and the grid buys nothing but a file dependency. A systematic disagreement growing towards one edge of the area means you are approaching a boundary or a deforming zone, and the grid is describing something the rotation cannot. Run this comparison once per working area rather than per job, and record the outcome — it is the justification for whichever model the pipeline then uses.

Frequently Asked Questions #

Which plate motion model should I use? #

Whichever the frame you are working in specifies. Models are defined relative to a frame, and mixing a model derived in one realisation with coordinates in another introduces a systematic velocity error of a millimetre or two per year. Where a national agency publishes both a frame and a companion model, use the pair — the consistency matters more than any accuracy difference between competing models.

Does the ellipsoidal height affect the velocity? #

Only slightly, and predictably: the cross product scales with the position vector, so a height of a few kilometres changes the velocity by a few parts in a thousand — well under a tenth of a millimetre a year. Getting the height roughly right is sufficient; getting the latitude and longitude right is not optional, because the direction of the velocity depends on them.

Can a plate rotation replace a velocity grid entirely? #

In a stable interior, generally yes, and it has the advantage of needing no data files or extent checks. What it cannot do is describe anything that is not rigid-plate motion: strain near a boundary, post-glacial rebound, subsidence over a pumped aquifer, or the offset left by an earthquake. A useful diagnostic is to compare the rigid velocity against a published grid at a handful of points across your working area — if they agree to a millimetre a year, the rigid model is adequate there.

Why does the model give a small vertical velocity? #

Because a rotation about an axis through the Earth’s centre is not exactly parallel to the local horizontal plane at every point on an ellipsoid. The resulting vertical component is a fraction of a millimetre a year and is a geometric artefact of the model, not a prediction of uplift. Real vertical motion — rebound, subsidence — comes from a velocity grid or a local time series, not from a plate rotation.

Related #

• Time-dependent transformations and plate motion — the parent topic
• Propagating coordinates between epochs with velocity grids
• ITRF to NAD83(2011) transformation in Python
• Geodetic conversion math: ellipsoid to Cartesian
• Helmert 7-parameter transformations in Python