Interpolating Gaps in AIS Vessel Trajectories

A reconstructed AIS track is only as honest as the space between its fixes. When a vessel sails out of terrestrial VHF range and reappears on a satellite pass tens of minutes later, a naive renderer draws a straight line between the last and first fix — a chord that can cut across a headland, slice over an island, or “teleport” the vessel dozens of kilometres in a single frame. This page is a focused fix for that artifact and belongs to the broader AIS Data Quality & Gap Filling workflow, which is the page to return to for the full quality-control context. Here we detect coverage gaps per MMSI, choose an interpolation method appropriate to the gap length, respect a maximum-gap policy that refuses to invent data across implausibly long gaps, and flag every synthetic fix so no downstream consumer mistakes it for an observation.

Root Cause: Why the Chord Is Wrong #

Coverage gaps have three ordinary causes, and none of them imply the vessel travelled in a straight line. Terrestrial receivers see roughly 40 nautical miles of VHF range, so any vessel beyond the coastal fringe simply is not heard. Satellite AIS fills the offshore gap but only on each orbital revisit, which can be tens of minutes apart at low latitudes. And any receiver can suffer downtime — power, antenna, or backhaul — leaving a hole even inshore.

The artifact appears the moment you connect the two fixes bracketing the gap with a Cartesian line in the projected or decimal-degree plane. That line is not the vessel’s path for two reasons: it ignores the great-circle geometry of a long over-water leg, and it ignores any turn the vessel made while unheard. Reproduce it with two fixes bracketing a wide gap around a headland:

import numpy as np
import pandas as pd

# Two real fixes bracketing a 42-minute coverage gap around a headland.
track = pd.DataFrame({
    "mmsi": [235099999, 235099999],
    "timestamp": pd.to_datetime(
        ["2026-03-01T09:00:00Z", "2026-03-01T09:42:00Z"], utc=True
    ),
    "lat": [50.20, 50.05],
    "lon": [-5.70, -5.00],   # straight chord cuts across Land's End
})

# Naive linear fill in lon/lat space: the vessel appears to sail over land.
dt = track["timestamp"].diff().dt.total_seconds().iloc[1]
print(f"gap = {dt/60:.0f} min, chord jumps "
      f"{track['lon'].iloc[1]-track['lon'].iloc[0]:.2f} deg lon in one step")

The chord above crosses the Land’s End peninsula. A great-circle path between the same endpoints bows differently, and if the vessel actually rounded the headland the true path is nothing like either straight line — which is exactly why long gaps must be flagged, not silently smoothed. Over-interpolation is its own hazard: a smooth curve laid over a multi-hour gap can hide a vessel that deliberately went dark, a signature the AIS spoofing and GPS jamming detectors are meant to catch.

Step-by-Step Fix #

Step 1 — Detect gaps by per-MMSI time delta #

Work one vessel at a time on a sorted, quality-controlled track. A gap is any inter-fix time delta above gap_threshold_s. Detecting on the delta — not on distance — is what makes the test robust to a vessel that stops: a stationary vessel with a long delta is still a gap in reporting, not a teleport.

from __future__ import annotations

import logging

import numpy as np
import pandas as pd

logger = logging.getLogger(__name__)


def detect_gaps(track: pd.DataFrame, gap_threshold_s: float) -> pd.DataFrame:
    """Label coverage gaps in a single-MMSI, time-sorted AIS track.

    Adds ``dt_s`` (inter-fix seconds) and ``is_gap`` columns. The gap flag
    marks the fix that *ends* a gap, i.e. the first fix after the silence.

    Raises:
        ValueError: if the track holds more than one MMSI.
    """
    if track["mmsi"].nunique() != 1:
        raise ValueError("detect_gaps expects a single-MMSI track")

    t = track.sort_values("timestamp").reset_index(drop=True)
    dt_s = t["timestamp"].diff().dt.total_seconds()
    t["dt_s"] = dt_s
    t["is_gap"] = (dt_s > gap_threshold_s).fillna(False)
    logger.info(
        "MMSI %d: %d fixes, %d coverage gaps > %.0f s",
        int(t["mmsi"].iloc[0]), len(t), int(t["is_gap"].sum()), gap_threshold_s,
    )
    return t

