Building an HMM-based Map Matcher with OSRM and Python

This page answers one specific engineering question: given raw GPS coordinates from a vehicle, how do you wire OSRM’s HTTP API into a complete Hidden Markov Model pipeline that outputs a topologically valid road sequence? It extends the Hidden Markov Model map matching in Python guide by turning the mathematical framework into a runnable, self-contained Python class. Naive nearest-road snapping fails on parallel streets, roundabouts, and multipath-heavy urban canyons because it considers each GPS ping in isolation. The HMM approach — rooted in Newson & Krumm’s 2009 formulation — encodes both geometric fit (how close is the ping to the road?) and network reachability (can a vehicle actually drive from the previous candidate to this one?), then uses the Viterbi algorithm to decode the globally optimal path through the road graph in a single forward pass.

Compatibility & Configuration Requirements

The following versions and flags are required. Mismatches produce silent failures — OSRM returns "code": "Ok" with an empty waypoint list rather than an error when coordinate order is wrong.

Start OSRM with the correct flags:

docker run -t -v "${PWD}:/data" osrm/osrm-backend \
    osrm-routed --algorithm mld --max-table-size 500 /data/map.osrm

For heavy freight or micromobility fleets, compile a custom .lua profile (osrm-extract -p /opt/car.lua) that encodes turn restrictions, axle-weight limits, and speed caps specific to your vehicle class. Using the default car profile for a 26-tonne HGV will produce routing distances that underestimate real journey times and contaminate transition probabilities.

Avoid pandas in the hot path: type coercion and index alignment add 15–20 ms per trajectory. Use plain numpy arrays throughout.

Production-Ready Implementation

The class below handles candidate fetching, log-space probability computation, Viterbi decoding, and path reconstruction in a single self-contained unit. Read the inline comments — every parameter choice reflects a deliberate trade-off between computational cost and matching accuracy.

import numpy as np
import requests
from typing import List, Tuple, Dict


