Coastal & Marine Spatial Analysis Pipelines

Computing Great-Circle Speed From AIS Fixes

Speed between two AIS fixes is distance over time, and both terms have a trap. The distance is on a curved, flattened surface and the obvious subtraction of coordinates is wrong by a factor that varies with latitude. The time comes from a reconstructed timestamp that can be equal, or negative, between consecutive rows. This page gets both right, within the speed and heading workflow.

Root cause: a degree of longitude shrinks with latitude

Subtracting latitude and longitude and applying Pythagoras treats the two axes as equivalent. They are not: one degree of latitude is about 111 km everywhere, while one degree of longitude is 111 km at the equator, 78 km at 45°, and 0 at the pole.

import numpy as np

def wrong_euclidean_km(lat1, lon1, lat2, lon2):
    """The mistake, written out. Wrong by cos(latitude) on the longitude axis."""
    return np.hypot(lat2 - lat1, lon2 - lon1) * 111.32


print(wrong_euclidean_km(60.0, 5.0, 60.0, 6.0))    # 111.32 km — claimed
# The true east-west distance at 60°N is 111.32 * cos(60°) = 55.66 km.
# The computed speed is therefore exactly 2× the real one at this latitude.

At 60° N — the Norwegian coast, the northern North Sea, the Gulf of Alaska — every east-west speed computed this way is double. The error is smooth, systematic and latitude-dependent, so it survives averaging and produces a fleet whose apparent speed rises with latitude.

Haversine fixes that by treating the earth as a sphere, and introduces a smaller error of its own: the earth is flattened by about 1/298, so a spherical distance differs from the true ellipsoidal one by up to roughly 0.5%. At AIS speeds and typical fix intervals that is a few hundredths of a knot — negligible for a traffic study, and not negligible for a fuel model where the whole effect being measured may be a few per cent.

Three Methods, Three Very Different Errors A chart of distance error as a percentage against latitude from the equator to eighty degrees. The Euclidean-in-degrees curve starts near zero at the equator, where a degree of longitude and a degree of latitude are similar, and rises steeply, passing forty per cent at sixty degrees and exceeding four hundred per cent near eighty. The haversine curve is a nearly flat band around one half of one per cent across the whole range, reflecting the earth's flattening rather than latitude. The geodesic inverse is a flat line at zero, exact by definition. A note records that the Euclidean error is smooth and systematic, so it survives averaging and produces a fleet whose apparent speed rises with latitude. latitude (degrees) distance error Euclidean in degrees haversine — flat, ≈ 0.5% geodesic inverse — exact 0406080 The first error is smooth and latitude-dependent, so it survives averaging and looks like a real geographic pattern.

Step-by-step with production code

1. Vectorised geodesic, with the time guard 

import logging

import numpy as np
import pandas as pd
from pyproj import Geod

logger = logging.getLogger(__name__)
GEOD = Geod(ellps="WGS84")
MS_TO_KN = 1.94384


def fix_to_fix_speed_kn(df: pd.DataFrame, *, method: str = "geodesic") -> pd.Series:
    """Speed between consecutive fixes of the same vessel, in knots.

    `df` must carry mmsi, lat, lon and a timezone-aware timestamp. Rows are
    sorted internally; the returned Series is aligned to the input index.
    """
    required = {"mmsi", "lat", "lon", "timestamp"}
    if missing := required - set(df.columns):
        raise ValueError(f"missing columns: {sorted(missing)}")

    d = df.sort_values(["mmsi", "timestamp"])
    same = d["mmsi"].to_numpy()[1:] == d["mmsi"].to_numpy()[:-1]
    dt_s = np.diff(d["timestamp"].astype("int64").to_numpy()) / 1e9

    lat, lon = d["lat"].to_numpy(), d["lon"].to_numpy()
    if method == "geodesic":
        _, _, dist_m = GEOD.inv(lon[:-1], lat[:-1], lon[1:], lat[1:])
    elif method == "haversine":
        dist_m = _haversine_m(lat[:-1], lon[:-1], lat[1:], lon[1:])
    else:
        raise ValueError(f"unknown method {method!r}")

    with np.errstate(divide="ignore", invalid="ignore"):
        kn = (np.abs(dist_m) / dt_s) * MS_TO_KN
    kn[~same] = np.nan
    bad_dt = dt_s <= 0
    if bad_dt.any():
        logger.warning("%d consecutive pair(s) with a non-positive time interval — "
                       "duplicates or clock skew upstream", int(bad_dt.sum()))
        kn[bad_dt] = np.nan

    out = pd.Series(np.concatenate([[np.nan], kn]), index=d.index, name="speed_kn")
    return out.reindex(df.index)


