Map-Matching to Road Networks

Map-matching is the process of snapping a sequence of noisy GPS observations onto the specific road-network edges a vehicle actually travelled, converting a scatter of imprecise coordinates into a connected, routable path along known geometry.

Prerequisites Checklist

Before running an HMM map-matcher, confirm the input trajectory and the road graph meet these requirements. This topic sits inside the Spatial Indexing & Density Aggregation section and assumes you have already produced a clean, projected trace.

• Python ≥ 3.10 with pandas ≥ 2.0, numpy ≥ 1.24, geopandas ≥ 0.14, and shapely ≥ 2.0 for vectorized geometry operations.
• osmnx ≥ 1.6 and networkx ≥ 3.0 to download and traverse the road-network graph, or a pre-built graph in an equivalent structure.
• pyproj ≥ 3.6 for datum transformations — imported with always_xy=True to avoid axis-order bugs.
• A cleaned input trajectory: each observation must have x_m, y_m (metric projected coordinates), ts (UTC datetime), and a segment_id. Drift, phantom velocity spikes, and duplicate timestamps must already be resolved by the GPS precision and error handling stage — map-matching amplifies residual noise rather than absorbing it.
• A metric projected CRS shared by both the trace and the road graph. Every distance in the emission and transition terms is metres; raw EPSG:4326 degrees are never used for scoring.
• Calibrated sigma and beta for the transport mode and sampling rate of your data (see the calibration table below).

Error and Problem Taxonomy

Map-matching failures are almost always failures of disambiguation: several candidate edges are geometrically plausible for a single observation, and a naive matcher picks the wrong one. Classifying the failure mode determines whether the fix belongs in candidate generation, emission scoring, or transition scoring. Getting this wrong corrupts everything downstream, from origin-destination flow matrices to turn-level congestion metrics.

Deterministic Pipeline Overview

The six stages below are the standard production order. Each stage is fully vectorized where possible, keeps candidate metadata for auditability, and never computes a distance outside the metric CRS.

1. Ingest the cleaned trace. Consume a drift-corrected trajectory already projected to a metric CRS, with monotonic UTC timestamps and stable segment_id boundaries from the upstream cleaning pipeline.
2. Build the road graph. Download the OSM drive network for the trajectory’s buffered bounding box, project every edge geometry to the same metric CRS, and build a spatial index (STR-tree) over the edges.
3. Generate candidate edges per point. For each observation, query the spatial index for edges within a search radius (typically three to four times sigma), project the point onto each candidate line, and record the perpendicular distance and the snapped position.
4. Score emission and transition probabilities. Emission scores each candidate by how close the observation lies to its edge (a zero-mean Gaussian on perpendicular distance). Transition scores each pair of consecutive candidates by comparing the straight-line inter-observation distance to the on-network shortest-path distance between the snapped points.
5. Decode with Viterbi. Run the Viterbi algorithm over the candidate lattice in log space to recover the single most probable sequence of edges — the joint optimum, not a chain of greedy local choices.
6. Emit the matched path. Return matched edge IDs, the snapped point geometry, and per-point match confidence, splitting the output at any lattice break caused by an unmatched observation.

Implementation Walkthrough

The core of a map-matcher is candidate scoring. The function below extracts candidate road segments for one observation and computes the Gaussian emission probability from the perpendicular distance. All geometry is handled in a metric projected CRS (EPSG:32633 as default) — never raw EPSG:4326 — so the perpendicular distances feeding the Gaussian are true metres.

from __future__ import annotations

import math
from dataclasses import dataclass

import geopandas as gpd
import numpy as np
from shapely.geometry import Point
from shapely.strtree import STRtree


@dataclass(frozen=True)
class Candidate:
    """A single candidate edge match for one GPS observation."""
    edge_id: int          # index into the edges GeoDataFrame
    snapped_x: float      # projected snap point (metres, metric CRS)
    snapped_y: float
    perp_dist_m: float    # perpendicular distance observation -> edge (metres)
    emission_logp: float  # log emission probability


