Segmenting Vessel Routes by Behavior #

Segmenting vessel routes by behavior is a deterministic, production-grade operation within the AIS Vessel Tracking & Route Automation pipeline that transforms continuous Automatic Identification System (AIS) telemetry into discrete, kinematically homogeneous track segments. Raw positional streams inherently mix transit, loitering, anchoring, and port maneuvering states, which corrupt downstream spatial indexing, regulatory compliance reporting, and habitat impact modeling. This workflow executes a memory-constrained, cloud-native pipeline that normalizes coordinate reference systems, derives vectorized kinematic derivatives, applies a finite state machine (FSM) for behavioral classification, and outputs standardized geospatial segments ready for spatial analysis and archival.

Reference Configuration and Specification #

Behavioral State Machine #

The pipeline evaluates each incoming ping against four primary behavioral states. The diagram below shows the allowed state transitions and their kinematic trigger conditions.

Solid arrows represent forward-progression transitions; dashed arrows represent recovery transitions back toward active navigation.

Memory-Constrained Ingestion and CRS Normalization #

Production AIS ingestion must handle fragmented payloads, inconsistent timestamp resolutions, and missing kinematic fields. Raw NMEA fixes should already be sanitized upstream by the AIS NMEA sentence parser so that malformed positional sentences never reach this stage. Raw data typically arrives as partitioned Parquet files with implicit EPSG:4326 coordinates. Distance-based behavioral thresholds require metric projections, so the pipeline reprojects to a locally appropriate UTM zone (or EPSG:3857 for global coverage) before computing inter-point distances and positional dispersion radii.

Memory constraints dictate out-of-core processing. Loading multi-terabyte AIS archives into a single GeoDataFrame triggers OOM failures on standard compute instances. Instead, the ingestion layer streams data via pyarrow.parquet.ParquetFile.iter_batches(), processes MMSI-grouped chunks sequentially, and flushes results to disk before advancing. This streaming architecture integrates with the upstream real-time AIS stream ingestion pipeline to guarantee temporal continuity, deduplicate retransmitted pings, and maintain strict ordering by BaseDateTime.

Coordinate transformation must be executed lazily and cached per spatial partition. Using pyproj with CRS.from_epsg() and Transformer.from_crs(always_xy=True) enforces thread-safe, repeatable projections and explicit longitude/latitude ordering per OGC standards. All distance calculations, buffer operations, and dispersion checks must occur in the projected metric space to prevent geodesic distortion artifacts.

Kinematic Derivation and Angular Unwrapping #

Behavioral classification cannot rely on raw coordinates alone. The segmentation engine derives speed over ground (SOG), course over ground (COG), and their temporal derivatives (ΔSOG, ΔCOG, acceleration) using rolling windows calibrated to vessel class metadata. COG discontinuities—for example, a transition from 359° to 1°—require angular unwrapping before differentiation to prevent false acceleration spikes.

Threshold calibration follows established maritime kinematic profiles. As documented in the speed and heading profiling workflow, transit states are defined by sustained SOG above 8 knots with heading variance below 15°, while loiter states exhibit SOG below 3 knots and positional dispersion exceeding 0.25 NM over a 30-minute window.

Derivative computation must account for irregular AIS transmission intervals. Instead of fixed-window rolling operations, the pipeline uses time-weighted differencing, normalizing every difference by the elapsed inter-ping interval:

$$\dot{v}_t = \frac{\text{SOG}_t - \text{SOG}_{t-1}}{\Delta t}, \qquad \Delta t = t - t_{t-1}\ \text{(seconds)}$$

Missing SOG/COG values are forward-filled only when the gap is less than 120 seconds; larger gaps trigger segment termination and UNKNOWN state assignment to preserve data integrity rather than fabricate kinematic continuity.

Memory-Constrained Python Implementation #

The implementation below demonstrates streaming ingestion, CRS projection, angular unwrapping, and FSM state assignment. It is structured for cloud execution with explicit chunking, vectorized operations, structured logging, and deterministic error handling.

import logging
import numpy as np
import pandas as pd
import pyarrow.parquet as pq
from pyproj import Transformer, CRS
from typing import Generator, Dict, Any, List

logger = logging.getLogger(__name__)

