Map Matching GPS Traces to OSM Networks

Almost everything useful you can learn from a fleet’s telemetry requires first knowing which road each GPS point was on, and a GPS point does not know that. Map matching is the step that converts a sequence of noisy coordinates into a sequence of graph edges, and it is the foundation under speed calibration, route-adherence analysis, and any attempt to compare planned against actual routes. This page covers the technique as part of Python Routing Engines & Isochrone Mapping, and its output feeds directly into the telematics work described in speed profile calibration for heavy vehicles.

The naive approach — snap each point to the nearest road — is correct often enough to look like it works and wrong often enough to be useless. The failure is systematic rather than random: it concentrates exactly where roads run close together, which is where the interesting questions are.

Prerequisites

Libraries. networkx or igraph for the graph, scipy for the spatial index, numpy for the probability arrays, and shapely for the perpendicular-distance geometry.

pip install networkx>=3.0 scipy>=1.11 numpy>=1.24 shapely>=2.0 geopandas>=0.14

A projected graph. Every distance in a matcher is a real distance in metres — emission distances, transition distances, search radii. Working in degrees makes all three meaningless, so project the graph and the trace to a local metric CRS before anything else.

A graph with turn restrictions. Optional but transformative. Without them the matcher will happily produce paths through prohibited movements, which then contaminate any behavioural analysis built on the output.

Traces with timestamps. Not just coordinates. The transition model needs elapsed time to distinguish a plausible movement from an implausible one, and a trace without timestamps can only be matched geometrically.

Conceptual architecture

Map matching is a sequence-labelling problem wearing a geospatial costume. Each GPS observation has several candidate road segments; the job is to pick one per observation such that the whole sequence is coherent.

The hidden Markov formulation gives each of those two intuitions a number. Emission probability says how likely an observation is given that the vehicle was on a particular segment — a function of perpendicular distance and the GPS noise model. Transition probability says how likely it is to move from one segment to another between consecutive observations — a function of how closely the on-road route distance matches the straight-line distance between the two fixes. A pair of points that are 200 metres apart in a straight line but 3 kilometres apart by road describes a transition that almost certainly did not happen.

Viterbi decoding then finds the single most likely sequence across the whole trace, which is what makes the result coherent rather than a series of independent guesses.

Step-by-step implementation

1. Filter the raw trace

Preprocessing is not optional and is where most of the practical accuracy comes from. Raw traces contain stationary jitter, single-fix outliers, and duplicated timestamps, all of which the matcher will faithfully try to explain.

# requires: numpy, pandas (pip install numpy pandas)
import numpy as np
import pandas as pd


def clean_trace(df: pd.DataFrame, *, max_speed_mps: float = 60.0,
                min_move_m: float = 5.0) -> pd.DataFrame:
    """Drop implausible jumps and collapse stationary jitter."""
    df = df.sort_values("timestamp").drop_duplicates("timestamp").reset_index(drop=True)

    dx = np.diff(df["x"].to_numpy(), prepend=df["x"].iloc[0])
    dy = np.diff(df["y"].to_numpy(), prepend=df["y"].iloc[0])
    step = np.hypot(dx, dy)
    dt = df["timestamp"].diff().dt.total_seconds().fillna(1.0).clip(lower=0.1)

    # A fix implying an impossible speed is a bad fix, not a fast vehicle
    df = df[step / dt <= max_speed_mps].reset_index(drop=True)

    # Collapse runs where the vehicle has not meaningfully moved
    keep, last = [0], 0
    xs, ys = df["x"].to_numpy(), df["y"].to_numpy()
    for i in range(1, len(df)):
        if np.hypot(xs[i] - xs[last], ys[i] - ys[last]) >= min_move_m:
            keep.append(i)
            last = i
    return df.iloc[keep].reset_index(drop=True)

Collapsing stationary runs matters more than it looks. A vehicle waiting at a signal for ninety seconds produces ninety fixes scattered over a few metres; left in, they dominate the trace’s point count and give the matcher ninety opportunities to wander onto a side street.

2. Generate candidates and score emissions

# requires: numpy, scipy (pip install numpy scipy)
import numpy as np
from scipy.spatial import cKDTree


def emission_logprob(dist_m: np.ndarray, sigma_m: float = 8.0) -> np.ndarray:
    """Log-likelihood of an observation given a candidate segment."""
    # Half-normal on perpendicular distance; sigma is the GPS noise scale
    return -0.5 * (dist_m / sigma_m) ** 2 - np.log(sigma_m)


def candidates(point_xy, seg_tree: cKDTree, seg_index, radius_m: float = 50.0):
    """Segments within radius of an observation, with perpendicular distances."""
    hits = seg_tree.query_ball_point(point_xy, r=radius_m)
    out = []
    for i in hits:
        seg = seg_index[i]
        d = seg.geometry.distance(point_xy)   # perpendicular, in metres
        out.append((seg, float(d)))
    return out

