Speed and Heading Profiling for Maritime Analytics

Speed and heading profiling is the foundational kinematic layer within the AIS Vessel Tracking & Route Automation pipeline: it transforms temporally irregular Automatic Identification System positional fixes into validated speed-over-ground (SOG) and course-over-ground (COG) trajectories that every downstream analytical stage depends on. This workflow is deterministic and compute-bound — the engineering problem is not modelling uncertainty, but eliminating it: ingest high-volume AIS telemetry, compute geodesically accurate kinematic derivatives on native WGS84 coordinates, apply boolean quality masks, and emit partitioned cloud-ready profiles without breaching memory budgets or accumulating projection-induced distortion.

Reference Configuration and Specification #

The CRS decision deserves emphasis: reprojecting AIS fixes to a planar UTM zone before computing speed and heading introduces seam discontinuities at zone boundaries, adds CRS transformation latency, and degrades accuracy at high latitudes. Computing directly on native WGS84 with pyproj.Geod eliminates all three pathologies. For background on CRS selection in coastal workflows, see CRS Alignment for Coastal GIS Projects.

Pipeline Architecture #

The diagram below shows the four-stage data flow from raw NMEA decode through to partitioned kinematic profiles:

Operational Schema and Ingestion Architecture #

Raw AIS streams arrive as decoded NMEA payloads. Before any kinematic computation, the ingestion layer must enforce a strict typed schema. The minimal viable schema includes mmsi (UInt32), timestamp (datetime with UTC timezone), lat/lon (Float64), sog (Float32), cog (Float32), and nav_status (UInt8). Missing or malformed fields must be coerced or raised as errors — never silently imputed — to preserve auditability across regulatory audits.

Full-trajectory loading into in-memory GeoDataFrames is unsustainable at regional or basin scale. Ingestion must use columnar, memory-mapped formats: Apache Parquet with Polars lazy frames provides the necessary throughput and deferred evaluation. When connecting to an upstream real-time AIS stream ingestion pipeline, the profiling layer should consume micro-batches of 5–15 minutes partitioned by mmsi and temporal index, guaranteeing bounded memory usage and strict temporal ordering for derivative calculations.

import logging
from pathlib import Path

import polars as pl

logger = logging.getLogger(__name__)

AIS_SCHEMA: dict[str, pl.PolarsDataType] = {
    "mmsi":       pl.UInt32,
    "timestamp":  pl.Datetime("us", "UTC"),
    "lat":        pl.Float64,
    "lon":        pl.Float64,
    "sog":        pl.Float32,   # AIS-reported SOG, knots
    "cog":        pl.Float32,   # AIS-reported COG, degrees true
    "nav_status": pl.UInt8,
}

REQUIRED_COLS = set(AIS_SCHEMA)


def load_ais_batch(parquet_path: Path) -> pl.LazyFrame:
    """
    Open a single Parquet micro-batch with schema validation.
    Raises ValueError if required columns are absent.
    """
    lf = pl.scan_parquet(parquet_path)
    present = set(lf.columns)
    missing = REQUIRED_COLS - present
    if missing:
        raise ValueError(
            f"AIS schema violation in {parquet_path}: missing columns {missing}"
        )
    return lf.cast(AIS_SCHEMA)  # type: ignore[arg-type]

Vectorized Geodesic Kinematics #

Kinematic profiling requires precise computation of inter-fix distance, elapsed time, and directional change. Row-wise iteration is unacceptable in production at basin scale; all operations must be vectorized across the entire column array. The implementation below uses pyproj.Geod for WGS84 orthodromic distance and forward azimuth, numpy for shift-and-mask operations, and polars for columnar execution. Operating on native WGS84 coordinates bypasses geometry-object instantiation overhead and maintains sub-metre positional fidelity at all latitudes — an approach consistent with the AIS NMEA sentence parsing standards that supply the upstream fixes.

import logging
from typing import Final

import numpy as np
import polars as pl
from pyproj import Geod

logger = logging.getLogger(__name__)

KNOTS_PER_MPS: Final[float] = 1.0 / 0.514444   # 1 m/s → knots
GEOD: Final[Geod] = Geod(ellps="WGS84")         # thread-safe, stateless


