Vessel Density & Traffic Raster Generation #

Operational intent: this topic turns cleaned AIS tracks into gridded traffic products — the density surfaces that underlie routeing studies, marine spatial planning, collision-risk models and environmental impact assessments. It sits inside AIS Vessel Tracking & Route Automation, downstream of data quality and gap filling and of the behaviour segmentation that decides which vessels count. The computation is a rasterisation; the difficulty is that the obvious rasterisation measures the wrong thing.

Counting AIS messages per cell measures the transponders, not the traffic. A Class A vessel reports every two seconds while manoeuvring and every three minutes at anchor — a factor of ninety, driven entirely by what the vessel is doing and what equipment it carries. A message-count raster therefore lights up harbour approaches and turning basins because vessels report faster there, and under-represents open-water transit lanes where the same vessels report slowly. The pattern looks plausible. It is largely an artefact of the reporting schedule described in the data quality topic.

Reference Configuration #

The CRS row is the one most often wrong. A density raster built in Web Mercator has cells whose true ground area varies by a factor of several between the north and south of a regional domain, so a per-cell count is not a density at all. Build in an equal-area projection and reproject for display only.

From Fixes to Segments: the Unit of Credit #

The fix that makes the whole product defensible is to stop treating AIS positions as points and start treating consecutive pairs as segments. A segment carries a duration and a length, so a cell can be credited with the time a vessel spent inside it or the distance it covered there — quantities that are properties of the voyage rather than of the transponder.

"""Build creditable track segments from cleaned AIS fixes."""
from __future__ import annotations

import logging

import geopandas as gpd
import numpy as np
import pandas as pd
from shapely.geometry import LineString

logger = logging.getLogger(__name__)

MAX_GAP_S = 900.0        # 15 min
MAX_SPEED_MS = 20.6      # ~40 kn


def build_segments(fixes: gpd.GeoDataFrame) -> gpd.GeoDataFrame:
    """Consecutive-fix segments with duration and length, per vessel.

    `fixes` must be sorted by (mmsi, timestamp), carry a projected CRS, and
    have already been through gap detection and outlier rejection.
    """
    if fixes.crs is None or fixes.crs.is_geographic:
        raise ValueError("segment lengths require a projected, equal-area CRS")

    f = fixes.sort_values(["mmsi", "timestamp"]).reset_index(drop=True)
    same_vessel = f["mmsi"].to_numpy()[1:] == f["mmsi"].to_numpy()[:-1]
    dt = np.diff(f["timestamp"].astype("int64").to_numpy()) / 1e9
    x, y = f.geometry.x.to_numpy(), f.geometry.y.to_numpy()
    dist = np.hypot(np.diff(x), np.diff(y))

    keep = same_vessel & (dt > 0) & (dt <= MAX_GAP_S) & (dist / np.maximum(dt, 1e-9) <= MAX_SPEED_MS)
    dropped = int((same_vessel & ~keep).sum())
    if dropped:
        logger.info("dropped %d segment(s) failing gap/speed limits", dropped)

    idx = np.flatnonzero(keep)
    geom = [LineString([(x[i], y[i]), (x[i + 1], y[i + 1])]) for i in idx]
    return gpd.GeoDataFrame(
        {
            "mmsi": f["mmsi"].to_numpy()[idx],
            "duration_s": dt[idx],
            "length_m": dist[idx],
            "start": f["timestamp"].to_numpy()[idx],
        },
        geometry=geom, crs=fixes.crs,
    )

The gap and speed limits are not tuning knobs to be relaxed when too much is dropped. A segment spanning forty minutes is a straight line drawn across water nobody observed the vessel in, and crediting cells along it manufactures traffic. Where long gaps are common — a satellite-only feed, a poorly covered area — the honest response is a coarser grid and a longer aggregation period, not a longer permitted gap.

Rasterising by Intersection, Not by Endpoint #