class OSRMMatcher:
    """
    HMM map matcher backed by a local OSRM routing engine.

    Parameters
    ----------
    osrm_url : str
        Base URL of your OSRM instance, e.g. "http://localhost:5000".
    sigma_z : float
        GPS accuracy (metres). Controls how strongly the emission term
        penalises candidates that are far from the raw ping.
        Typical range: 5–15 m. Raise for urban-canyon devices; lower
        only with survey-grade RTK receivers.
    beta : float
        Routing-distance decay constant (metres). Controls how strongly
        the transition term penalises long routing detours between
        consecutive candidates.
        Typical range: 200–500 m. Raise for sparse sampling (>30 s
        intervals) or highway corridors; lower for dense urban delivery
        routes with frequent turns.
    """

    def __init__(
        self,
        osrm_url: str = "http://localhost:5000",
        sigma_z: float = 10.0,
        beta: float = 200.0,
    ):
        self.osrm_url = osrm_url
        self.sigma_z = sigma_z
        self.beta = beta
        # Reuse the TCP connection across all API calls — critical for
        # low-latency batch processing.
        self.session = requests.Session()

    def _fetch_candidates(
        self,
        coords: List[Tuple[float, float]],
        radius: int = 50,
        k: int = 3,
    ) -> List[List[Dict]]:
        """
        Query OSRM /nearest for each GPS ping.

        Parameters
        ----------
        coords : list of (lon, lat) tuples in WGS 84 / EPSG:4326.
            NOTE: longitude first — reversed from (lat, lon) conventions.
        radius : int
            Search radius in metres. Increase to 100+ for sparse road
            networks (rural or motorway-only areas).
        k : int
            Number of candidate road nodes per ping. k=3 is a good
            default; increase to 5 for complex junctions or roundabouts.

        Returns
        -------
        List of candidate lists; each candidate is
        {"loc": [lon, lat], "dist": float (metres)}.
        """
        all_candidates = []
        for lon, lat in coords:
            url = (
                f"{self.osrm_url}/nearest/v1/driving/{lon},{lat}"
                f"?number={k}&radius={radius}"
            )
            resp = self.session.get(url, timeout=5).json()
            if resp.get("code") != "Ok" or not resp.get("waypoints"):
                # Append empty list so index alignment is preserved.
                # Viterbi will treat this step as a gap.
                all_candidates.append([])
                continue
            all_candidates.append(
                [
                    {"loc": wp["location"], "dist": wp["distance"]}
                    for wp in resp["waypoints"]
                ]
            )
        return all_candidates

    def _emission_log_prob(self, distance: float) -> float:
        """
        Log-space emission probability.

        P(obs | state) ∝ exp(-0.5 * (d / σ_z)²)
        log P = -0.5 * (d / σ_z)²

        Using log-space throughout avoids floating-point underflow when
        multiplying many small probabilities across long trajectories.
        """
        return -0.5 * (distance / self.sigma_z) ** 2

    def _transition_log_prob(self, route_dist: float) -> float:
        """
        Log-space transition probability.

        P(s_t | s_{t-1}) ∝ exp(-route_dist / β)
        log P = -route_dist / β

        route_dist is the OSRM shortest-path distance in metres.
        np.inf input (no valid route) correctly maps to -inf, which
        Viterbi treats as an impossible transition.
        """
        return -route_dist / self.beta

    def _pairwise_route_distances(
        self,
        candidates_a: List[Dict],
        candidates_b: List[Dict],
    ) -> np.ndarray:
        """
        Fetch routing distances between every pair of candidates at
        two consecutive time steps.

        Uses the OSRM /table endpoint for O(1) HTTP overhead instead of
        O(N²) individual /route calls. Falls back to np.inf on API failure
        so the Viterbi step still runs — the path simply avoids that pair.

        Returns
        -------
        dist_matrix : ndarray of shape (len(candidates_a), len(candidates_b))
            Routing distances in metres. np.inf where no route exists.
        """
        n_a, n_b = len(candidates_a), len(candidates_b)
        if n_a == 0 or n_b == 0:
            return np.full((max(n_a, 1), max(n_b, 1)), np.inf)

        # Build the combined coordinate list: sources first, then
        # destinations. OSRM /table uses 0-based indices.
        all_coords = candidates_a + candidates_b
        coord_str = ";".join(
            f"{c['loc'][0]},{c['loc'][1]}" for c in all_coords
        )
        src_indices = ",".join(str(i) for i in range(n_a))
        dst_indices = ",".join(str(i) for i in range(n_a, n_a + n_b))

        url = (
            f"{self.osrm_url}/table/v1/driving/{coord_str}"
            f"?sources={src_indices}&destinations={dst_indices}"
            "&annotations=distance"
        )
        try:
            resp = self.session.get(url, timeout=10).json()
            if resp.get("code") != "Ok":
                return np.full((n_a, n_b), np.inf)
            # /table returns distances in metres (when annotations=distance).
            matrix = np.array(resp["distances"], dtype=float)
            # OSRM encodes unreachable pairs as null (Python None).
            matrix = np.where(matrix is None, np.inf, matrix)
            return matrix
        except Exception:
            return np.full((n_a, n_b), np.inf)

    def match(self, coords: List[Tuple[float, float]]) -> List[Dict]:
        """
        Run the full HMM pipeline and return the matched road sequence.

        Parameters
        ----------
        coords : list of (lon, lat) tuples, WGS 84 / EPSG:4326.

        Returns
        -------
        List of dicts, one per input coordinate:
        {"lon": float, "lat": float, "distance_to_road": float}
        where lon/lat is the matched road-node position and
        distance_to_road is the perpendicular distance in metres.
        """
        candidates = self._fetch_candidates(coords)
        n_steps = len(coords)
        # Maximum candidate count across all steps — used to size arrays.
        k_max = max((len(c) for c in candidates), default=1)

        # delta[i] holds the log-probability of the best path ending at
        # candidate i of the current step. Initialise from step 0 emissions.
        delta = np.full(k_max, -np.inf)
        for i, c in enumerate(candidates[0]):
            delta[i] = self._emission_log_prob(c["dist"])

        # psi[t, i] stores the predecessor candidate index at step t-1 that
        # maximises the path reaching candidate i at step t.
        psi = np.zeros((n_steps, k_max), dtype=int)

        for t in range(1, n_steps):
            new_delta = np.full(k_max, -np.inf)
            cands_prev = candidates[t - 1]
            cands_curr = candidates[t]

            if not cands_prev or not cands_curr:
                # Gap in the trace — preserve previous delta, skip step.
                psi[t] = np.arange(k_max)
                delta = new_delta
                continue

            dist_matrix = self._pairwise_route_distances(cands_prev, cands_curr)

            for j in range(len(cands_curr)):
                emission = self._emission_log_prob(cands_curr[j]["dist"])
                # Vectorised: score all predecessors simultaneously.
                trans = np.array(
                    [self._transition_log_prob(dist_matrix[i, j])
                     for i in range(len(cands_prev))]
                )
                scores = delta[: len(cands_prev)] + trans + emission
                best_prev = int(np.argmax(scores))
                psi[t, j] = best_prev
                new_delta[j] = scores[best_prev]

            delta = new_delta

        # --- Backtrack ---
        path = np.zeros(n_steps, dtype=int)
        last_k = len(candidates[-1]) if candidates[-1] else 1
        path[-1] = int(np.argmax(delta[:last_k]))
        for t in range(n_steps - 2, -1, -1):
            path[t] = psi[t + 1, path[t + 1]]

        # --- Reconstruct output ---
        result = []
        for t in range(n_steps):
            cands = candidates[t]
            idx = path[t]
            if cands and idx < len(cands):
                c = cands[idx]
                result.append(
                    {
                        "lon": c["loc"][0],
                        "lat": c["loc"][1],
                        "distance_to_road": c["dist"],
                    }
                )
            else:
                # Gap: fall back to the raw ping.
                result.append(
                    {"lon": coords[t][0], "lat": coords[t][1], "distance_to_road": None}
                )
        return result