Step 2 — Choose the interpolation method by gap length #

Two regimes need two methods. Short gaps are best filled with a monotone cubic that respects the local curvature of the surrounding fixes so the reconstructed turn stays smooth; SciPy’s PchipInterpolator is shape-preserving and will not overshoot into an implausible loop. Long gaps span enough distance that Earth curvature dominates, so the honest fill is a great-circle path computed with pyproj’s geodesic solver. The Geod.npts method returns evenly spaced points along the true geodesic between two coordinates.

from pyproj import Geod
from scipy.interpolate import PchipInterpolator

_GEOD = Geod(ellps="WGS84")


def geodesic_fill(
    lat0: float, lon0: float, lat1: float, lon1: float, n_points: int
) -> list[tuple[float, float]]:
    """Return n_points evenly spaced (lat, lon) fixes along the great circle.

    Raises:
        ValueError: if n_points is negative.
    """
    if n_points < 0:
        raise ValueError(f"n_points must be non-negative; got {n_points}")
    if n_points == 0:
        return []
    pts = _GEOD.npts(lon0, lat0, lon1, lat1, n_points)  # returns (lon, lat)
    return [(lat, lon) for lon, lat in pts]


def pchip_fill(
    t_known: np.ndarray, lat: np.ndarray, lon: np.ndarray, t_query: np.ndarray
) -> np.ndarray:
    """Shape-preserving cubic fill for short gaps; returns (lat, lon) rows.

    Raises:
        ValueError: if fewer than two anchor fixes are supplied.
    """
    if t_known.size < 2:
        raise ValueError("PCHIP needs at least two anchor fixes")
    lat_i = PchipInterpolator(t_known, lat)(t_query)
    lon_i = PchipInterpolator(t_known, lon)(t_query)
    return np.column_stack([lat_i, lon_i])

Step 3 — Fill gaps, respecting the max-gap threshold, flag every synthetic fix #

The orchestrator walks each gap, picks the method by duration, refuses to fill anything longer than max_interp_gap_s, and stamps every inserted row interpolated=True. The cadence of synthetic fixes is set so an interpolated leg matches the vessel’s typical reporting interval.

from dataclasses import dataclass


@dataclass(frozen=True)
class InterpConfig:
    """Policy for gap interpolation."""

    gap_threshold_s: float = 600.0        # below this, not a gap
    short_gap_s: float = 1800.0           # PCHIP up to here; geodesic beyond
    max_interp_gap_s: float = 10_800.0    # 3 h: refuse to fill longer gaps
    synth_interval_s: float = 300.0       # cadence of inserted fixes