def compute_kinematics(df: pl.DataFrame) -> pl.DataFrame:
    """
    Vectorized SOG/COG computation for one MMSI micro-batch.

    The input DataFrame MUST be sorted by (mmsi, timestamp) before
    calling this function. group_by(..., maintain_order=True) satisfies
    this contract when applied after .sort(["mmsi", "timestamp"]).

    Returns the original frame augmented with:
      sog_calc  – computed speed-over-ground, knots (Float32)
      cog_calc  – computed course-over-ground, degrees true [0, 360) (Float32)
      dt_sec    – elapsed seconds between consecutive fixes (Float32)
    """
    lat_arr: np.ndarray = df["lat"].to_numpy()
    lon_arr: np.ndarray = df["lon"].to_numpy()
    # Polars Datetime cast to Int64 yields microseconds since epoch
    ts_us:   np.ndarray = df["timestamp"].cast(pl.Int64).to_numpy()

    n = len(lat_arr)
    if n < 2:
        # Single-fix trajectory: kinematics are undefined
        nan_col = np.full(n, np.nan, dtype=np.float32)
        return df.with_columns([
            pl.Series("sog_calc", nan_col, dtype=pl.Float32),
            pl.Series("cog_calc", nan_col, dtype=pl.Float32),
            pl.Series("dt_sec",   nan_col, dtype=pl.Float32),
        ])

    # Build shifted predecessor arrays (index 0 is always NaN/undefined)
    lat_prev = np.empty(n, dtype=np.float64)
    lon_prev = np.empty(n, dtype=np.float64)
    ts_prev  = np.empty(n, dtype=np.int64)

    lat_prev[0] = np.nan
    lon_prev[0] = np.nan
    ts_prev[0]  = ts_us[0]
    lat_prev[1:] = lat_arr[:-1]
    lon_prev[1:] = lon_arr[:-1]
    ts_prev[1:]  = ts_us[:-1]

    # Temporal delta: microseconds → seconds
    dt_sec = (ts_us - ts_prev) / 1_000_000.0

    # Geodesic inverse: forward azimuth (°) and ellipsoidal distance (m)
    # pyproj.Geod.inv accepts (lon1, lat1, lon2, lat2) in degrees
    fwd_az, _, dist_m = GEOD.inv(lon_prev, lat_prev, lon_arr, lat_arr)

    # Guard against zero/negative time deltas (duplicate timestamps, clock skew)
    dt_safe = np.where(dt_sec > 0.0, dt_sec, np.nan)

    # SOG in knots
    sog_calc = (dist_m / dt_safe) * KNOTS_PER_MPS

    # COG normalized to [0, 360) — numpy mod handles negative azimuths correctly
    cog_calc = np.mod(fwd_az, 360.0)

    # Null out the seed row and any invalid delta rows
    valid = dt_sec > 0.0
    sog_out = np.where(valid, sog_calc, np.nan).astype(np.float32)
    cog_out = np.where(valid, cog_calc, np.nan).astype(np.float32)
    dt_out  = dt_sec.astype(np.float32)

    logger.debug(
        "compute_kinematics: mmsi=%s  fixes=%d  valid_rows=%d",
        df["mmsi"][0], n, int(valid.sum()),
    )

    return df.with_columns([
        pl.Series("sog_calc", sog_out, dtype=pl.Float32),
        pl.Series("cog_calc", cog_out, dtype=pl.Float32),
        pl.Series("dt_sec",   dt_out,  dtype=pl.Float32),
    ])


def profile_batch(parquet_path: "Path") -> pl.DataFrame:
    """
    End-to-end kinematic profiling for one Parquet micro-batch.
    Suitable for direct invocation inside a Dask or multiprocessing worker.
    """
    from pathlib import Path  # noqa: PLC0415
    lf = pl.scan_parquet(Path(parquet_path))
    df = (
        lf
        .sort(["mmsi", "timestamp"])
        .collect()
    )
    return (
        df
        .group_by("mmsi", maintain_order=True)
        .map_groups(compute_kinematics)
    )

Validation Gates and Quality Control #

Raw AIS telemetry contains GPS drift, multipath errors, and transmission artefacts that corrupt kinematic derivatives. Quality control applies four deterministic boolean masks; none introduce stochastic behaviour, ensuring full reproducibility across pipeline runs.

Stationary suppression #

