Origin-Destination Flow Matrices

An origin-destination flow matrix is a square, zone-by-zone tally of how many trips move from each origin area to each destination area over a defined time window — the workhorse structure that turns millions of raw movement records into a compact, analysable model of where people and vehicles actually go.

Prerequisites Checklist

Before building OD matrices, confirm your inputs and environment are ready. This topic sits within the Spatial Indexing & Density Aggregation section and consumes the output of upstream trajectory processing.

• Python ≥ 3.10 with pandas ≥ 2.0 for tabular aggregation and numpy ≥ 1.24 for array work.
• scipy ≥ 1.11 for scipy.sparse matrix construction — the OD matrix must never be materialised densely at fine zoning.
• geopandas ≥ 0.14 and shapely ≥ 2.0 for polygon-based zone joins.
• h3 ≥ 4.0 (the H3 v4 Python bindings) if you assign zones via the discrete global grid rather than administrative polygons.
• pyproj ≥ 3.6 with always_xy=True for any endpoint distance or centroid work in a metric CRS.
• Input schema: a trip table where each row has origin_lon, origin_lat, dest_lon, dest_lat (float64), start_ts, end_ts (UTC datetime64), and weight (float, defaulting to 1.0). These trips come from upstream segmentation and stop detection — do not build OD matrices directly on raw pings.
• Projection context: endpoints arrive as EPSG:4326. Any distance, centroid, or buffering computation must first project to a metric CRS; zone lookups against a pre-indexed grid may stay in geographic coordinates because H3 indexes lat/lon directly.

Trip and Zone Error Taxonomy

The quality of an OD matrix is determined almost entirely by how you define a trip and where you draw zone boundaries. The failures below are structural, not stochastic — each one silently biases the matrix rather than throwing an error, so each needs an explicit guard.

Deterministic Pipeline Overview

The five-stage sequence below is the standard production order. Each stage preserves trip lineage (keep a trip_id), records the zoning scheme and time window in metadata, and is fully vectorized — no Python-level row iteration.

1. Define trips — reduce each segmented trajectory to one trip record with origin/destination coordinates, start_ts, end_ts, and weight. A trip is the unit of flow; get its definition right before anything else.
2. Assign zones — map each origin and destination point to a zone identifier, either by H3 indexing at a chosen resolution or by a polygon spatial join in a metric CRS. Trips with an unresolvable endpoint go to an explicit unknown bucket.
3. Aggregate to a sparse matrix — group trips by (origin_zone, dest_zone), sum the weights, and emit both a scipy.sparse.coo_matrix for computation and a labelled pandas pivot for inspection.
4. Normalise & threshold — apply row normalisation (outflow shares), column normalisation (inflow shares), or global normalisation as the question requires, then suppress or aggregate cells below the minimum trip count.
5. Slice temporally & export — partition the trip table by time window before aggregation to produce comparable hourly, daily, or seasonal matrices, then export the sparse matrix, the zone index, and the metadata together.

Implementation Walkthrough

The function below implements stages 1 through 4 on a table of trips. It assumes zones are already assigned as integer-encoded origin_zone / dest_zone labels (the companion deep-dive covers zone assignment from raw endpoints end to end). Any distance work here happens in a metric projected CRS, never on raw EPSG:4326 degrees.

from __future__ import annotations

import logging
from dataclasses import dataclass

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

logger = logging.getLogger(__name__)


@dataclass
class ODResult:
    """Container for an origin-destination aggregation result."""
    matrix: sparse.csr_matrix          # zone x zone trip counts (sparse)
    zone_index: pd.Index               # row/col position -> zone label
    pivot: pd.DataFrame                # labelled dense view (small slices only)
    diagonal: pd.Series                # self-loop counts per zone
    suppressed_cells: int              # cells removed by the min-trip threshold


