Coastal & Marine Spatial Analysis Pipelines

Building a Vessel Density Grid With GeoPandas

This is the end-to-end implementation of the product the vessel density topic specifies: a raster whose cell values are vessel-hours, in an equal-area projection, normalised for reception coverage, and carrying enough metadata to be compared with next month’s. Every step below exists because its absence produces a raster that looks right and means something else.

Root cause: the obvious implementation counts transponders

gdf.geometry from an AIS feed is a set of points, and rasterising points counts messages. As the parent topic sets out, message rate varies by a factor of ninety with vessel state, so the resulting map is a picture of reporting behaviour with traffic as a confounder. Everything here follows from replacing points with segments.

import geopandas as gpd

# The one-liner that produces the wrong product:
counts = gpd.sjoin(fixes, grid, predicate="within").groupby("index_right").size()

That is three lines shorter than the correct version and answers a different question.

Counting Endpoints Leaves Holes in Every Track A vessel track crosses a row of five grid cells diagonally. Its reported fixes are unevenly spaced, with four clustered in the first cell where the vessel was manoeuvring slowly and only two spread across the remaining four cells during a fast transit. Under endpoint counting the first cell receives four credits, one further cell receives two, and three cells receive nothing at all despite the vessel having crossed them, leaving a dotted line. Under intersection-length accumulation each cell is credited with the duration of the track segment falling inside it, so all five cells receive a value and the fast cells receive proportionally less than the slow one. 4 fixes 1 fix 1 fix none none endpoint counting 3 of 5 crossed cells get nothing intersection length every crossed cell credited, in proportion to time spent the slow cell earns most of the time, and the fast ones still earn some The dotted line is what a point-based density map looks like at any resolution finer than the fix spacing.

Step-by-step implementation

1. Equal-area CRS, asserted rather than assumed 

import logging

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

logger = logging.getLogger(__name__)

EQUAL_AREA = "EPSG:3035"          # Europe; use 6933 or a local Albers elsewhere


def to_equal_area(fixes: gpd.GeoDataFrame, crs: str = EQUAL_AREA) -> gpd.GeoDataFrame:
    """Reproject and confirm the target really is equal-area."""
    from pyproj import CRS
    target = CRS.from_user_input(crs)
    if target.is_geographic:
        raise ValueError(f"{crs} is geographic; cell areas would vary across the grid")
    name = (target.coordinate_operation.method_name or "").lower()
    if not any(k in name for k in ("equal area", "albers", "lambert azimuthal")):
        logger.warning("%s may not be equal-area (%s) — per-cell density will be biased",
                       crs, target.coordinate_operation.method_name)
    return fixes.to_crs(target)

The warning rather than a hard failure is deliberate: some legitimate national grids are equal-area under names the check does not recognise, and a false failure would push people to remove the check entirely.

2. Segments, filtered 

MAX_GAP_S = 900.0
MAX_SPEED_MS = 20.6


def build_segments(fixes: gpd.GeoDataFrame) -> gpd.GeoDataFrame:
    """Consecutive-fix segments carrying duration and length."""
    f = fixes.sort_values(["mmsi", "timestamp"]).reset_index(drop=True)
    same = 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))

    observed = f.get("source", pd.Series("observed", index=f.index)).to_numpy()
    both_observed = (observed[1:] == "observed") & (observed[:-1] == "observed")

    keep = same & both_observed & (dt > 0) & (dt <= MAX_GAP_S) & \
           (dist / np.maximum(dt, 1e-9) <= MAX_SPEED_MS)
    logger.info("kept %d of %d candidate segments", int(keep.sum()), int(same.sum()))

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

The both_observed condition is the one most implementations omit. If the gap-filling stage inserted interpolated positions, a segment between two of them is a straight line through water nobody observed, and crediting it manufactures traffic exactly where the evidence is weakest.

3. Accumulate by densified intersection 