def extract_candidates(
    obs_x: float,
    obs_y: float,
    edges: gpd.GeoDataFrame,
    edge_index: STRtree,
    sigma_m: float = 8.0,
    search_radius_m: float = 30.0,
    max_candidates: int = 8,
) -> list[Candidate]:
    """
    Return scored candidate edges for a single GPS observation.

    Emission probability models the perpendicular GPS error as a
    zero-mean Gaussian with standard deviation sigma_m:

        p(obs | edge) = (1 / (sqrt(2*pi) * sigma)) * exp(-0.5 * (d/sigma)**2)

    Parameters
    ----------
    obs_x, obs_y : Observation coordinates in a METRIC projected CRS
                   (e.g. UTM). NEVER pass EPSG:4326 degrees — the
                   perpendicular distance must be in metres.
    edges : GeoDataFrame of LineString road segments in the SAME
            metric CRS as the observation.
    edge_index : STRtree built over edges.geometry for fast radius query.
    sigma_m : GPS perpendicular-error standard deviation (metres).
    search_radius_m : Maximum snap distance; edges beyond this are ignored.
    max_candidates : Keep only the closest N candidates per observation.

    Returns
    -------
    list[Candidate], sorted by ascending perpendicular distance. Empty
    if no edge falls within search_radius_m (an off-road / gap point).
    """
    if edges.empty:
        raise ValueError("edges GeoDataFrame is empty; build the road graph first.")
    if math.isnan(obs_x) or math.isnan(obs_y):
        return []  # masked / interpolated point with no valid position

    obs = Point(obs_x, obs_y)
    norm = 1.0 / (math.sqrt(2.0 * math.pi) * sigma_m)

    # STRtree.query returns integer positions of candidate geometries whose
    # bounding boxes intersect the buffer — a coarse metric-space filter.
    hit_positions = edge_index.query(obs.buffer(search_radius_m))

    candidates: list[Candidate] = []
    for pos in np.atleast_1d(hit_positions):
        geom = edges.geometry.iloc[int(pos)]
        # project() + interpolate() give the perpendicular foot of the
        # observation on the edge, all in metric units.
        proj_dist = geom.project(obs)
        snap = geom.interpolate(proj_dist)
        perp = obs.distance(geom)  # perpendicular distance in metres
        if perp > search_radius_m:
            continue
        emission_logp = math.log(norm) - 0.5 * (perp / sigma_m) ** 2
        candidates.append(
            Candidate(
                edge_id=int(edges.index[int(pos)]),
                snapped_x=snap.x,
                snapped_y=snap.y,
                perp_dist_m=perp,
                emission_logp=emission_logp,
            )
        )

    candidates.sort(key=lambda c: c.perp_dist_m)
    return candidates[:max_candidates]

The companion transition score compares how far the vehicle appears to have moved (straight-line distance between two observations) against how far it would have to drive on the network (shortest-path route distance between the two snapped points). When these agree the transition is likely; when the route is far longer than the straight line, the candidate pair is implausible.

import networkx as nx


def transition_logp(
    great_circle_m: float,
    route_dist_m: float,
    beta_m: float = 5.0,
) -> float:
    """
    Log transition probability between two consecutive candidates.

    Follows the Newson & Krumm (2009) formulation: the difference
    between the straight-line inter-observation distance and the
    on-network route distance is exponentially distributed with
    mean beta_m.

        p(edge_j | edge_i) = (1 / beta) * exp(-|gc - route| / beta)

    Both distances MUST be in metres (metric CRS). A route distance
    much larger than the straight-line distance signals an implausible
    jump (wrong turn, parallel road, or graph gap).
    """
    if not math.isfinite(route_dist_m):
        return -math.inf  # no route -> disconnected candidates
    delta = abs(great_circle_m - route_dist_m)
    return math.log(1.0 / beta_m) - delta / beta_m


def route_distance_m(
    graph: nx.MultiDiGraph,
    node_from: int,
    node_to: int,
    cutoff_m: float,
) -> float:
    """
    Shortest on-network route length (metres) between two graph nodes,
    bounded by cutoff_m to keep sparse-sampling searches tractable.
    Returns +inf if no path exists within the cutoff (disconnected).
    """
    try:
        return nx.shortest_path_length(
            graph, node_from, node_to, weight="length", method="dijkstra"
        )
    except (nx.NetworkXNoPath, nx.NodeNotFound):
        return math.inf