# Production constants — adjust per operational area and vessel class
TRANSIT_SOG_MIN: float = 8.0        # knots
LOITER_SOG_MAX: float  = 3.0        # knots
ANCHOR_RADIUS_M: float = 185.2      # 0.1 NM — anchoring dispersion threshold
COG_HYSTERESIS_TRANSIT: float = 15.0  # degrees — heading variance ceiling
STATE_TRANSITION_WINDOW: int = 3    # consecutive pings for state confirmation
ANCHOR_TRANSITION_WINDOW: int = 5   # stricter counter for ANCHORED
SOG_GAP_FILL_SECONDS: float = 120.0 # maximum gap before segment break
EPSG_WGS84: int  = 4326
EPSG_UTM_LOCAL: int = 32618         # WGS 84 / UTM Zone 18N — set per area


def unwrap_cog(degrees: np.ndarray) -> np.ndarray:
    """Angular unwrapping to prevent 359→1 discontinuities in differentiated COG."""
    return np.unwrap(np.deg2rad(degrees)) * (180.0 / np.pi)


def compute_kinematics(df: pd.DataFrame) -> pd.DataFrame:
    """
    Vectorized SOG/COG derivative computation with time-weighted differencing.

    Sorts by BaseDateTime, computes inter-ping time deltas, unwraps COG before
    differentiation, and flags gaps exceeding SOG_GAP_FILL_SECONDS as segment
    break points to prevent kinematic fabrication across data voids.
    """
    df = df.sort_values("BaseDateTime").reset_index(drop=True)
    df["dt_s"] = df["BaseDateTime"].diff().dt.total_seconds().fillna(0)

    # Mark large temporal gaps — downstream FSM forces UNKNOWN on these rows
    df["gap_break"] = df["dt_s"] > SOG_GAP_FILL_SECONDS

    df["sog_diff"] = df["SOG"].diff().fillna(0)
    df["cog_rad"]  = unwrap_cog(df["COG"].fillna(0).values)
    df["cog_diff"] = np.degrees(
        np.diff(df["cog_rad"], prepend=df["cog_rad"].iloc[0])
    )
    df["accel"] = df["sog_diff"] / df["dt_s"].replace(0, np.nan)
    df["heading_var"] = (
        df["cog_diff"]
        .rolling(window=STATE_TRANSITION_WINDOW, min_periods=1)
        .std()
    )
    return df


def assign_fsm_states(df: pd.DataFrame) -> pd.DataFrame:
    """
    Deterministic FSM state assignment with asymmetric hysteresis.

    Transitions require sustained threshold breaches across N consecutive
    valid pings (STATE_TRANSITION_WINDOW for TRANSIT, ANCHOR_TRANSITION_WINDOW
    for ANCHORED). Gap-break rows are immediately reset to UNKNOWN and restart
    the sustained counter, ensuring data voids never silently inherit a prior state.
    """
    states: List[str] = []
    current_state: str = "UNKNOWN"
    sustained_count: int = 0
    target_state: str = "UNKNOWN"

    for _, row in df.iterrows():
        # Temporal gap: reset state and counter
        if row.get("gap_break", False):
            current_state = "UNKNOWN"
            sustained_count = 0
            states.append(current_state)
            continue

        hv  = float(row.get("heading_var", 0) or 0)
        sog = float(row.get("SOG", 0) or 0)

        if sog >= TRANSIT_SOG_MIN and hv < COG_HYSTERESIS_TRANSIT:
            candidate = "TRANSIT"
            required  = STATE_TRANSITION_WINDOW
        elif sog < 1.0 and hv < 5.0:
            candidate = "ANCHORED"
            required  = ANCHOR_TRANSITION_WINDOW
        elif sog <= LOITER_SOG_MAX:
            candidate = "LOITER"
            required  = STATE_TRANSITION_WINDOW
        else:
            candidate = "MANEUVERING"
            required  = STATE_TRANSITION_WINDOW

        if candidate == target_state:
            sustained_count += 1
        else:
            target_state  = candidate
            sustained_count = 1

        if sustained_count >= required:
            current_state = target_state

        states.append(current_state)

    df["behavior_state"] = states
    return df


