The theory behind hidden Markov map matching fits in a paragraph; the implementation has half a dozen places where a reasonable-looking choice produces silently degraded output. This page is the working implementation, with those places called out — it sits under map matching GPS traces to OSM networks in Python Routing Engines & Isochrone Mapping and assumes the trace has already been through the preprocessing in filtering GPS noise before map matching.

When to use this approach

Write the matcher yourself when you need something a hosted endpoint will not give you: per-observation candidate scores for a confidence measure, a custom emission model for an unusual sensor, or matching against a graph no engine hosts. For ordinary vehicle traces against a standard OSM graph, the hosted endpoint described in using the OSRM match endpoint for trace snapping is faster and better tested.

The implementation is still worth understanding even when you use the endpoint, because every parameter the endpoint exposes maps onto something below.

The Viterbi lattice across three observations Each observation column holds three candidate segments. Edges between columns carry transition scores; one edge is pruned entirely because a turn restriction makes the movement illegal. The best-scoring path through the lattice is the matched sequence. Three observations, three candidates each, one illegal transition obs t−1 obs t obs t+1 a b c d e f g h i −inf: no_left_turn The bold path a → d → g wins on the sum of emission and transition scores, not on proximity at any single observation.

Implementation

Candidate generation first. The spatial index is built once for the whole graph, not per trace.

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


class SegmentIndex:
    """Radius queries over road segments, built once and reused."""

    def __init__(self, segments) -> None:
        self.segments = list(segments)
        # Index midpoints; the radius is widened to cover segment half-length
        mids = np.array([s.geometry.interpolate(0.5, normalized=True).coords[0]
                         for s in self.segments])
        self.tree = cKDTree(mids)
        self.half = np.array([s.geometry.length / 2 for s in self.segments])

    def near(self, xy, radius_m: float, cap: int = 10):
        """Candidate segments with perpendicular distance, nearest first."""
        pt = Point(xy)
        # Widen by the longest half-segment so long ways are not missed
        idx = self.tree.query_ball_point(xy, r=radius_m + float(self.half.max()))
        scored = []
        for i in idx:
            d = self.segments[i].geometry.distance(pt)
            if d <= radius_m:
                scored.append((self.segments[i], float(d)))
        scored.sort(key=lambda t: t[1])
        return scored[:cap]

Indexing segment midpoints and widening the query radius is the detail that avoids a subtle recall failure. A KD-tree over midpoints alone misses a long segment whose midpoint is far away but whose nearest point is close, and long segments are exactly the motorway links a trace spends most of its time on.

Then the two scoring functions, both in log space:

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

NEG_INF = -np.inf


def emission_logp(dist_m: float, sigma_m: float) -> float:
    """Half-normal on perpendicular distance."""
    return -0.5 * (dist_m / sigma_m) ** 2 - np.log(sigma_m)


def transition_logp(straight_m: float, route_m: float | None,
                    beta_m: float) -> float:
    """Exponential on how far the on-road distance departs from straight-line."""
    if route_m is None or not np.isfinite(route_m):
        return NEG_INF                # no legal path — including turn restrictions
    return -abs(route_m - straight_m) / beta_m - np.log(beta_m)

Returning -inf rather than a large negative number is deliberate. It makes an illegal transition genuinely impossible rather than merely expensive, so no amount of favourable emission evidence elsewhere can buy it. That is how the turn-restriction set from handling turn restrictions in routing graphs enters the matcher: if the routing query respects restrictions, illegal movements return no route and the transition is -inf automatically.

Finally the decode, with the recovery case handled explicitly:

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


