Building origin-destination matrices from trajectory data

This guide walks the full path from a set of segmented trajectories to a validated, exportable origin-destination matrix: derive trip endpoints → assign zones → build a sparse matrix → normalise and export. The result is a compact scipy.sparse structure that records how many trips ran between every pair of zones in a time window, keyed on a zone index you can map back to real geography. Every distance calculation happens in a metric projected CRS, the matrix is never densified at scale, and each stage is guarded against the edge cases — single-point segments, missing zones, empty frames — that otherwise corrupt the totals silently.

Why this starts with well-defined trips

An OD matrix is only ever as good as its trip definition, and a trip is defined upstream of this page. The endpoints you aggregate come from splitting each device’s path into discrete journeys, which is the job of trajectory segmentation: a segmenter that fires too eagerly manufactures short phantom trips, while one that merges a whole day into one path collapses everything onto the matrix diagonal. The conceptual model behind the matrix itself — row and column sums, self-loops, normalisation, and the Modifiable Areal Unit Problem — is covered in the parent origin-destination flow matrices reference and in the broader Spatial Indexing & Density Aggregation section. Here we focus on the mechanics of turning trajectories into the matrix.

Core pipeline

The four stages below run in order: you cannot assign a zone before you have an endpoint, and you cannot aggregate before every endpoint has a zone.

OD Matrix Construction Pipeline Four sequential stages: Derive Trip Endpoints, Assign Zones, Build Sparse Matrix, and Normalise and Export, connected by arrows. 1. Derive Trip Endpoints 2. Assign Zones 3. Build Sparse Matrix 4. Normalise & Export
  1. Derive trip endpoints. Reduce each segment to its first and last valid fix — the origin and destination — carrying timestamps and a weight of 1. Drop segments with fewer than two distinct points.
  2. Assign zones. Map each endpoint to a zone. With H3 this is a direct lat/lon-to-cell lookup; with polygons it is a spatial join performed in a metric CRS. Endpoints outside every zone go to an explicit unknown bucket.
  3. Build the sparse matrix. Group trips by (origin_zone, dest_zone), sum weights, and construct a scipy.sparse.coo_matrix plus a labelled pandas pivot for inspection.
  4. Normalise and export. Optionally convert counts to outflow shares by row-normalising, then write the matrix, the zone index, and the run metadata together.

Production-ready Python implementation

The function below runs stages 1 through 3 and optional normalisation. It supports both zoning strategies through a single zone_fn callback, so the aggregation logic never needs to know whether zones came from H3 or a polygon join. The H3 path indexes lat/lon directly; the polygon path is shown in the helper and does its join in a metric CRS. The endpoint derivation is fully vectorized — no iterrows.

PYTHON 
from __future__ import annotations

from typing import Callable

import numpy as np
import pandas as pd
from scipy import sparse

# Optional imports guarded so the module loads without H3/GeoPandas installed.
try:
    import h3
except ImportError:  # pragma: no cover
    h3 = None


def derive_trip_endpoints(
    points: pd.DataFrame,
    segment_col: str = "segment_id",
    time_col: str = "ts",
    lat_col: str = "lat",
    lon_col: str = "lon",
) -> pd.DataFrame:
    """
    Collapse a table of segmented trajectory points into one trip per segment.

    Each trip takes the first fix of its segment as the origin and the last
    fix as the destination. Segments with fewer than two distinct points are
    stationary dwells, not trips, and are dropped.

    Returns a frame with origin_lat/lon, dest_lat/lon, start_ts, end_ts, weight.
    """
    required = {segment_col, time_col, lat_col, lon_col}
    missing = required - set(points.columns)
    if missing:
        raise ValueError(f"Points table missing columns: {missing}")
    if points.empty:
        return pd.DataFrame(
            columns=["origin_lat", "origin_lon", "dest_lat", "dest_lon",
                     "start_ts", "end_ts", "weight"]
        )

    df = points.sort_values([segment_col, time_col])
    grp = df.groupby(segment_col, observed=True)

    # Vectorized first/last fix per segment via idxmin/idxmax on time.
    counts = grp[time_col].transform("size")
    df = df.loc[counts >= 2]
    if df.empty:
        return pd.DataFrame(
            columns=["origin_lat", "origin_lon", "dest_lat", "dest_lon",
                     "start_ts", "end_ts", "weight"]
        )

    grp = df.groupby(segment_col, observed=True)
    first = grp.first()
    last = grp.last()

    trips = pd.DataFrame({
        "origin_lat": first[lat_col].to_numpy(),
        "origin_lon": first[lon_col].to_numpy(),
        "dest_lat": last[lat_col].to_numpy(),
        "dest_lon": last[lon_col].to_numpy(),
        "start_ts": first[time_col].to_numpy(),
        "end_ts": last[time_col].to_numpy(),
        "weight": 1.0,
    })
    return trips