def process_ais_chunk(
    chunk: pd.DataFrame,
    transformer: Transformer,
) -> Generator[Dict[str, Any], None, None]:
    """
    Stream processing pipeline for a single MMSI chunk.

    Projects to metric CRS, computes kinematics, assigns FSM states, and
    yields one output record per homogeneous behavioral segment. Segments
    with fewer than two pings are discarded — they cannot form a valid LineString.
    """
    if chunk.empty or len(chunk) < 2:
        logger.debug("Skipping MMSI chunk with insufficient pings: %d", len(chunk))
        return

    lons = chunk["Longitude"].values
    lats = chunk["Latitude"].values
    xs, ys = transformer.transform(lons, lats)

    chunk = chunk.copy()
    chunk["x_m"] = xs
    chunk["y_m"] = ys

    chunk = compute_kinematics(chunk)
    chunk = assign_fsm_states(chunk)

    # Assign segment IDs on each behavioral state change
    state_changes      = chunk["behavior_state"].ne(chunk["behavior_state"].shift())
    chunk["segment_id"] = state_changes.cumsum()

    for sid, group in chunk.groupby("segment_id"):
        if len(group) < 2:
            continue

        dx = np.diff(group["x_m"].values)
        dy = np.diff(group["y_m"].values)
        track_len = float(np.sum(np.hypot(dx, dy)))

        wkt_coords = ", ".join(
            f"{x} {y}"
            for x, y in zip(group["Longitude"].values, group["Latitude"].values)
        )
        yield {
            "segment_id":       f"{group['MMSI'].iloc[0]}_{sid}",
            "mmsi":             int(group["MMSI"].iloc[0]),
            "start_ts":         group["BaseDateTime"].min(),
            "end_ts":           group["BaseDateTime"].max(),
            "behavior_state":   group["behavior_state"].mode().iloc[0],
            "track_length_m":   track_len,
            "duration_s":       float(
                (group["BaseDateTime"].max() - group["BaseDateTime"].min())
                .total_seconds()
            ),
            "mean_sog_kts":     float(group["SOG"].mean()),
            "cog_variance_deg": float(group["heading_var"].mean()),
            "geometry":         f"LINESTRING ({wkt_coords})",
        }


def run_segmentation_pipeline(input_path: str, output_path: str) -> None:
    """
    Entry point: streams Parquet in 500 k-row batches, groups by MMSI,
    and writes validated GeoParquet output with explicit CRS metadata.
    Raises ValueError on schema drift before any rows are written.
    """
    transformer = Transformer.from_crs(
        CRS.from_epsg(EPSG_WGS84),
        CRS.from_epsg(EPSG_UTM_LOCAL),
        always_xy=True,
    )

    parquet_file = pq.ParquetFile(input_path)
    results: List[Dict[str, Any]] = []

    for batch_idx, batch in enumerate(parquet_file.iter_batches(batch_size=500_000)):
        df = batch.to_pandas()
        logger.info("Processing batch %d: %d rows", batch_idx, len(df))

        for mmsi, group in df.groupby("MMSI"):
            results.extend(process_ais_chunk(group.copy(), transformer))

    if not results:
        raise ValueError(f"No valid segments produced from {input_path}")

    out_df = pd.DataFrame(results)
    out_df.to_parquet(output_path, index=False, engine="pyarrow")
    logger.info("Wrote %d segments to %s", len(out_df), output_path)

Key implementation decisions:

• always_xy=True enforces longitude/latitude ordering per OGC standards. The Transformer is instantiated once and reused across all MMSI groups to avoid per-chunk initialization overhead.
• The FSM uses separate counters for target and required thresholds. For high-frequency pings (>1 Hz), replace the row-level loop with a vectorized state-machine evaluation using numpy structured arrays or polars expressions for a 10–50× speedup.
• Port maneuvering can be isolated by intersecting projected coordinates with official port boundary polygons (IHO S-57 ENC layers). If a vessel remains within port limits and exhibits COG variance above 45°/min with SOG below 5 knots, override the FSM to MANEUVERING regardless of kinematic thresholds.

Validation Gates and Quality Control #

These checkpoints must pass before any segmented output is promoted to downstream stages:

Geometry validity check. Every geometry field must parse as a valid WKT LineString with at least two distinct coordinate pairs. Reject and log any null geometries or single-point degenerate segments.

import re

def validate_linestring(wkt: str) -> bool:
    """Returns True only when WKT contains at least two coordinate pairs."""
    coords = re.findall(r"-?\d+\.?\d*\s+-?\d+\.?\d*", wkt)
    return len(coords) >= 2

State distribution audit. After processing each Parquet batch, compute the proportion of pings assigned to each behavioral state. If UNKNOWN exceeds 15% of total pings, this typically indicates excessive temporal gaps in the source data or miscalibrated SOG thresholds for the vessel class being processed.