Finally, the Viterbi decoder walks the candidate lattice and recovers the globally optimal edge sequence. Working in log space avoids floating-point underflow when many small probabilities multiply along a long trace.

def viterbi_decode(
    emission_logp: list[list[float]],
    transition_logp: list[list[list[float]]],
) -> list[int]:
    """
    Decode the most probable candidate sequence over a lattice.

    Parameters
    ----------
    emission_logp : emission_logp[t][i] = log emission of candidate i at step t.
    transition_logp : transition_logp[t][i][j] = log transition from
                      candidate i at step t to candidate j at step t+1.

    Returns
    -------
    list[int] : chosen candidate index at each time step. Empty if the
                lattice has no steps.
    """
    n_steps = len(emission_logp)
    if n_steps == 0:
        return []

    # Initialise with the first step's emissions.
    viterbi = [list(emission_logp[0])]
    backpointer: list[list[int]] = [[-1] * len(emission_logp[0])]

    for t in range(1, n_steps):
        prev = viterbi[t - 1]
        step_scores: list[float] = []
        step_bp: list[int] = []
        for j in range(len(emission_logp[t])):
            best_prev, best_score = -1, -math.inf
            for i in range(len(prev)):
                score = prev[i] + transition_logp[t - 1][i][j]
                if score > best_score:
                    best_score, best_prev = score, i
            step_scores.append(best_score + emission_logp[t][j])
            step_bp.append(best_prev)
        viterbi.append(step_scores)
        backpointer.append(step_bp)

    # Back-trace from the best final state.
    best_last = int(np.argmax(viterbi[-1]))
    path = [best_last]
    for t in range(n_steps - 1, 0, -1):
        best_last = backpointer[t][best_last]
        path.append(best_last)
    path.reverse()
    return path

Geometric and Mathematical Grounding

The hidden Markov formulation. Map-matching casts the road-network edges as hidden states and the GPS observations as emissions of those states. At each time step t the vehicle occupies some true edge s_t, which we cannot observe directly; we observe only a noisy position z_t. The Hidden Markov Model factorises the joint probability of a state sequence and the observations as a product of emission terms p(z_t | s_t) and transition terms p(s_t | s_{t-1}), and map-matching seeks the state sequence that maximises this product.

Emission probability. The GPS error perpendicular to the road is modelled as a zero-mean Gaussian. For a candidate edge at perpendicular distance d_t from observation z_t:

p(z_t | s_t) = (1 / (√(2π)·σ)) · exp(−½ · (d_t / σ)²)

Here σ is the standard deviation of the GPS perpendicular error in metres. Small d_t (the observation lands almost on the edge) gives high emission probability; a candidate 25 m away with σ = 8 is exponentially penalised.

Transition probability. The transition term encodes the physical intuition that the on-network route distance between two consecutive snapped points should closely match the straight-line distance between the raw observations. Following the Newson and Krumm formulation, the absolute difference is exponentially distributed:

p(s_t | s_{t-1}) = (1 / β) · exp(−|‖z_t − z_{t-1}‖ − route(s_{t-1}, s_t)| / β)

where route(·) is the Dijkstra shortest-path length on the graph and β controls how sharply route-vs-straight-line mismatch is penalised. A candidate pair whose route detours far beyond the straight-line gap — the signature of a wrong turn or a jump to a parallel road — receives a vanishing transition probability.

Viterbi decoding. Multiplying emission and transition terms over the whole trace and maximising is exactly the problem the Viterbi algorithm solves in O(T · K²) time for T observations and K candidates per observation. Each lattice cell stores the best log-probability of reaching that candidate and a back-pointer to its predecessor; a single backward pass recovers the optimal edge sequence. Working in log space turns the product into a sum and prevents underflow on long traces.

Calibration & Parameter Tuning

The two dominant parameters are the GPS error standard deviation σ (emission sharpness) and the transition decay β (route-mismatch tolerance). Both depend on transport mode and, critically, on sampling rate — as intervals lengthen the transition term must tolerate larger route-vs-straight-line gaps.