sigma_m is the single most consequential parameter in the matcher. Set it too low and the model refuses to accept a fix that is genuinely 20 metres off; too high and it stops discriminating between adjacent roads. Eight metres suits consumer GPS in open conditions; urban canyons justify 15 to 20.

3. Score transitions

# requires: numpy (pip install numpy)
import numpy as np


def transition_logprob(great_circle_m: float, route_m: float,
                       beta_m: float = 30.0) -> float:
    """Log-likelihood of moving between two candidate segments."""
    if not np.isfinite(route_m):
        return -np.inf                      # no legal path — impossible transition
    # Exponential on the discrepancy between straight-line and on-road distance
    return -abs(route_m - great_circle_m) / beta_m - np.log(beta_m)

Returning -inf for an unreachable pair is what lets turn restrictions enter the model. If the routing query that produces route_m respects the restriction set — as the edge-keyed search in handling turn restrictions in routing graphs does — then illegal movements simply have no finite transition and are never chosen.

4. Decode with Viterbi

# requires: numpy (pip install numpy)
import numpy as np


def viterbi(obs_candidates, trans_fn):
    """Most likely candidate sequence across all observations."""
    scores = {c: emission for c, emission in obs_candidates[0]}
    back = [{}]

    for t in range(1, len(obs_candidates)):
        nxt, ptr = {}, {}
        for cand, emission in obs_candidates[t]:
            best, best_prev = -np.inf, None
            for prev, prev_score in scores.items():
                s = prev_score + trans_fn(prev, cand, t) + emission
                if s > best:
                    best, best_prev = s, prev
            if best > -np.inf:
                nxt[cand], ptr[cand] = best, best_prev
        if not nxt:
            # No survivable transition — break the trace and restart here
            return None
        scores, back = nxt, back + [ptr]

    path, cur = [], max(scores, key=scores.get)
    for t in range(len(back) - 1, 0, -1):
        path.append(cur)
        cur = back[t][cur]
    path.append(cur)
    return list(reversed(path))

The empty-nxt case deserves attention rather than an exception. Traces legitimately contain gaps — a tunnel, a dead battery, a ferry crossing — and the right response is to split the trace and match the segments independently rather than to force a path across the gap.

Configuration reference

Production optimization and scaling

Precompute the segment index once. The KD-tree over segment geometries is expensive to build and constant across traces. Build it at service start, not per trace, and reuse it for every match.

Cap and prune candidates aggressively. Viterbi cost grows with the square of candidates per observation. Keeping the ten nearest rather than everything within the radius typically halves runtime with no measurable accuracy change, because the eleventh candidate essentially never wins.

Cache route distances. The transition term needs an on-road distance between candidate pairs, and consecutive observations reuse the same pairs constantly. A dictionary keyed on the segment pair plus a rounded time bucket converts most of those queries into lookups.

Batch by region. Traces from the same area touch the same part of the graph, so grouping them keeps the working set small. This matters most when matching against a graph too large to hold entirely in memory.

Decide early whether to write a matcher at all. OSRM’s /match and Valhalla’s /trace_attributes implement the same hidden Markov formulation, are heavily exercised, and run at speeds a Python implementation will not reach. The reasons to write your own are specific: you need the per-observation candidate scores rather than only the winning path, you are matching against a graph the engine does not host, or your sensor has a noise profile the engine’s fixed emission model handles badly. Each of those is real, and none of them is “we wanted to understand it” — the understanding is worth having, but it does not need to reach production.

Match against the contemporaneous graph. Version the graph alongside the traces and match each trace against the build that was live when it was recorded, using the artefact retention described in incremental graph rebuilds without full reprocessing. Otherwise every road built since the recording shows up as driver deviation.

Validation and testing

The only trustworthy validation is a hand-matched sample. Everything else measures self-consistency rather than correctness.

# requires: numpy (pip install numpy)
import numpy as np


def match_quality(matched: list, truth: list) -> dict:
    """Compare a matched edge sequence against a hand-verified one."""
    m, t = set(matched), set(truth)
    correct = len(m & t)
    return {
        "edge_precision": round(correct / max(len(m), 1), 4),
        "edge_recall": round(correct / max(len(t), 1), 4),
        "length_error_pct": round(
            100 * abs(sum(e.length for e in matched) - sum(e.length for e in truth))
            / max(sum(e.length for e in truth), 1e-9), 2
        ),
    }

A hundred hand-matched traces is enough to calibrate sigma_m and beta_m properly and to establish a regression baseline. Below about ninety-five percent edge recall on that sample, downstream speed statistics are being computed on the wrong roads often enough to matter.

Troubleshooting

Related

• NetworkX shortest-path algorithms for logistics — the routing queries that supply the transition term’s on-road distances
• Handling turn restrictions in routing graphs — making illegal movements impossible transitions rather than merely expensive ones
• Speed profile calibration for heavy vehicles — the telematics ground-truthing that consumes matched traces
• Deploying OSRM with Docker for local routing — hosting the engine whose match endpoint is the alternative to a hand-written matcher
• Comparing Python graph libraries for OSM — choosing the graph backend that makes the inner routing queries fast enough