Once segments exist, each one is split at cell boundaries and each cell credited with the piece that falls inside it. For a 100 m grid and 10-second fixes at transit speed, most segments cross one or two cells and the exact split matters little; for a 1 km grid or a satellite feed with minute-scale fixes, a segment can cross ten cells and endpoint-only crediting misses eight of them entirely, leaving a dotted line where a continuous track should be.

The cheap approximation that is usually good enough is to densify each segment to sub-cell spacing and credit the resulting points with a proportional share; the exact method rasterises the segment analytically. Densification is a few lines and is easy to reason about:

def accumulate_time(
    segments: gpd.GeoDataFrame, transform, shape: tuple[int, int], *, step_m: float
) -> np.ndarray:
    """Vessel-seconds per cell, by densifying each segment to `step_m`."""
    grid = np.zeros(shape, dtype="float64")
    inv = ~transform
    for geom, dur, length in zip(segments.geometry, segments.duration_s, segments.length_m):
        n = max(int(np.ceil(length / step_m)), 1)
        share = dur / n
        for k in range(n):
            px, py = geom.interpolate((k + 0.5) / n, normalized=True).coords[0]
            col, row = inv * (px, py)
            r, c = int(row), int(col)
            if 0 <= r < shape[0] and 0 <= c < shape[1]:
                grid[r, c] += share
    return grid

Set step_m to about a third of the cell size. Smaller wastes time; larger reintroduces the dotted-line artefact at exactly the scale the grid is meant to resolve.

Coverage Normalisation Is Not Optional #

A raw density raster is the product of two fields: how much traffic there was, and how well it was heard. The coverage geometry section of the parent overview explains why the second varies enormously across any real domain — terrestrial reception falls off at the horizon, satellite reception degrades where traffic is densest.

Dividing by a per-cell coverage layer — the fraction of the aggregation period during which that cell was within reception of a working receiver — turns an “observed traffic” raster into an estimate of actual traffic. It also makes the uncertainty explicit: cells with low coverage have their values inflated by a large factor and should be masked or flagged rather than published as confident estimates.

Grid or Kernel? Two Different Products #

There are two families of density product and they answer different questions. A binned grid credits each cell with what happened inside it: the value is an exact count of vessel-hours in that square, and summing the raster reproduces the fleet’s total time. A kernel density surface spreads each observation over a smooth bandwidth: the value is an estimate of an underlying intensity, the surface is continuous, and summing it reproduces the total only if the kernel is normalised and nothing left the domain.

The binned grid is the right default for anything that will be summed, compared or reported: total vessel-hours in a management zone, traffic through a proposed development area, exposure of a habitat polygon. Its values are conserved and auditable, and a zonal statistic over it is exact.

Kernel density is the right choice when the product is a map for human judgement and the binning artefacts are distracting — the checkerboard that appears when a lane runs diagonally across a coarse grid, or the speckle of a short aggregation period. It buys visual smoothness at the cost of introducing a bandwidth, which is a second scale parameter that has to be justified and published, and of moving traffic into cells where none occurred. A kernel surface should never be the input to a zonal sum.

The failure mode worth naming is producing a kernel surface because it looks better and then using it as if it were binned. Smoothed traffic crosses land, crosses closed areas, and appears inside exclusion zones the vessels never entered, all at low but non-zero intensity. Where such a product feeds a compliance question — “did any traffic enter the protected area” — the answer it gives is yes, always, everywhere, in proportion to the bandwidth.

Temporal Aggregation and the Comparability Trap #

Traffic products are almost always produced as a series — monthly grids, annual grids — and the series is usually the point. That makes comparability between periods more important than accuracy within any one of them, and comparability is fragile in ways that accuracy is not.

Four things drift between periods and each one produces a spurious trend:

The receiver network changes. A new base station comes online and the observed density in its footprint rises sharply, with no change in traffic. This is exactly what the coverage normalisation above corrects, and it only works if the coverage layer is recomputed per period rather than computed once and reused.

The fleet’s equipment changes. Class B carriage grows year on year in most jurisdictions, so a filter that includes Class B sees a rising vessel count that is partly regulatory. Splitting the product by class, and publishing both, makes the effect visible rather than mysterious.