def _haversine_m(lat1, lon1, lat2, lon2, radius_m: float = 6_371_008.8):
    p1, p2 = np.radians(lat1), np.radians(lat2)
    dp, dl = p2 - p1, np.radians(lon2 - lon1)
    a = np.sin(dp / 2) ** 2 + np.cos(p1) * np.cos(p2) * np.sin(dl / 2) ** 2
    return 2 * radius_m * np.arcsin(np.sqrt(a))

GEOD.inv accepts arrays, which is the whole reason to use pyproj rather than a per-row loop: a million fixes take about a second vectorised and several minutes row by row. The non-positive interval warning is not decorative — duplicate timestamps after deduplication, and clock skew between a receiver and an aggregator, both produce them regularly on a real feed.

2. Choose the method from the use, and record it

Use Method Why
Traffic density, routeing studies haversine 0.5% is far below the other uncertainties
Outlier rejection haversine the threshold is an order of magnitude above the error
Fuel and emissions modelling geodesic the effect being measured is often a few per cent
Regulatory speed compliance geodesic a disputed figure must not depend on the earth model
Anything crossing 70° latitude geodesic convergence makes spherical error less predictable

The recording matters as much as the choice. Two years of speeds computed by haversine with six months of geodesic spliced in has a discontinuity of a few hundredths of a knot at the join — invisible, and enough to produce a spurious trend in a fleet-efficiency study.

3. Cross-check against the reported speed

The transponder reports its own speed over ground, derived from the GNSS receiver’s Doppler solution rather than from position differences. It is an independent measurement, and comparing the two is nearly free:

def speed_residual_kn(df: pd.DataFrame) -> pd.Series:
    """Reported speed minus computed speed. Should be small and unbiased."""
    computed = fix_to_fix_speed_kn(df)
    resid = df["sog_kn"] - computed
    finite = resid[np.isfinite(resid)]
    if len(finite) > 100:
        logger.info("speed residual: median %.2f kn, p95 |%.2f| kn",
                    float(finite.median()), float(finite.abs().quantile(0.95)))
    return resid

A residual centred on zero with a spread of a few tenths of a knot is a healthy feed. A residual with a bias means one of the two is systematically wrong: a positive bias, where reported exceeds computed, usually means position fixes are being dropped so the computed speed averages over a longer, straighter path than the vessel took.

The Residual Between Two Independent Speeds Is a Free Diagnostic Three residual distributions of reported speed minus computed speed, in knots. The healthy feed is a narrow distribution centred on zero with a spread of a few tenths of a knot. The dropped-fixes feed is shifted positive by about half a knot, because missing intermediate fixes make the computed path straighter and therefore shorter than the vessel actually travelled, so the computed speed is too low. The clock-skew feed is centred on zero but far wider, because errors in the time denominator scatter the computed speed in both directions. A note records that the bias and the spread are independent signals pointing at different upstream faults. reported − computed (knots) 0 healthy dropped fixes biased positive: the computed path is straighter than the real one clock skew centred, but wide Bias and spread are independent signals: one points at the position stream, the other at the timestamps.

Verification and acceptance test 

import pytest


def test_a_known_leg_recovers_its_speed():
    """A 60 nm leg in exactly 4 hours is 15 knots, by definition."""
    lon1, lat1 = -75.0, 38.0
    lon2, lat2, _ = GEOD.fwd(lon1, lat1, 90.0, 60.0 * 1852.0)
    df = pd.DataFrame({
        "mmsi": [1, 1],
        "lat": [lat1, lat2], "lon": [lon1, lon2],
        "timestamp": pd.to_datetime(["2026-03-14T00:00Z", "2026-03-14T04:00Z"]),
    })
    kn = fix_to_fix_speed_kn(df).iloc[1]
    assert abs(kn - 15.0) < 0.01, f"recovered {kn:.3f} kn"