def assign_h3_zones(
    trips: pd.DataFrame, resolution: int = 8
) -> pd.DataFrame:
    """Assign origin/dest H3 cells. H3 indexes lat/lon directly (no CRS needed)."""
    if h3 is None:
        raise ImportError("The 'h3' package (v4) is required for H3 zone assignment.")
    out = trips.copy()
    # h3.latlng_to_cell is scalar; vectorize with a list comprehension over arrays.
    out["origin_zone"] = [
        h3.latlng_to_cell(la, lo, resolution)
        for la, lo in zip(out["origin_lat"], out["origin_lon"])
    ]
    out["dest_zone"] = [
        h3.latlng_to_cell(la, lo, resolution)
        for la, lo in zip(out["dest_lat"], out["dest_lon"])
    ]
    return out


def assign_polygon_zones(
    trips: pd.DataFrame, zones_gdf, zone_id_col: str = "zone_id",
    metric_epsg: int = 32633,   # UTM 33N — pick the zone covering your data
):
    """
    Assign zones by spatial join. The join runs in a metric projected CRS
    (never raw EPSG:4326) so any nearest-neighbour tolerance is in metres.
    Endpoints outside every polygon receive the sentinel zone 'unknown'.
    """
    import geopandas as gpd

    zones_m = zones_gdf.to_crs(epsg=metric_epsg)
    out = trips.copy()
    for role in ("origin", "dest"):
        pts = gpd.GeoDataFrame(
            out,
            geometry=gpd.points_from_xy(out[f"{role}_lon"], out[f"{role}_lat"]),
            crs="EPSG:4326",
        ).to_crs(epsg=metric_epsg)   # project before any spatial predicate
        joined = gpd.sjoin(
            pts, zones_m[[zone_id_col, "geometry"]],
            how="left", predicate="within",
        )
        joined = joined[~joined.index.duplicated(keep="first")]
        out[f"{role}_zone"] = joined[zone_id_col].fillna("unknown").to_numpy()
    return out


def build_od_from_trips(
    trips: pd.DataFrame,
    zone_fn: Callable[[pd.DataFrame], pd.DataFrame],
    min_trips: int = 5,
    normalise_rows: bool = False,
) -> tuple[sparse.csr_matrix, pd.Index]:
    """
    End-to-end: assign zones via zone_fn, then aggregate to a sparse OD matrix.

    Returns (matrix, zone_index). zone_index maps each matrix row/column
    position back to its zone label. Cells below min_trips are suppressed.
    """
    if trips.empty:
        return sparse.csr_matrix((0, 0)), pd.Index([], name="zone")

    zoned = zone_fn(trips)
    zoned = zoned.dropna(subset=["origin_zone", "dest_zone"])

    grouped = (
        zoned.groupby(["origin_zone", "dest_zone"], observed=True)["weight"]
        .sum()
        .reset_index()
    )
    grouped = grouped.loc[grouped["weight"] >= min_trips]
    if grouped.empty:
        return sparse.csr_matrix((0, 0)), pd.Index([], name="zone")

    zones = pd.Index(
        sorted(set(grouped["origin_zone"]) | set(grouped["dest_zone"])),
        name="zone",
    )
    pos = {z: i for i, z in enumerate(zones)}
    rows = grouped["origin_zone"].map(pos).to_numpy()
    cols = grouped["dest_zone"].map(pos).to_numpy()
    data = grouped["weight"].to_numpy(dtype=float)

    n = len(zones)
    mat = sparse.coo_matrix((data, (rows, cols)), shape=(n, n)).tocsr()

    if normalise_rows:
        row_sums = np.asarray(mat.sum(axis=1)).ravel()
        inv = np.divide(1.0, row_sums, out=np.zeros_like(row_sums), where=row_sums > 0)
        mat = sparse.diags(inv) @ mat

    return mat, zones