Vessels at anchor or mooring exhibit GPS jitter that generates false SOG/COG spikes. Apply a compound filter: sog_calc < 0.5 and nav_status == 1 (at anchor per ITU-R M.1371). Kinematic updates during stationary periods are nullified, not dropped, so the record remains in the audit trail.

def apply_stationary_mask(df: pl.DataFrame) -> pl.DataFrame:
    stationary = (df["sog_calc"] < 0.5) & (df["nav_status"] == 1)
    return df.with_columns([
        pl.when(stationary).then(None).otherwise(pl.col("sog_calc")).alias("sog_calc"),
        pl.when(stationary).then(None).otherwise(pl.col("cog_calc")).alias("cog_calc"),
    ])

Temporal gap enforcement #

AIS dropout or satellite handoff creates discontinuous trajectories. Any fix where dt_sec > 3600 terminates the current trajectory segment; all kinematic fields for that row are set to null, and a new segment identifier (seg_id) is incremented. Downstream behavioral route segmentation relies on these segment boundaries to avoid computing metrics across multi-hour gaps.

def assign_segment_ids(df: pl.DataFrame) -> pl.DataFrame:
    gap = (df["dt_sec"] > 3600) | df["dt_sec"].is_null()
    seg_id = gap.cast(pl.UInt32).cum_sum()
    return df.with_columns(seg_id.alias("seg_id"))

Heading rate-of-change rejection #

Commercial vessels cannot execute instantaneous course reversals. A rolling check rejects fixes where |cog_calc − lag(cog_calc, 1)| > 45° over a 5-minute window — violations indicate GPS multipath or sensor fault. Rejected rows receive q_flag = 2.

SOG cross-validation #

The computed SOG from geodesic distance is cross-checked against the AIS-reported sog field. Discrepancies exceeding 20 % of the reported value (accounting for GNSS positioning latency and encoding quantisation) trigger q_flag = 1 for downstream anomaly review. This validation catches GNSS receiver clock errors and erroneous NMEA sentence re-broadcasts.

The single most common quality-gate defect is the circular wrap-around in the heading-rate check. Course-over-ground is an angle on a circle, so the arithmetic difference between two bearings is not the angular distance between them. A vessel steering steadily through true north reports cog_calc values such as 358° then 002° — an arithmetic delta of 356°, but a true turn of only 4°. The wrap-aware reduction collapses the difference onto the shorter arc, so a genuine 4° correction is never misclassified as an implausible 356° reversal. Formally, the angular distance between two true bearings $\theta_1$ and $\theta_2$ is:

$$\Delta\theta = \min\bigl(\,|\theta_1 - \theta_2|,\; 360^\circ - |\theta_1 - \theta_2|\,\bigr)$$

The heading-rate gate then rejects a fix only when $\Delta\theta > 45^\circ$ within the 5-minute window:

import polars as pl


def apply_quality_flags(df: pl.DataFrame) -> pl.DataFrame:
    """
    Assign q_flag:
      0 = clean
      1 = SOG cross-validation failure (>20% deviation)
      2 = implausible heading rate-of-change (>45° / fix)
    """
    sog_deviation = (
        (df["sog_calc"] - df["sog"]).abs() / (df["sog"].abs() + 1e-6)
    )
    sog_fail = sog_deviation > 0.20

    cog_delta = (df["cog_calc"] - df["cog_calc"].shift(1)).abs()
    # Wrap-around: 359° → 1° is a 2° change, not 358°
    cog_delta_wrapped = pl.min_horizontal(cog_delta, 360.0 - cog_delta)
    heading_fail = cog_delta_wrapped > 45.0

    q_flag = (
        pl.when(heading_fail).then(pl.lit(2))
        .when(sog_fail).then(pl.lit(1))
        .otherwise(pl.lit(0))
        .cast(pl.UInt8)
        .alias("q_flag")
    )
    return df.with_columns(q_flag)

Common Failure Modes and Diagnosis #

Phantom speed spikes at trajectory start. The first row of each MMSI group always has lat_prev = NaN; if the NaN guard is absent, pyproj.Geod.inv returns inf or raises. Fix: confirm that valid = dt_sec > 0.0 suppresses row 0, and that the single-fix branch returns an all-NaN kinematic output.

