Hidden Markov Model Map Matching in Python

When raw GPS coordinates drift, skip, or suffer multipath interference in dense urban canyons, simple nearest-road snapping fails catastrophically — a single outlier ping can snap a vehicle onto a parallel street one block over, corrupting every downstream metric from dwell time to fuel consumption. Hidden Markov Model (HMM) map matching solves this by treating the road network as a set of latent states and GPS observations as noisy emissions from those states. By jointly optimising spatial proximity, topological connectivity, and temporal plausibility, HMM matching delivers production-grade route reconstruction even under severe noise.

This guide provides a complete, step-by-step implementation workflow for mobility engineers, Python GIS developers, and logistics platform builders. For the broader context on spatial trajectory processing, see the parent overview at Trajectory Analysis & Map Matching Techniques.

Prerequisites

Before implementing an HMM matcher, ensure your environment and data pipeline meet these baseline requirements:

• Python 3.10+ with numpy, pandas, scipy, shapely, and requests
• Routing engine: OSRM, Valhalla, or GraphHopper running locally or via managed API
• Input data: Time-ordered GPS traces in GeoJSON, CSV, or Parquet format containing lat, lon, timestamp, and optionally speed/heading
• Network data: OSM-derived road graph with topology (nodes, edges, speed limits, turn restrictions)
• Mathematical background: Markov chains, emission/transition probabilities, and the Viterbi algorithm at a conceptual level

HMM map matching does not require deep learning frameworks. It relies on probabilistic graph traversal, making it deterministic, auditable, and suitable for compliance-heavy fleet operations.

Step-by-Step Workflow

1. Data Ingestion and Temporal Alignment

Raw telemetry arrives asynchronously from telematics control units (TCUs) or mobile SDKs. Sort by timestamp, remove stationary pings (consecutive points within 2 m of each other), and extract kinematic features that later refine probability matrices.

Temporal alignment also demands clock drift handling and timezone normalisation. Always convert timestamps to UTC and compute inter-point deltas (Δt) in seconds. Use timestamp synchronisation across mixed OBD-II and mobile devices as a reference if your fleet mixes hardware types. If Δt exceeds a configurable threshold (e.g. 30 seconds), treat the trace as two independent segments — this prevents the Viterbi decoder from forcing implausible long-range transitions.

import pandas as pd
import numpy as np

def prepare_trace(df: pd.DataFrame, gap_threshold_s: float = 30.0) -> list[pd.DataFrame]:
    """
    Sort, deduplicate, and split a GPS trace on temporal gaps.

    Parameters
    ----------
    df : DataFrame with columns lat, lon, timestamp (UTC ISO-8601 or epoch)
    gap_threshold_s : seconds; gaps larger than this start a new segment

    Returns
    -------
    List of DataFrames, each a contiguous trace segment
    """
    df = df.sort_values("timestamp").reset_index(drop=True)
    df["timestamp"] = pd.to_datetime(df["timestamp"], utc=True)
    df["dt"] = df["timestamp"].diff().dt.total_seconds().fillna(0.0)

    # Drop near-stationary duplicates (< 2 m apart)
    from shapely.geometry import Point
    coords = df[["lon", "lat"]].values
    mask = np.ones(len(df), dtype=bool)
    for i in range(1, len(df)):
        if Point(coords[i]).distance(Point(coords[i - 1])) * 111_320 < 2.0:
            mask[i] = False
    df = df[mask].reset_index(drop=True)
    df["dt"] = df["timestamp"].diff().dt.total_seconds().fillna(0.0)

    # Split on gaps
    split_points = df.index[df["dt"] > gap_threshold_s].tolist()
    boundaries = [0] + split_points + [len(df)]
    return [df.iloc[boundaries[i]:boundaries[i + 1]] for i in range(len(boundaries) - 1)]

Expected output shape: a list of DataFrames, each with columns lat, lon, timestamp, dt. A trace of 500 pings typically splits into 1–5 segments depending on signal quality.

2. Candidate Road Segment Generation