def build_od_matrix(
    trips: pd.DataFrame,
    origin_col: str = "origin_zone",
    dest_col: str = "dest_zone",
    weight_col: str = "weight",
    min_trips: int = 5,
    keep_self_loops: bool = True,
    normalise: str | None = None,       # None | "row" | "col" | "global"
) -> ODResult:
    """
    Aggregate a trip table into a sparse origin-destination flow matrix.

    Parameters
    ----------
    trips : DataFrame with integer/string zone labels in origin_col and
            dest_col, plus an optional numeric weight_col (defaults to 1.0
            per trip). Rows with a missing zone in either endpoint are
            dropped and counted; assign them to an explicit 'unknown' zone
            upstream if you need to retain them.
    min_trips : Cells with a summed weight below this are suppressed
                (set to zero) for statistical reliability and privacy.
    keep_self_loops : If False, same-zone trips are removed from the matrix
                      body; their counts are always returned in `diagonal`.
    normalise : Row-normalise (outflow shares), column-normalise (inflow
                shares), global-normalise, or leave as raw counts (None).

    Returns
    -------
    ODResult with a sparse CSR matrix, the zone index that maps matrix
    positions back to zone labels, a labelled pivot for inspection, the
    self-loop diagonal, and the count of suppressed cells.

    Raises
    ------
    ValueError if required columns are missing.
    """
    required = {origin_col, dest_col}
    missing = required - set(trips.columns)
    if missing:
        raise ValueError(f"Trip table missing columns: {missing}")

    if trips.empty:
        logger.warning("build_od_matrix received an empty trip table.")
        empty = sparse.csr_matrix((0, 0))
        return ODResult(empty, pd.Index([]), pd.DataFrame(), pd.Series(dtype=float), 0)

    df = trips.copy()
    if weight_col not in df.columns:
        df[weight_col] = 1.0

    # Drop trips with an unresolved endpoint; never coerce NaN to a real zone.
    n_before = len(df)
    df = df.dropna(subset=[origin_col, dest_col])
    n_dropped = n_before - len(df)
    if n_dropped:
        logger.info("Dropped %d trips with a missing origin/dest zone.", n_dropped)

    if df.empty:
        empty = sparse.csr_matrix((0, 0))
        return ODResult(empty, pd.Index([]), pd.DataFrame(), pd.Series(dtype=float), 0)

    # ── Stage 3: aggregate by (origin, dest) ─────────────────────────────────
    grouped = (
        df.groupby([origin_col, dest_col], observed=True)[weight_col]
        .sum()
        .reset_index()
    )

    # Split off the self-loop diagonal before it can dominate the body.
    is_self = grouped[origin_col] == grouped[dest_col]
    diagonal = (
        grouped.loc[is_self]
        .set_index(origin_col)[weight_col]
        .rename("self_loops")
        .sort_index()
    )
    if not keep_self_loops:
        grouped = grouped.loc[~is_self].copy()

    # ── Stage 4a: threshold small cells ──────────────────────────────────────
    below = grouped[weight_col] < min_trips
    suppressed = int(below.sum())
    if suppressed:
        logger.info("Suppressing %d cells below min_trips=%d.", suppressed, min_trips)
    grouped = grouped.loc[~below].copy()

    # Build a stable, sorted zone index spanning every zone that still appears.
    zones = pd.Index(
        sorted(set(grouped[origin_col]) | set(grouped[dest_col])),
        name="zone",
    )
    pos = {z: i for i, z in enumerate(zones)}
    rows = grouped[origin_col].map(pos).to_numpy()
    cols = grouped[dest_col].map(pos).to_numpy()
    data = grouped[weight_col].to_numpy(dtype=float)

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

    # ── Stage 4b: optional normalisation ─────────────────────────────────────
    if normalise == "row":
        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
    elif normalise == "col":
        col_sums = np.asarray(mat.sum(axis=0)).ravel()
        inv = np.divide(1.0, col_sums, out=np.zeros_like(col_sums), where=col_sums > 0)
        mat = mat @ sparse.diags(inv)
    elif normalise == "global":
        total = mat.sum()
        if total > 0:
            mat = mat.multiply(1.0 / total).tocsr()

    # Labelled pivot for inspection — safe to densify only because thresholding
    # and (optionally) self-loop removal keep the surviving zone set small.
    pivot = pd.DataFrame(mat.toarray(), index=zones, columns=zones)

    return ODResult(
        matrix=mat,
        zone_index=zones,
        pivot=pivot,
        diagonal=diagonal,
        suppressed_cells=suppressed,
    )

The ODResult bundles everything a downstream consumer needs: the sparse matrix for computation, the zone_index to map matrix positions back to real zone labels, the labelled pivot for eyeballing, the self-loop diagonal kept out of the body, and an audit count of suppressed cells. Because the pivot is only densified after thresholding, it stays small enough to inspect; at production scale you carry the sparse matrix and the index, not the dense array.

Geometric and Mathematical Grounding

The OD matrix definition. Let the study area be partitioned into a set of zones Z = {z_1, ..., z_N}. The origin-destination matrix T is the N × N array whose entry T[i, j] is the total weight of trips whose origin lies in z_i and whose destination lies in z_j:

T[i, j] = Σ_k w_k · 1[origin(k) ∈ z_i] · 1[dest(k) ∈ z_j]

where the sum runs over all trips k, w_k is the trip weight (usually 1), and 1[·] is the indicator function. The row sum Σ_j T[i, j] is the total outflow (trip production) from zone i; the column sum Σ_i T[i, j] is the total inflow (trip attraction) to zone j; the diagonal T[i, i] holds intra-zone self-loops.

Row and column normalisation. Raw counts answer “how many trips?”; normalised matrices answer “what share?”. Row normalisation divides each row by its outflow, P[i, j] = T[i, j] / Σ_j T[i, j], giving the conditional probability that a trip leaving z_i arrives at z_j — the rows then sum to 1 and form a transition matrix. Column normalisation divides each column by its inflow and answers the reverse question: given an arrival in z_j, where did it come from. Global normalisation divides by the grand total and yields the joint distribution over all zone pairs. Choose the normalisation to match the question, and always keep the raw counts alongside the shares, because a 100 percent share backed by two trips is not the same finding as one backed by two thousand.