SOG in m/s instead of knots. Omitting the KNOTS_PER_MPS conversion factor produces SOG values approximately 1.94× too high. The cross-validation gate catches this: virtually every row will fail the ±20 % tolerance against AIS-reported SOG. Fix: verify the unit constant at the top of compute_kinematics.

COG wrap-around artefacts. When a vessel crosses true north (359° → 1°), naive subtraction yields Δ = 358° and triggers the heading-rate rejection filter. Fix: apply the wrap-aware delta min(|Δ|, 360 − |Δ|), as shown in apply_quality_flags above.

Memory overflow on basin-scale Parquet scans. Calling .collect() on an un-partitioned basin-scale lazy frame materialises hundreds of gigabytes. Fix: partition input by MMSI prefix or date before scanning, and use pl.scan_parquet with predicate pushdown (filter before collect). Confirm the execution plan with lf.explain(optimized=True) before running.

Pipeline Integration and Downstream Handoff #

Once kinematic profiles pass quality validation, they are serialised to cloud-optimised Parquet with polars.LazyFrame.sink_parquet, which streams results directly to object storage without materialising intermediate DataFrames. The output schema appends sog_calc, cog_calc, dt_sec, seg_id, and q_flag to the original record set.

import logging
from pathlib import Path

import polars as pl

logger = logging.getLogger(__name__)


def write_kinematic_profiles(
    profiled: pl.DataFrame,
    output_root: Path,
) -> None:
    """
    Write validated kinematic profiles to Hive-partitioned Parquet.
    Partition key: mmsi / year / month / day.
    Streaming sink keeps memory footprint constant for basin-scale datasets.
    """
    profiled = profiled.with_columns([
        pl.col("timestamp").dt.year().alias("year"),
        pl.col("timestamp").dt.month().alias("month"),
        pl.col("timestamp").dt.day().alias("day"),
    ])

    output_root.mkdir(parents=True, exist_ok=True)
    lf = profiled.lazy()

    lf.sink_parquet(
        output_root / "part-{mmsi}-{year}-{month}-{day}.parquet",
        maintain_order=False,   # permit parallel partition writes
    )
    logger.info("Kinematic profiles written to %s", output_root)

Every output partition carries a small sidecar manifest so downstream consumers can validate provenance before reading: the source micro-batch URI, the geodesic ellipsoid (WGS84), the conversion constant KNOTS_PER_MPS, the active threshold set (gap limit, stationary cut-off, heading-rate ceiling, SOG tolerance), and a profiler version string. Pinning the threshold set into the manifest makes a re-run reproducible even after defaults change, and lets the segmentation and anomaly stages reject profiles produced under incompatible filter parameters rather than silently mixing them.

import json
from pathlib import Path

PROFILER_VERSION = "1.4.0"


def write_partition_manifest(output_root: Path, source_uri: str) -> None:
    """Persist a provenance manifest beside the profiled partitions."""
    manifest = {
        "profiler_version": PROFILER_VERSION,
        "source_uri": source_uri,
        "ellipsoid": "WGS84",
        "knots_per_mps": 1.0 / 0.514444,
        "thresholds": {
            "temporal_gap_sec": 3600,
            "stationary_sog_kn": 0.5,
            "heading_rate_deg": 45.0,
            "sog_tolerance_frac": 0.20,
        },
    }
    (output_root / "_manifest.json").write_text(json.dumps(manifest, indent=2))

Raw telemetry must be retained alongside profiled outputs for forensic auditing; schema evolution should be managed with explicit Parquet metadata and a version field to support backward-compatible reads. The partitioned kinematic profiles feed directly into the behavioral route segmentation stage, where SOG and COG signatures are mapped to operational states (transit, maneuvering, fishing, drifting). They also provide the kinematic baseline for anomaly detection in AIS trajectories, where persistent q_flag patterns indicate either sensor malfunction or intentional position spoofing. The same seg_id boundaries computed here are reused as the temporal index for any real-time AIS stream ingestion replay, so segment numbering must remain stable across reprocessing.

---

Related #

• AIS Vessel Tracking & Route Automation — parent pipeline overview
• Real-Time AIS Stream Ingestion Pipelines — upstream micro-batch source
• Segmenting Vessel Routes by Behavior — consumes kinematic profiles
• Anomaly Detection in AIS Trajectories — downstream quality-flag consumer
• CRS Alignment for Coastal GIS Projects — WGS84 datum decisions