For each GPS observation, query the routing engine or a spatial index (e.g. scipy.spatial.cKDTree over OSM edge centroids) to retrieve k nearest road segments within radius R (typically 50–150 m). Each candidate records:

• Edge ID and full geometry
• Perpendicular distance from the GPS point to its projected point on the edge
• Road class, speed limit, and one-way status

When projecting points onto linestrings, use Shapely’s nearest_points to compute exact snap locations and bearing vectors. If your fleet includes mixed vehicle types, filter candidates against vehicle-specific constraints (bridge clearance, turn radius) — the multi-modal considerations are covered in Multi-Modal Route Matching for Mixed Fleets. For strategies on aligning sensor-derived heading with road network directionality, see Directionality & Heading Synchronisation.

from scipy.spatial import cKDTree
from shapely.geometry import Point, LineString
from shapely.ops import nearest_points

def build_candidate_index(edges: list[dict]) -> cKDTree:
    """Build a KD-tree over edge midpoints for fast radius queries."""
    centroids = [
        LineString(e["geometry"]).interpolate(0.5, normalized=True).coords[0]
        for e in edges
    ]
    return cKDTree(centroids)

def get_candidates(
    point: tuple[float, float],
    edges: list[dict],
    tree: cKDTree,
    radius_m: float = 100.0,
) -> list[dict]:
    """
    Return candidate edges within radius_m of the GPS point.

    Parameters
    ----------
    point   : (lon, lat) in WGS-84 — note longitude-first order
    radius_m: search radius in metres (approximate; 1 deg lat ≈ 111 320 m)
    """
    radius_deg = radius_m / 111_320
    idxs = tree.query_ball_point(point, radius_deg)
    candidates = []
    for i in idxs:
        edge = edges[i]
        geom = LineString(edge["geometry"])
        snap, _ = nearest_points(geom, Point(point))
        dist_m = Point(point).distance(snap) * 111_320
        candidates.append({
            "edge_id": edge["id"],
            "geometry": geom,
            "snap_point": snap,
            "distance_m": dist_m,
            "speed_limit_ms": edge.get("speed_limit_kmh", 50) / 3.6,
        })
    return sorted(candidates, key=lambda c: c["distance_m"])

3. Probability and Math Model

The HMM framework requires two probability distributions. Working in log-space throughout is non-negotiable — floating-point underflow silently produces all-zero probability columns, which causes the Viterbi decoder to emit arbitrary paths with no warning.

Emission Probability

Models how likely a GPS ping is, given the vehicle occupies a specific road segment. Modelled as a zero-mean Gaussian over the perpendicular snap distance d:

P_emission(d | σ) = (1 / (σ√(2π))) · exp(−0.5 · (d/σ)²)

σ represents GPS measurement error — typically 5–15 m, but calibrate it from the device’s reported HDOP if available (σ ≈ 5 + HDOP × 3 metres).

Transition Probability

Models how likely it is that the vehicle moves from road segment i at time t−1 to segment j at time t. The routing engine returns the network distance D_route(i→j); the expected travel distance is v_avg × Δt. The probability decays exponentially with the absolute mismatch:

P_transition(i→j) = exp(−β · |D_route(i→j) − v_avg · Δt|)

β (default 2.0 m⁻¹) controls sensitivity. Speed profiling from raw GPS coordinates can supply per-segment v_avg estimates, tightening these bounds considerably.

import numpy as np

def compute_log_emission(distances_m: np.ndarray, sigma: float = 10.0) -> np.ndarray:
    """
    Log-space Gaussian emission probabilities.

    Parameters
    ----------
    distances_m : (N,) perpendicular snap distances in metres
    sigma       : GPS error standard deviation in metres

    Returns
    -------
    (N,) log-probabilities; higher (less negative) is more likely
    """
    return -0.5 * (distances_m / sigma) ** 2 - np.log(sigma * np.sqrt(2 * np.pi))