To run the whole thing you compose the pieces: endpoints = derive_trip_endpoints(points), then mat, zones = build_od_from_trips(endpoints, lambda t: assign_h3_zones(t, resolution=8)). Swap the lambda for assign_polygon_zones when your zones are administrative. Export is a scipy.sparse.save_npz(path, mat) for the matrix plus a zones.to_series().to_json() for the index — always ship the two together, because the matrix positions are meaningless without the label map.

Validation block

Run these checks before the matrix leaves your pipeline. The central invariant is conservation: no trip may vanish or be counted twice between the trip table and the matrix.

PYTHON 
def validate_od_matrix(
    mat: sparse.csr_matrix,
    zones: pd.Index,
    zoned_trips: pd.DataFrame,
    min_trips: int = 5,
) -> None:
    """
    Reconcile a sparse OD matrix against the trip table it came from.
    Raises AssertionError on any inconsistency; prints a summary on success.
    """
    assert mat.shape[0] == mat.shape[1] == len(zones), (
        "Matrix is not square or zone_index length does not match its dimensions"
    )

    valid = zoned_trips.dropna(subset=["origin_zone", "dest_zone"])

    # Per-pair counts that SHOULD survive the min_trips threshold.
    pair_counts = valid.groupby(
        ["origin_zone", "dest_zone"], observed=True
    )["weight"].sum()
    surviving = pair_counts[pair_counts >= min_trips].sum()

    # Conservation: matrix total equals the sum of un-suppressed pair weights.
    matrix_total = float(mat.sum())
    assert np.isclose(matrix_total, surviving), (
        f"Conservation failed: matrix sums to {matrix_total} "
        f"but surviving trips sum to {surviving}"
    )

    # Row sums must equal each zone's un-suppressed outflow (raw counts only).
    pos = {z: i for i, z in enumerate(zones)}
    row_sums = np.asarray(mat.sum(axis=1)).ravel()
    outflow = (
        pair_counts[pair_counts >= min_trips]
        .groupby(level=0).sum()
    )
    for zone, expected in outflow.items():
        if zone in pos:
            assert np.isclose(row_sums[pos[zone]], expected), (
                f"Row sum mismatch for zone {zone}: "
                f"{row_sums[pos[zone]]} != {expected}"
            )

    print(
        f"Validation passed. {len(zones)} zones, "
        f"{mat.nnz} non-zero flows, {matrix_total:.0f} trips retained, "
        f"density {mat.nnz / (len(zones) ** 2):.4%}."
    )

The conservation check catches the most common corruption modes at once: a dropped trip lowers the surviving total below the matrix total, a double-count raises it, and a NaN zone coerced to a real label shifts weight between cells. Run this on the raw-count matrix before any normalisation, because row-normalised rows sum to 1 and no longer reconcile against trip counts.

Common mistakes and gotchas

  • Assigning zones in EPSG:4326 for the polygon path. A within predicate is topological and works in any CRS, but the moment you add a distance tolerance, a nearest-zone fallback, or a buffer, degree units silently break it. Project to a metric CRS before any spatial predicate that involves distance, exactly as the polygon helper above does.

  • Snapping unmatched endpoints to the nearest polygon. A trip that ends outside every zone is telling you about a coverage gap. Snapping it to the nearest zone invents a flow that never happened; dropping it biases totals downward. Send it to an explicit unknown zone so the gap stays visible in validation and downstream reporting.

  • Densifying the matrix to inspect it. Calling .toarray() on a national-scale matrix allocates an N² dense array and can exhaust memory instantly. Keep the scipy.sparse form, and only densify a small named slice — a handful of zones of interest — when you actually need to read cells by eye.

  • Treating a suppressed cell as a true zero. After small-cell suppression, a missing (i, j) entry means “below the threshold”, not “no trips”. Downstream code that interprets absence as zero will underestimate low-volume flows. Carry the suppression rule in the export metadata and, where it matters, keep the pre-suppression totals.

  • Mixing time windows. Aggregating a full week of trips into one matrix averages the morning peak, the evening peak, and the quiet small hours into a composite that describes no real moment. Slice the trip table by window — hour of day, day type, or season — before aggregation, and align that slicing with your aggregating GPS pings into H3 hexagon grids resolution so zones and windows stay comparable across runs.

Back to Origin-Destination Flow Matrices