def viterbi(columns, trans_logp):
    """Best path through the lattice; returns None if the lattice dies."""
    scores = {cand: emit for cand, emit in columns[0]}
    back = []

    for t in range(1, len(columns)):
        nxt, ptr = {}, {}
        for cand, emit in columns[t]:
            best, arg = -np.inf, None
            for prev, prev_score in scores.items():
                if prev_score == -np.inf:
                    continue
                s = prev_score + trans_logp(prev, cand, t) + emit
                if s > best:
                    best, arg = s, prev
            if arg is not None and best > -np.inf:
                nxt[cand], ptr[cand] = best, arg
        if not nxt:
            return None               # dead lattice — caller splits the trace
        scores, back = nxt, back + [ptr]

    cur = max(scores, key=scores.get)
    path = [cur]
    for ptr in reversed(back):
        cur = ptr[cur]
        path.append(cur)
    return list(reversed(path))

Skipping predecessors already at -inf is a small optimisation with a large effect on traces that contain restricted junctions, because dead branches otherwise stay in the frontier for the rest of the decode.

Key parameters and tuning

Parameter Recommended value Notes
sigma_m 8 m open, 15–20 m urban The dominant accuracy parameter; prefer per-fix accuracy where reported
beta_m 20–40 m Lower values enforce route plausibility harder and prune more transitions
Candidate radius 50 m Roughly 3σ; wider adds candidates the emission term rejects anyway
Candidate cap 10 Decode cost is quadratic in this; the tail never wins
Index geometry segment midpoints + widened radius Midpoint-only indexing silently misses long segments
Impossible transition -inf A large negative number can still be bought by favourable emissions
Dead-lattice policy split and restart Forcing a path fabricates a route

Integration points

Routing queries in the inner loop. The transition term needs an on-road distance between candidate pairs, which is a routing query per pair per step. This is the dominant cost of the whole matcher, so cache aggressively on the segment pair, and choose a graph backend that makes the query cheap — the comparison in comparing Python graph libraries for OSM is directly relevant, since igraph’s sub-millisecond queries make an otherwise impractical matcher usable.

Confidence output. Because the lattice holds a score for every candidate at every step, the ratio between the best and second-best path score is a usable per-observation confidence. Hosted endpoints do not expose this, and it is the main reason to run your own matcher.

Downstream consumers. The matched edge sequence feeds speed calibration and route-adherence analysis. Both benefit from the confidence measure, because a low-confidence stretch should be excluded from calibration rather than averaged in.

Where a Python matcher spends its time Runtime for a thousand-observation trace is dominated by routing queries for the transition term at 74 percent. Candidate generation is 14 percent and the Viterbi decode itself only 9 percent. Caching route distances by segment pair cuts the dominant stage by roughly three quarters. Runtime by stage, 1 000-observation trace routing queries 74 % candidate generation 14 % Viterbi decode 9 % with pair cache total falls to 26 % of the original wall clock Consecutive observations reuse the same candidate pairs constantly, so a plain dictionary cache converts most queries into lookups. Optimising the decode before the routing queries is the classic wasted effort here — it is under a tenth of the time. Candidate cap against decode time Viterbi decode cost grows with the square of candidates per observation: five candidates take under a second per thousand observations, ten take three, and an uncapped search averaging thirty-four takes thirty-one seconds for identical matched output. Decode cost against candidates kept per observation cap 5 0.9 s per 1 000 obs cap 10 3.1 s — the default cap 20 11.8 s uncapped, ~34 mean 31 s, identical output The output is identical from about cap 10 upward on ordinary vehicle traces, because the emission term makes the remaining candidates unselectable. Verify that on your own hand-matched sample rather than assuming it — a sensor with unusually wide noise moves the point where the cap starts to bite.

Validation checklist

  • Log-space arithmetic never overflows. Run a thousand-observation trace and confirm no score becomes NaN. A NaN means an -inf was added to a +inf somewhere.
  • A restricted junction is respected. Construct a trace that would take an illegal turn under nearest-segment matching and confirm the matcher routes around it.
  • A long segment is found. Place an observation near the end of a two-kilometre way and confirm that way appears among the candidates. If it does not, the index radius is not widened.
  • The dead-lattice path returns None. Feed a trace with an impossible jump and confirm the matcher returns None rather than a fabricated path.
  • Candidate cap does not change output. Compare matched sequences at cap 10 and cap 30 on the hand-matched sample. They should be identical; if not, the cap is too tight for your radius.