Gravity-model intuition. OD flows are not random. The classic gravity model predicts that the flow between two zones scales with the product of their “masses” (population, jobs, activity) and decays with the cost of travel between them: T[i, j] ≈ G · (M_i · M_j) / f(c_ij), where M_i and M_j are the production and attraction masses, c_ij is the inter-zone travel cost (distance or time computed in a metric CRS), and f is a deterrence function such as c_ij^β or exp(β · c_ij). You do not need to fit a gravity model to build an OD matrix, but the intuition is diagnostic: a measured matrix that shows strong flows between distant, low-mass zones usually signals a zoning artifact or a segmentation bug rather than genuine long-distance demand.

Calibration & Parameter Tuning

Zone size and the minimum trip threshold are the two dials that most change the output. The table pairs a zoning scheme with the use case it suits, and gives starting thresholds by mode. Distances that inform these choices are always computed in a metric CRS.

Two rules of thumb keep the matrix meaningful. First, keep the zone diameter well below the median trip length: when zones approach trip length, most trips collapse onto the diagonal as self-loops. Second, set min_trips from both reliability and privacy: 5 is a common floor for person-based data, higher when the matrix is published externally. Calibrate against a held-out subset rather than theory, and re-run at a second resolution to gauge MAUP sensitivity before trusting any single matrix.

Integration & Compatibility

OD matrices sit at a junction between trajectory analysis and spatial aggregation, so they connect to several neighbouring topics:

• Discrete Global Grid Systems: the cleanest source of zone identifiers. Assigning origins and destinations to H3 cells gives you a globally consistent, hierarchical zoning scheme that sidesteps ad hoc polygon boundaries and makes MAUP sensitivity testing a one-line resolution change. See aggregating GPS pings into H3 hexagon grids for the indexing mechanics.
• Trajectory Segmentation: the definition of a trip. Where one trajectory is split into trips determines every downstream OD entry — a too-eager splitter manufactures short phantom trips, a too-lazy one merges a home-work-home tour into a single self-loop.
• Stay-Point Detection Algorithms: the natural source of trip endpoints. The stops a device dwells at are the origins and destinations between which trips run, so stay-point output feeds directly into stage 1 of this pipeline.
• Coordinate Reference System Mapping: any endpoint distance, centroid, or polygon buffering must happen in a metric CRS. Reuse the same projected CRS across your whole pipeline so travel-cost figures in a gravity check line up with the zoning.
• Kernel Density Surfaces: the continuous counterpart. Where an OD matrix models discrete zone-to-zone flows, a density surface models where movement concentrates; the two are complementary views of the same trips.

FAQ

What zone size should I use for an OD flow matrix?

Match the zone size to the trip length distribution and the decision you are supporting. For city-scale mobility, H3 resolution 7 or 8 balances detail against sparsity. For regional or inter-city flows, administrative polygons or H3 resolution 5 to 6 are more stable. Zones much smaller than the median trip length inflate the matrix and produce mostly empty cells; zones much larger than trip length collapse everything onto the diagonal as self-loops.

How do I handle trips whose origin and destination fall in the same zone?

Same-zone trips populate the matrix diagonal as self-loops. They are legitimate intra-zone movement, but they dominate the matrix when zones are large relative to trips and often mask the inter-zone structure analysts care about. Keep them in a separate diagonal vector rather than deleting them, and decide per use case whether to visualise them. For flow maps, self-loops are usually rendered as a node symbol sized by intra-zone volume rather than as an arc.

Why is my OD matrix almost entirely zeros?

An N-zone study area produces an N-by-N matrix with N² cells, but real mobility only ever fills a tiny fraction of them, because most zone pairs never exchange a trip. This is expected: never materialise the full dense matrix. Store the aggregate as a scipy.sparse matrix keyed on the (origin, destination) pairs that actually occur, and only densify a small slice for display. Sparsity above 99 percent is normal at fine zoning resolutions.

How does the Modifiable Areal Unit Problem affect OD matrices?

The Modifiable Areal Unit Problem means the flows you measure depend on where you draw zone boundaries. Shifting or resizing zones changes which trips count as inter-zone versus intra-zone and can even reverse apparent flow directions. Mitigate it by choosing a zoning scheme tied to the analytical question, testing sensitivity across at least two resolutions, and preferring a consistent global grid such as H3 over ad hoc polygons when comparability matters.

Should I suppress small cells before publishing an OD matrix?

Yes, when the data derives from individuals. Cells backed by only one or two trips are both statistically unreliable and a privacy re-identification risk. Apply a minimum trip threshold (commonly 5 or more) and either drop, aggregate to a coarser zone, or add calibrated noise to cells below it. Record the suppression rule as metadata so downstream consumers know the matrix is thresholded and do not treat missing cells as true zeros.

Related

• Building Origin-Destination Matrices from Trajectory Data
• Discrete Global Grid Systems
• Trajectory Segmentation
• Stay-Point Detection Algorithms
• Kernel Density Surfaces

Back to Spatial Indexing & Density Aggregation