def compute_log_transition(
    route_distances_m: np.ndarray,
    delta_t: float,
    avg_speed_ms: float,
    beta: float = 2.0,
) -> np.ndarray:
    """
    Log-space transition probabilities based on routing vs. kinematic distance.

    Parameters
    ----------
    route_distances_m : (N_prev, N_curr) shortest-path distances from routing engine
    delta_t           : elapsed time in seconds between observations
    avg_speed_ms      : expected average speed in m/s for this edge pair
    beta              : decay rate in m⁻¹ (higher = stricter distance matching)

    Returns
    -------
    (N_prev, N_curr) log-transition matrix
    """
    expected_dist_m = avg_speed_ms * delta_t
    diff = np.abs(route_distances_m - expected_dist_m)
    return -beta * diff

Routing Engine Integration

Use OSRM’s /table endpoint to batch-request all N_prev × N_curr route distances in a single HTTP call rather than making individual /route requests:

import requests

def osrm_distance_matrix(
    sources: list[tuple[float, float]],
    destinations: list[tuple[float, float]],
    base_url: str = "http://localhost:5000",
) -> np.ndarray:
    """
    Fetch a routing distance matrix from a local OSRM instance.

    Parameters
    ----------
    sources      : list of (lon, lat) snap points — OSRM expects lon/lat order
    destinations : list of (lon, lat) snap points

    Returns
    -------
    (len(sources), len(destinations)) distance matrix in metres;
    unreachable pairs are set to 1e9 (effectively -inf in log-space)
    """
    coords_str = ";".join(
        f"{lon},{lat}"
        for lon, lat in (sources + destinations)
    )
    src_idxs = ";".join(str(i) for i in range(len(sources)))
    dst_idxs = ";".join(str(i + len(sources)) for i in range(len(destinations)))

    resp = requests.get(
        f"{base_url}/table/v1/driving/{coords_str}",
        params={"sources": src_idxs, "destinations": dst_idxs, "annotations": "distance"},
        timeout=5,
    )
    resp.raise_for_status()
    data = resp.json()

    matrix = np.array(data["distances"], dtype=float)
    matrix[matrix is None] = 1e9  # unreachable
    matrix = np.where(matrix == 0.0, 1e9, matrix)  # zero-distance self-loops
    return matrix

Coordinate order gotcha: OSRM uses (longitude, latitude) throughout its API — the opposite of most GeoJSON tools. Passing (lat, lon) silently produces wrong distances with no error.

4. Viterbi Decoding and Path Reconstruction

With log-probability matrices populated, the Viterbi algorithm finds the most likely sequence of road segments. The algorithm runs in O(T × N²) time, where T is the number of GPS observations and N is the average candidate count per observation. Avoid nested Python loops; use numpy broadcasting throughout.

def viterbi_decode(
    log_emissions: np.ndarray,
    log_transitions: list[np.ndarray],
) -> tuple[np.ndarray, float]:
    """
    Vectorised Viterbi decoder for HMM map matching.

    Parameters
    ----------
    log_emissions   : (T, N) — log emission probability for each observation/candidate
    log_transitions : list of T-1 arrays, each (N_prev, N_curr)
                      — log transition probability between consecutive observations

    Returns
    -------
    path            : (T,) array of candidate indices (the optimal state sequence)
    log_prob        : scalar log-probability of the optimal path
    """
    T = log_emissions.shape[0]
    N = log_emissions.shape[1]
    delta = np.full((T, N), -np.inf)
    psi = np.zeros((T, N), dtype=int)

    # Initialise from first observation
    delta[0] = log_emissions[0]

    # Forward pass
    for t in range(1, T):
        # log_transitions[t-1] has shape (N_prev, N_curr)
        # delta[t-1] has shape (N_prev,) — broadcast over columns
        scores = delta[t - 1, :, np.newaxis] + log_transitions[t - 1]  # (N_prev, N_curr)
        psi[t] = np.argmax(scores, axis=0)                              # (N_curr,)
        delta[t] = np.max(scores, axis=0) + log_emissions[t]           # (N_curr,)

    # Backtrack
    path = np.zeros(T, dtype=int)
    path[-1] = int(np.argmax(delta[-1]))
    for t in range(T - 2, -1, -1):
        path[t] = psi[t + 1, path[t + 1]]

    return path, float(delta[-1, path[-1]])