def audit_state_distribution(df: pd.DataFrame) -> None:
    dist = df["behavior_state"].value_counts(normalize=True)
    unknown_frac = dist.get("UNKNOWN", 0.0)
    if unknown_frac > 0.15:
        logger.warning(
            "UNKNOWN state fraction %.1f%% exceeds 15%% threshold — "
            "check source temporal gaps or SOG calibration",
            unknown_frac * 100,
        )

Schema drift detection. Before writing output Parquet, validate that the output DataFrame contains exactly the expected columns with correct dtypes. Raise ValueError on any schema mismatch to prevent silent corruption from propagating to downstream analytics.

Temporal monotonicity. Within each segment group, start_ts must be strictly less than end_ts. Zero-duration segments with equal timestamps indicate timestamp resolution problems in the source AIS feed and must be filtered before archival.

Track length floor. Segments with track_length_m below 10 m are almost certainly GPS jitter artifacts rather than genuine behavioral events. Flag and exclude these before handoff to the DBSCAN-based vessel track clustering stage.

Common Failure Modes and Diagnosis #

State flapping on noisy SOG feeds. Vessels operating near behavioral thresholds (e.g., repeatedly hovering between 7.5 and 8.5 knots) oscillate between TRANSIT and MANEUVERING every few pings when the hysteresis window is set too low. Symptom: the output contains hundreds of sub-minute segments for a single MMSI. Fix: increase STATE_TRANSITION_WINDOW from 3 to 5–7, or apply a 30-second median SOG smoothing pass before FSM evaluation.

Angular unwrapping failure on sparse pings. When pings are separated by more than 5–10 minutes, np.unwrap() may produce incorrect phase corrections because the vessel may have completed more than half a rotation between observations. Symptom: cog_diff values exceeding ±180° that trigger MANEUVERING on vessels in open-water transit. Fix: reset COG unwrapping at each gap-break boundary (already enforced by the gap_break flag in compute_kinematics). Verify that all gap-break rows yield UNKNOWN state before re-evaluation.

OOM on large MMSI groups. A single high-frequency fishing vessel transmitting at 2-second intervals over 12 hours generates ~21,600 rows. If many such vessels land in a single Parquet batch, the in-memory df.groupby("MMSI") operation can exceed available RAM. Fix: pre-partition input Parquet by MMSI using pyarrow.dataset so that each worker processes only a bounded per-vessel slice.

Port polygon intersection failures. When IHO S-57 ENC boundary polygons are loaded in geographic CRS (EPSG:4326) but vessel positions are projected (EPSG:32618), the spatial intersection silently produces empty results rather than raising a CRS mismatch error. Fix: always reproject port polygons to the working metric CRS before any spatial predicate operation; verify with assert port_gdf.crs == vessel_crs at pipeline startup.

Output Standardization and Downstream Handoff #

Segmented outputs are serialized to GeoParquet with explicit CRS metadata, schema validation, and partitioning by mmsi and segment_date. Each segment record contains: segment_id, mmsi, start_ts, end_ts, behavior_state, track_length_m, duration_s, mean_sog_kts, cog_variance_deg, and geometry (WKT LineString). The schema enforces strict typing and rejects null geometries or malformed timestamps before write.

Standardized segments feed directly into the DBSCAN-based vessel track clustering stage for automated shipping lane delineation, habitat corridor mapping, and regulatory compliance auditing. The deterministic segmentation ensures downstream algorithms operate on kinematically coherent trajectories rather than mixed-state noise.

A metadata manifest accompanies each GeoParquet output, recording: source Parquet path and SHA-256 checksum, processing timestamp, EPSG codes for input and working CRS, FSM threshold configuration, state distribution audit results, and total segment count. This lineage record is required for anomaly detection in AIS trajectories, which relies on known-clean behavioral segments as its baseline reference population.

All outputs must pass schema validation against OGC GeoParquet specifications before archival. Implement checksum verification and pipeline rollback triggers on schema drift so that a corrupted batch does not silently propagate through the analytics stack.

---

Related

• AIS Vessel Tracking & Route Automation — parent topic overview
• Real-Time AIS Stream Ingestion Pipelines — upstream ingestion stage
• Speed and Heading Profiling for Maritime Analytics — kinematic threshold calibration
• Clustering Vessel Tracks with DBSCAN — downstream spatial analysis
• Anomaly Detection in AIS Trajectories — behavioral baseline consumer