Calibrate σ and β against a small hand-matched ground-truth subset rather than trusting defaults. A practical approach: fix σ from the observed perpendicular residuals of an obviously-correct straight segment, then grid-search β to maximise matched-edge accuracy on the labelled set. For multi-modal fleets, infer mode from speed and acceleration profiling first and run a separate parameter set per mode — a single “unknown” configuration will over-penalise pedestrians or under-constrain vehicles.

Integration & Compatibility

Map-matching sits between trajectory cleaning and network-level aggregation, so it both consumes and feeds several adjacent stages:

• GPS Precision & Error Handling: the matcher consumes the drift-corrected, metric-projected output of this cleaning stage. Feeding raw traces directly into an HMM inflates perpendicular residuals and forces an artificially large σ, blurring the emission term.
• Handling GPS Drift in Raw Trajectory Logs: the specific step-by-step wiring of a cleaned trace into a matcher is covered in the companion deep-dive, Feeding Cleaned GPS Output into an HMM Map-Matcher.
• Coordinate Reference System Mapping: the trace and the road graph must share one metric CRS. Use the same UTM zone for both to avoid re-projection drift between the snap geometry and the network edges.
• Origin-Destination Flow Matrices: matched edge sequences are the ideal input for building network-aware OD matrices, because trips are now expressed as routable paths rather than point clouds.
• Sampling Rate Optimization: knowing your effective sampling interval drives the β choice; excessively downsampled traces need gap-filling before matching stays reliable.

FAQ

Why use a Hidden Markov Model instead of nearest-edge snapping?

Nearest-edge snapping decides each point independently, so at a junction or between parallel roads it hops between the wrong edges and produces a physically impossible path. An HMM couples an emission term (how close the observation is to a candidate edge) with a transition term (whether consecutive candidates are connected by a plausible on-network route). The Viterbi decode chooses the whole sequence jointly, resolving ambiguity that no per-point rule can.

How do I set the GPS sigma and transition beta parameters?

σ is the standard deviation of the perpendicular GPS error, typically 4–10 m for open-sky vehicle data and 15–25 m in urban canyons or for pedestrians. β controls how strongly route-vs-straight-line mismatch is penalised and scales with sampling interval: roughly 1–3 for dense 1 Hz traces and 5–15 for sparse taxi data at 30–60 s intervals. Calibrate both against a small hand-matched ground-truth set.

What sampling rate does HMM map-matching need?

The algorithm degrades gracefully as intervals lengthen. For 1–5 s data the emission term dominates and matching is near-exact. Between 30 and 60 s the transition term does the heavy lifting and β must be raised. Beyond about two minutes, gap-fill the trajectory first so the route-distance comparison stays informative.

Must I project to a metric CRS before map-matching?

Yes. The emission Gaussian is defined on perpendicular distance in metres and the transition term compares metric route lengths, so both the trace and the road graph must live in the same metric projected CRS such as a UTM zone. Point-to-line projection and shortest-path lengths computed in raw EPSG:4326 degrees are not metrically consistent and bias every score toward higher-latitude distortion.

How do I handle an observation with no candidate edge in range?

Widen the search radius adaptively before giving up — start at three to four σ and expand once if the candidate set is empty. If still empty, the point is likely an off-road excursion. Mark it unmatched and break the Viterbi lattice at that index so the decoder restarts a fresh sub-path afterwards rather than forcing a spurious long-distance transition across the gap.

Related

• Feeding Cleaned GPS Output into an HMM Map-Matcher — the end-to-end wiring deep-dive with runnable osmnx and Viterbi code.
• GPS Precision & Error Handling — the upstream cleaning stage that produces the drift-corrected trace this matcher consumes.
• Origin-Destination Flow Matrices — the natural downstream consumer of matched, routable edge sequences.
• Coordinate Reference System Mapping — shared metric CRS discipline for both trace and road graph.
• Speed & Acceleration Profiling — infers transport mode to drive per-mode σ and β calibration.

Back to Spatial Indexing & Density Aggregation