Execution & Tuning Guidelines

Instantiate and call match() with a list of (longitude, latitude) tuples sourced directly from your telematics device or preprocessed GPS stream:

matcher = OSRMMatcher(
    osrm_url="http://localhost:5000",
    sigma_z=10.0,   # metres — adjust to your device's stated CEP
    beta=300.0,     # metres — adjust to your average inter-ping distance * 10
)

# coords must be (lon, lat), WGS 84 — not (lat, lon).
raw_trace = [
    (13.388860, 52.517037),
    (13.397634, 52.529407),
    (13.428555, 52.523219),
]
matched = matcher.match(raw_trace)
for point in matched:
    print(point)
# {"lon": 13.38885, "lat": 52.51703, "distance_to_road": 3.2}
# {"lon": 13.39763, "lat": 52.52940, "distance_to_road": 1.8}
# ...

sigma_z tuning:

• Above 15 m: the matcher trusts OSRM’s road topology over raw ping proximity. Use when your fleet devices report HDOP > 3 or operate in dense urban canyons where multipath inflates reported accuracy by 10–30 m.
• Below 5 m: emission dominates and the decoder snaps to whichever node is geometrically closest. Reserve for survey-grade RTK receivers (CEP < 2 m). Using this setting with a standard 5 Hz OBD-II dongle will break topology at parallel roads.

beta tuning:

• Above 500 m: allows long routing detours between consecutive pings. Required when sampling rate drops below one ping per 30 seconds — a common failure mode in tunnels or when power-saving firmware buffers points. Combine with GPS preprocessing outlier removal to prevent spurious candidates from being selected across those long gaps.
• Below 100 m: imposes strict topological continuity. Best for high-frequency (1–5 Hz) delivery routes inside dense urban grids where every ping should be reachable from the previous one via a short route.

Replacing per-pair /route calls with /table: The implementation above already uses /table with annotations=distance. If you are on OSRM 5.24 or earlier where annotations=distance is not available, fall back to annotations=duration and multiply by the road speed limit — an approximation sufficient for transition scoring. Never use individual /route calls in a nested loop for production workloads; the O(k²) HTTP overhead makes a 10-ping trace with k=3 candidates issue 90 separate requests.

Numerical stability: When merging evidence from multiple hypotheses (for example, combining a Kalman-smoothed position estimate from Kalman filtering for GPS noise reduction with the emission score), use numpy.logaddexp to sum log-probabilities without leaving log-space. Direct addition in log-space is correct for the Viterbi max-product step, but summation across hypotheses requires logaddexp to avoid underflow.

Common Pitfalls

• Coordinate order reversal: OSRM’s HTTP API expects longitude,latitude throughout — in the URL path, in /table sources/destinations, and in responses. Swapping to latitude,longitude (the default in many Python GIS libraries such as geopandas) places your queries in the ocean off the west coast of Africa. The API returns "code": "Ok" with distance values in the thousands of kilometres, which silently poisons every transition probability.

• max-table-size exhaustion: The default OSRM max-table-size is 100, which limits a /table call with k=3 candidates to a maximum of 16 time-step pairs before the matrix exceeds capacity ((sources + destinations) × k ≤ 100). For trajectories longer than ~15 pings at k=3, either restart OSRM with --max-table-size 500 or implement sliding-window chunking with a two-step overlap so Viterbi can stitch adjacent windows at their shared boundary.

• Empty candidate set propagation: When OSRM returns no candidates for a ping — because the vehicle was off-road, in a car park not in the OSM extract, or beyond the search radius — the current implementation preserves the raw ping and moves on. If two or more consecutive pings have empty candidate sets, the backtrack step will produce a run of None distances. Detect these gaps before passing output to downstream stop detection or speed profiling stages; treating a gap-filled segment as a continuous trajectory inflates dwell times and corrupts speed calculations.

Related

• Parent: Hidden Markov Model Map Matching in Python — full mathematical derivation of the emission and transition models, log-space normalisation, and state-space pruning strategies
• Trajectory Analysis & Map Matching Techniques — overview of geometric, probabilistic, and ML-based matching approaches and when to choose each
• GPS Data Preprocessing & Cleaning Fundamentals — outlier removal, Kalman filtering, and CRS normalisation that feed clean coordinates into this matcher
• Speed Profiling from Raw GPS Coordinates — downstream analysis that consumes the matched road sequence produced here
• Stop Detection & Dwell Time Analytics — uses the matched trajectory to identify stationary periods and classify visited locations