def accumulate_vessel_seconds(segments: gpd.GeoDataFrame, transform, shape,
                              *, step_m: float | None = None) -> np.ndarray:
    """Vessel-seconds per cell, by densifying each segment to sub-cell spacing."""
    cell = abs(transform.a)
    step = step_m or cell / 3.0
    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)), 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
    logger.info("accumulated %.1f vessel-hours over %d populated cell(s)",
                grid.sum() / 3600.0, int((grid > 0).sum()))
    return grid

Setting step to a third of the cell size is the balance point: finer wastes time without changing the result, coarser reintroduces the dotted-line artefact at exactly the scale the grid resolves.

4. Normalise, then write the units into the file 

import rasterio


def write_density_cog(grid_s: np.ndarray, coverage: np.ndarray, transform, crs,
                      path: str, *, period_start: str, period_end: str,
                      vessel_filter: str, min_coverage: float = 0.3) -> None:
    """Coverage-normalised vessel-hours, as a COG that states what it is."""
    if grid_s.shape != coverage.shape:
        raise ValueError("grid and coverage layer must have the same shape")
    with np.errstate(invalid="ignore", divide="ignore"):
        hours = (grid_s / 3600.0) / coverage
    hours[coverage < min_coverage] = np.nan

    profile = {
        "driver": "COG", "dtype": "float32", "nodata": np.nan, "count": 2,
        "height": hours.shape[0], "width": hours.shape[1],
        "transform": transform, "crs": crs, "compress": "DEFLATE", "blocksize": 512,
    }
    with rasterio.open(path, "w", **profile) as dst:
        dst.write(hours.astype("float32"), 1)
        dst.write(coverage.astype("float32"), 2)
        dst.set_band_description(1, "vessel_hours_per_cell")
        dst.set_band_description(2, "reception_coverage_fraction")
        dst.update_tags(
            quantity="vessel-hours", units="h", period_start=period_start,
            period_end=period_end, vessel_filter=vessel_filter,
            coverage_normalised="true", min_coverage_published=str(min_coverage),
        )
    logger.info("wrote %s (%s to %s, filter %s)", path, period_start, period_end, vessel_filter)

Two bands rather than one, and six tags. The coverage band lets a consumer re-derive the raw counts or apply a stricter mask; the tags are what make two months comparable. A file without them is a picture.

Five Stages, and What Each One Guarantees A left-to-right pipeline of five stages with record counts. Fixes, 41 million, in a geographic CRS. Reprojected to an equal-area CRS, still 41 million, establishing that per-cell area is constant. Segments, 38.6 million after rejecting long gaps, impossible speeds and interpolated positions, establishing that every credited unit is an observation. Accumulated, giving 2.1 million populated cells holding vessel-seconds, establishing that credit follows time rather than message count. Normalised and written, a two-band cloud-optimized GeoTIFF carrying vessel-hours and the coverage fraction, plus tags for quantity, units, period, vessel filter and the minimum published coverage, establishing that the product can be compared with another month. fixes 41.0 M equal-area cell area constant segments 38.6 M observed accumulate 2.1 M cells 2-band COG hours + coverage 6 metadata tags Each stage establishes one property, and skipping any of them produces a raster that renders identically and answers a different question. The last stage is the one most often skipped, and it is the one that decides whether the product can be compared with the same product from a different month or a different team. 2.4 M segments rejected: long gaps, impossible speeds, interpolated positions.

Verification and acceptance test 

import pytest


def test_the_raster_conserves_the_observed_time(segments, grid_s):
    total_h = float(segments["duration_s"].sum()) / 3600.0
    grid_h = float(grid_s.sum()) / 3600.0
    assert abs(grid_h - total_h) / total_h < 0.01, (
        f"grid holds {grid_h:.0f} vessel-hours but segments total {total_h:.0f} — "
        "segments are being dropped at the grid edge or credited twice")


def test_cells_are_equal_area(transform, crs):
    from pyproj import CRS
    assert not CRS.from_user_input(crs).is_geographic
    assert abs(abs(transform.a) - abs(transform.e)) < 1e-6, "cells are not square"


def test_the_product_states_what_it_is(path):
    import rasterio
    with rasterio.open(path) as src:
        for tag in ("quantity", "units", "period_start", "vessel_filter"):
            assert src.tags().get(tag), f"missing required tag: {tag}"