The reporting schedule changes. Firmware updates and regulatory changes alter reporting intervals. A product built on message counts inherits every one of these; a product built on vessel-hours is immune, which is one more reason to prefer it.

Your own pipeline changes. A retuned outlier threshold, a corrected segment-gap limit, an upgraded coordinate transformation — each shifts the numbers. The defence is the reprocessing rule from the deployment section: when the pipeline changes, recompute the whole series with the new code and republish it as a new version, rather than appending a period produced by different code to a series produced by the old one.

That last rule is the expensive one and the one most often skipped, because recomputing five years of monthly grids is a real cost. It is also what makes a trend statement defensible, and the alternative is a series with an undocumented discontinuity at whatever date the deployment happened to land on.

Validation Gates #

Gate 1 — conservation. The sum of the raster, in vessel-hours, must equal the total observed duration of all retained segments, to within the rasterisation’s discretisation error. If it does not, segments are being dropped at the grid edge or credited twice.

Gate 2 — equal-area sanity. Total raster area must equal n_cells × cell_area to within a fraction of a per cent of the domain area. A conformal projection fails this immediately and obviously, which is the point.

Gate 3 — inter-period comparability. Recompute the previous period with the current code and compare. A change in the pipeline that alters a historical figure is a finding, not a nuisance: either the old number was wrong, or the new one is, and publishing a trend across the change without knowing which is not defensible.

Publishing an uncertainty band alongside the density #

A density value in a well-observed cell crossed by four hundred segments and a value in a cell crossed by two are not the same kind of number, and a single raster cannot say which is which. The cheapest useful companion band is the segment count per cell: it costs one extra accumulator in the same loop, it needs no statistical assumption, and it lets a consumer mask thin cells or weight by confidence.

Where a stronger statement is needed, a bootstrap over vessels — not over segments — gives an interval that respects the real unit of independence. Segments from one vessel are heavily correlated: a ship that transits a lane contributes forty consecutive segments that all say the same thing about that lane. Resampling vessels with replacement and recomputing the raster a few dozen times gives a per-cell spread that reflects how much of the value rests on a handful of ships. In coastal cells the interval is usually narrow; in offshore cells crossed by two vessels a month it is wide, and that width is the honest answer to “how much traffic is out there”.

Publish whichever of these the product can support, and say in the metadata which it is. A raster with no companion band is being asserted as uniformly reliable, which for a traffic product is never true.

Common Failure Modes #

Mixing Class A and Class B without saying so. Class B units are lower-powered and report less often, so their inclusion changes both the coverage field and the density. A product that silently changed its ship-type filter between years has a trend that is entirely an artefact.

Counting a moored vessel as traffic. A ship alongside for six days contributes 144 vessel-hours to one cell and swamps the berth’s neighbourhood. Whether that is wanted depends on the product — an emissions inventory wants it, a collision-risk model does not — and the behaviour segmentation stage is what makes the choice possible.

Interpolated positions credited as observations. The source column from the gap-filling stage must reach this stage. Crediting interpolated segments produces smooth, confident traffic through waters where nothing was observed.

Pipeline Integration and Downstream Handoff #

The output is a COG per period per quantity, in the equal-area CRS, with quantity, units, period_start, period_end, vessel_filter, coverage_normalised and min_coverage_published in the metadata, plus the coverage layer itself as a companion band. Downstream consumers — routeing studies, cumulative-impact models, port planning — join on the period, and the metadata is what allows two periods produced by different versions of this pipeline to be compared, or to be correctly refused.

Related #

• Building a Vessel Density Grid With GeoPandas — the segment-to-raster step end to end
• Interpolating Gaps in AIS Vessel Trajectories — why interpolated segments must not be credited as observations
• AIS Data Quality & Gap Filling — the stage that decides which segments may be credited
• Segmenting Vessel Routes by Behavior — how moored, fishing and transiting time are separated before aggregation
• Publishing COG & Zarr to Object Storage — how these rasters are published and versioned

Up: AIS Vessel Tracking & Route Automation