def interpolate_track(track: pd.DataFrame, cfg: InterpConfig) -> pd.DataFrame:
    """Fill coverage gaps in one vessel track, flagging synthetic fixes.

    Returns a track with an added boolean ``interpolated`` column. Gaps
    longer than ``max_interp_gap_s`` are left unfilled and remain flagged
    as real gaps for downstream awareness.
    """
    t = detect_gaps(track, cfg.gap_threshold_s)
    t["interpolated"] = False
    lat = t["lat"].to_numpy()
    lon = t["lon"].to_numpy()
    ts = t["timestamp"].astype("int64").to_numpy() / 1e9  # epoch seconds

    inserted: list[dict[str, object]] = []
    for i in np.flatnonzero(t["is_gap"].to_numpy()):
        dt = float(t["dt_s"].iloc[i])
        if dt > cfg.max_interp_gap_s:
            logger.warning(
                "MMSI %d: gap of %.0f s exceeds max %.0f s — left unfilled",
                int(t["mmsi"].iloc[0]), dt, cfg.max_interp_gap_s,
            )
            continue

        n_synth = max(int(dt // cfg.synth_interval_s) - 1, 0)
        if n_synth == 0:
            continue
        t_query = np.linspace(ts[i - 1], ts[i], n_synth + 2)[1:-1]

        if dt <= cfg.short_gap_s:
            anchors = slice(max(i - 2, 0), min(i + 2, len(t)))
            coords = pchip_fill(ts[anchors], lat[anchors], lon[anchors], t_query)
            synth = [(float(la), float(lo)) for la, lo in coords]
        else:
            synth = geodesic_fill(lat[i - 1], lon[i - 1], lat[i], lon[i], n_synth)

        for (la, lo), tq in zip(synth, t_query):
            inserted.append({
                "mmsi": int(t["mmsi"].iloc[0]),
                "timestamp": pd.Timestamp(tq, unit="s", tz="UTC"),
                "lat": la, "lon": lo, "interpolated": True,
            })

    if inserted:
        t = pd.concat([t, pd.DataFrame(inserted)], ignore_index=True)
        t = t.sort_values("timestamp").reset_index(drop=True)
    logger.info(
        "MMSI %d: inserted %d synthetic fixes", int(t["mmsi"].iloc[0]), len(inserted)
    )
    return t.drop(columns=["dt_s", "is_gap"], errors="ignore")

Verification and Acceptance Test #

An interpolated track is only acceptable if every synthetic leg is kinematically plausible and, when a land mask is available, stays on water. Assert both.

from pyproj import Geod

_GEOD_V = Geod(ellps="WGS84")


def verify_interpolation(
    track: pd.DataFrame,
    max_sog_kn: float = 60.0,
    land_mask=None,   # optional shapely (Multi)Polygon of land
) -> None:
    """Acceptance gate for an interpolated vessel track.

    Raises:
        AssertionError: if any synthetic leg is too fast or crosses land.
    """
    t = track.sort_values("timestamp").reset_index(drop=True)
    lat = t["lat"].to_numpy()
    lon = t["lon"].to_numpy()
    dt = t["timestamp"].diff().dt.total_seconds().to_numpy()

    _, _, dist = _GEOD_V.inv(lon[:-1], lat[:-1], lon[1:], lat[1:])
    with np.errstate(invalid="ignore", divide="ignore"):
        sog_kn = (dist / dt[1:]) * 3600.0 / 1852.0
    bad = np.nan_to_num(sog_kn, nan=0.0) > max_sog_kn
    if bad.any():
        raise AssertionError(
            f"{int(bad.sum())} interpolated legs exceed {max_sog_kn} kn"
        )

    if land_mask is not None:
        from shapely.geometry import Point
        synth = t.loc[t.get("interpolated", False)]
        on_land = [land_mask.contains(Point(lo, la))
                   for la, lo in zip(synth["lat"], synth["lon"])]
        if any(on_land):
            raise AssertionError(f"{sum(on_land)} interpolated fixes fall on land")
    logger.info("Interpolation verified: %d fixes, all legs plausible", len(t))

A green run guarantees the fill introduced no teleport and, if the mask was supplied, no over-land leg. The land test is a note rather than a solver here: true land-avoidance routing (bending the geodesic around a coastline) is a routing problem beyond simple interpolation, so the practical policy is to detect an over-land fill and either leave the gap unfilled or defer to a coastal routing graph.

Edge Cases and Gotchas #

• Interpolating across port stops. A vessel that berths, goes quiet, and later departs produces a long gap whose endpoints are the same position. Geodesic fill collapses to a stationary point (correct), but PCHIP over near-identical anchors can ring. Detect near-zero endpoint separation and skip interpolation entirely — the vessel did not move, so no synthetic track is warranted.

• Over-interpolation hiding dark-ship gaps. Filling every gap smooths away exactly the signal an investigator needs. Keep max_interp_gap_s conservative and always preserve the interpolated flag; a gap that a vessel created by switching off its transponder is evidence for detecting AIS spoofing and GPS jamming, and a fill must never overwrite it.

• Datum and CRS consistency. All geodesic math here assumes WGS84 decimal degrees (EPSG:4326), which is the AIS native datum. If upstream fixes were reprojected to a local metric CRS for another stage, convert them back to geographic coordinates before calling pyproj.Geod, or the great-circle distances will be silently wrong. Keep interpolation in geographic space and reproject only afterward.

---

Related #

• AIS Data Quality & Gap Filling — parent workflow: the QC stage that marks the gaps this page fills
• Segmenting Vessel Routes by Behavior — downstream stage that consumes gap-filled, flagged tracks
• Speed & Heading Profiling for Maritime Analytics — must exclude interpolated legs when deriving true kinematics
• Detecting AIS Spoofing & GPS Jamming — why over-interpolation can erase a dark-ship signal

Up: AIS Data Quality & Gap Filling