def test_euclidean_would_fail_the_same_test():
    """Guard: the wrong method must not pass, or the test proves nothing."""
    err = wrong_euclidean_km(60.0, 5.0, 60.0, 6.0) / (111.32 * np.cos(np.radians(60.0)))
    assert err > 1.9, "the Euclidean error at 60°N should be roughly a factor of two"


def test_duplicate_timestamps_do_not_produce_infinity(dupe_track):
    kn = fix_to_fix_speed_kn(dupe_track)
    assert np.isfinite(kn.dropna()).all(), "a zero time interval leaked an infinity"

The middle test is the fleet habit again: feeding the machinery input it must reject, so a check that quietly stopped working cannot pass by finding nothing.

Fix interval sets the ceiling on what speed can mean

Speed derived from two positions is an average over the interval between them, and the interval on an AIS feed varies by a factor of ninety. That makes “the speed” a quantity whose meaning changes across the same dataset, and no amount of care in the arithmetic fixes it.

At a two-second interval the derived speed is close to instantaneous and is dominated by position noise: a five-metre GNSS error over two seconds is 4.9 knots of spurious speed. At a three-minute interval the same error contributes 0.05 knots, and the derived speed is a genuine average — but it has smoothed away every manoeuvre shorter than three minutes.

def noise_speed_kn(position_sigma_m: float, interval_s: float) -> float:
    """Spurious speed produced by position error alone, at this fix interval."""
    if interval_s <= 0:
        raise ValueError("interval must be positive")
    return (position_sigma_m * np.sqrt(2) / interval_s) * MS_TO_KN


for dt in (2, 10, 60, 180):
    print(f"{dt:>4} s: {noise_speed_kn(5.0, dt):.2f} kn of noise")
#    2 s: 6.87 kn
#   10 s: 1.37 kn
#   60 s: 0.23 kn
#  180 s: 0.08 kn

Nearly seven knots of noise on a two-second interval is not a hypothetical: it is why raw fix-to-fix speed on a fast-reporting vessel is visibly noisier than the transponder’s own reported figure, which comes from a Doppler solution and does not have this problem.

The practical response is to derive speed over a distance rather than over a fixed number of fixes — accumulate consecutive segments until the vessel has moved far enough that position noise is a small fraction of the distance, then divide. That makes the averaging window adaptive: short in time for a fast vessel, long for a slow one, and constant in the quantity that actually governs the error.

At Two Seconds, Position Noise Is Seven Knots A curve of spurious speed in knots, caused by five metres of position error alone, against the fix interval in seconds. At two seconds the value is 6.9 knots, falling to 1.4 at ten seconds, 0.23 at sixty and 0.08 at one hundred and eighty. A shaded region marks intervals below about fifteen seconds, where the noise exceeds a half-knot analysis tolerance. An annotation notes that the fastest-reporting vessels are the ones with the shortest intervals, so the noisiest derived speeds belong to exactly the vessels a traffic study cares most about, and that averaging over a fixed distance rather than a fixed number of fixes makes the window adaptive. fix interval (seconds) noise (knots) 6.9 kn 0.5 kn analysis tolerance shaded: noise exceeds tolerance 21560180 The fastest-reporting vessels have the shortest intervals, so the noisiest speeds belong to the vessels of most interest.

Edge cases and gotchas

  • The antimeridian. A vessel crossing 180° has a longitude difference of about 360°, and every naive method produces a distance of half the earth’s circumference. Geod.inv handles it correctly; a hand-rolled haversine that does not normalise the longitude difference does not.
  • Timestamps must be timezone-aware. Mixing naive and aware datetimes raises in pandas; mixing two naive series in different local times does not, and produces intervals wrong by whole hours. Normalise to UTC at ingestion, as the timestamp drift guide sets out.
  • Speed and course must be computed from the same pair. Deriving speed from fixes i and i+1 while deriving course from i-1 and i+1 produces a kinematic record that is internally inconsistent, and any model relating the two will find structure that is an artefact of the offset.

Up: Speed and Heading Profiling for Maritime Analytics