The conservation test is the one that catches almost every implementation bug at once: a densification error, an off-by-one in the cell index, a segment credited to two cells, or a filter applied after accumulation instead of before.

Scaling the accumulation beyond one machine

The loop above is correct and single-threaded, and a month of a regional feed is tens of millions of segments. Three changes take it to production volume without changing what it computes.

Vectorise the densification. Instead of interpolating point by point, generate all sample points for all segments at once and accumulate with np.bincount on the flattened cell index. The result is identical and roughly two orders of magnitude faster:

def accumulate_vectorised(segments: gpd.GeoDataFrame, transform, shape,
                          *, step_m: float) -> np.ndarray:
    coords = np.array([g.coords for g in segments.geometry])       # (n, 2, 2)
    n_steps = np.maximum(np.ceil(segments.length_m / step_m).astype("int64"), 1)
    share = (segments.duration_s / n_steps).to_numpy()

    reps = np.repeat(np.arange(len(segments)), n_steps)
    within = np.concatenate([(np.arange(k) + 0.5) / k for k in n_steps])
    p0, p1 = coords[reps, 0, :], coords[reps, 1, :]
    pts = p0 + (p1 - p0) * within[:, None]

    inv = ~transform
    col, row = inv * (pts[:, 0], pts[:, 1])
    r, c = np.floor(row).astype("int64"), np.floor(col).astype("int64")
    ok = (r >= 0) & (r < shape[0]) & (c >= 0) & (c < shape[1])
    flat = np.bincount((r[ok] * shape[1] + c[ok]), weights=share[reps][ok],
                       minlength=shape[0] * shape[1])
    return flat.reshape(shape)

Tile the domain. Accumulation is a sum, so tiles are independent and their partial grids add. Split by geographic tile, process in parallel, and sum — which also bounds memory, because no worker holds more than one tile’s segments.

Partition the input by time, not by vessel. Segments are built from consecutive fixes of one vessel, so a split by MMSI keeps every segment intact while a split by time cuts tracks at the boundary and loses one segment per vessel per boundary. For a month-long aggregation that loss is negligible; for a daily product it is a systematic under-count at midnight, which shows up as a faint daily periodicity in any time-series built from the grids.

The conservation test in the previous section is what makes all three safe to adopt: if a refactor drops or double-counts segments, the vessel-hours total stops matching and the test fails before the product ships.

Partition by Vessel, Not by Time Two partitioning schemes applied to the same set of vessel tracks. In the by-vessel scheme, each worker receives all the fixes for a subset of vessels, so every consecutive pair remains within one worker and no segment is lost. In the by-time scheme, each worker receives all vessels for one time window, so the pair of fixes straddling each boundary is split across two workers and the segment between them is never built; the diagram marks one lost segment per vessel per boundary. A chart beneath plots daily vessel-hours from grids built the second way, showing a small but regular dip at each day boundary that reads as a diurnal traffic pattern and is an artefact of the partitioning. by vessel — nothing lost each worker holds whole tracks by time — one segment lost per boundary the straddling pair is split in two daily vessel-hours from grids built the second way a small regular dip at every day boundary — read as a diurnal pattern, produced by the partitioning The artefact is a fraction of a per cent and perfectly periodic, which is exactly what makes it convincing.

Edge cases and gotchas

  • A segment longer than the grid. A vessel that reports either side of a fifteen-minute gap at 20 knots produces a 9 km segment. Densification handles it correctly and slowly; where such segments are common the gap threshold is too generous for the chosen cell size.
  • Densification in a loop is the slow part. For production volumes, replace the Python loop with a vectorised interpolation over all segments at once, or accept the cost and run it per tile. The version above is written for clarity; profile before scaling it.
  • min_coverage is a publication decision, not a technical one. Masking below 0.3 hides the noisiest cells and also hides the fact that the survey extends there. Publishing the coverage band means the consumer can choose their own threshold, which is why both bands go in the file.

Up: Vessel Density & Traffic Raster Generation