The algorithm’s formal foundation is the Newson & Krumm paper “Hidden Markov Map Matching Through Noise and Sparseness” (2009), summarised on Wikipedia’s map matching overview. The implementation in Building an HMM-based map matcher with OSRM and Python integrates these components into a complete, runnable class.

5. Post-Processing and Production Validation

The decoded path yields a sequence of edge IDs. Convert these into actionable telemetry:

1. Geometric smoothing: Interpolate along matched edges using cumulative distance ratios to reconstruct a continuous polyline.
2. Timestamp alignment: Reassign timestamps proportionally along the matched geometry to preserve temporal fidelity for downstream stop detection and dwell-time analytics.
3. Confidence scoring: Compute the normalised log-probability of the optimal path per observation: score = log_prob / T. Paths with per-observation scores below −15.0 indicate severe GPS degradation or unmapped infrastructure — flag them for manual QA rather than silently propagating bad data.
4. Outlier removal validation: Cross-check that the matched speed profile is consistent with speed_limit × 1.3 along each edge. Matched segments that imply impossible speeds are a reliable signal that the trace had an uncorrected outlier upstream.

Always validate against ground-truth routes or known depot-to-destination corridors. Implement automated regression tests that feed synthetic traces — with known noise profiles — into your pipeline to verify decoder stability before deploying to production fleets.

Routing Engine Integration Notes

Config flags to verify before deployment:

• OSRM: compile with --algorithm MLD for live traffic support; --max-table-size 1000 prevents request timeouts with large candidate grids.
• Valhalla: set "use_restrictions": 1 in the costing options; default is 0.
• GraphHopper: enable turn_costs=true in config.yml; without it, illegal U-turns are not penalised.

Cache D_route values keyed by (edge_id_prev, edge_id_curr, vehicle_type) tuples. For a fleet of 100 vehicles on an urban grid, the warm-cache hit rate typically exceeds 70 % within a day of operation.

Operational Troubleshooting

Deployment Checklist

Validate all GPS coordinates are in WGS-84 (EPSG:4326) before spatial indexing
Set σ dynamically from device-reported HDOP; fall back to 10 m if unavailable
Verify routing engine coordinate order (lon,lat for OSRM; lat,lon for Valhalla/GraphHopper)
Confirm turn restrictions are active in the routing engine config
Use log-space arithmetic and apply delta rescaling to prevent underflow
Implement sliding-window Viterbi for traces longer than 500 observations
Add circuit breakers and fallback distances for routing API failures
Log unmatched segments with per-observation confidence scores for QA review
Benchmark decoder latency against your SLA — target: under 50 ms per 100-point trace on standard cloud instances
Run regression tests with synthetic traces at multiple noise levels before production deployment

By combining Gaussian emission probabilities, routing-derived transition costs, and log-space Viterbi decoding, this workflow delivers deterministic, auditable route reconstruction at sub-10-metre accuracy across thousands of fleet vehicles. The probabilistic foundation gives compliance teams a transparent paper trail — every matched segment has a quantified confidence score — making HMM matching the preferred approach for regulated logistics, insurance telematics, and municipal mobility planning.

Related

• Building an HMM-based map matcher with OSRM and Python — complete, runnable implementation integrating all components above
• Speed Profiling from Raw GPS Coordinates — supplying per-segment speed priors for tighter transition probabilities
• Directionality & Heading Synchronisation — aligning heading vectors with road network directionality before candidate filtering
• Stop Detection & Dwell-Time Analytics — downstream consumer of matched trajectories; dwell accuracy depends directly on matching quality
• GPS Data Preprocessing & Cleaning Fundamentals — upstream CRS normalisation, outlier removal, and timestamp alignment that feed into the matcher

Parent: Trajectory Analysis & Map Matching Techniques