Feeding cleaned GPS output into an HMM map-matcher

This guide wires a drift-corrected trajectory into a Hidden Markov Model map-matcher end to end: you load the cleaned trace, download and project an osmnx drive graph, run an emission-plus-transition HMM decoded by Viterbi, and get back matched edge IDs and snapped geometry. The recurring failure is a mismatch of coordinate systems — the trace, the road graph, and every distance in the scoring functions must share one metric CRS, or the emission Gaussian and the route-distance transition both silently distort. Everything below keeps that discipline explicit.

Why this happens

A map-matcher does not clean data; it assumes the data is already clean. If you pipe a raw GPS log straight into an HMM, residual multipath jitter inflates the perpendicular distance from each observation to its true edge, which forces you to raise the GPS standard deviation σ to compensate. A large σ flattens the emission term, so the matcher can no longer distinguish the correct edge from a parallel road or an overpass — exactly the ambiguity the HMM exists to resolve.

The correct input is a trajectory that has already passed through drift correction: resampled to a stable cadence, kinematically filtered, and projected to a metric CRS. The full statistical machinery of candidate scoring and Viterbi decoding is covered in the parent Map-Matching to Road Networks reference, which itself sits inside the broader Spatial Indexing & Density Aggregation section. This page focuses on the wiring: getting a clean frame and an OSM graph into a matcher without breaking metric-CRS discipline.

Core matching pipeline

The four steps below take you from a cleaned DataFrame to a matched path. Each step is guarded against the edge cases that break real fleet data — empty frames, points with no candidate, and disconnected graphs.

HMM Map-Matcher Wiring Four sequential steps: Load Cleaned Trace, Project OSM Graph, HMM Viterbi Match, and Return & Validate, connected by arrows. 1. Load Cleaned Trace 2. Project OSM Drive Graph 3. HMM Viterbi Match 4. Return & Validate
  1. Load the cleaned trace. Read the drift-corrected trajectory and confirm it carries metric x_m/y_m coordinates (or WGS84 lat/lon you will project once), UTC timestamps, and a segment_id so you match each contiguous sub-trip independently.
  2. Download and project the road graph. Use osmnx to fetch the drive network for the trajectory’s buffered bounding box, then project it to the trace’s metric CRS so every distance is metres.
  3. Match with HMM Viterbi. Generate candidate edges per point, score the Gaussian emission and route-distance transition, and decode the most probable edge sequence.
  4. Return and validate. Emit matched edge IDs and snapped geometry, then validate connectivity, match rate, and residual perpendicular distances before the result flows downstream.

Production-ready Python implementation

The function below is the complete wiring. It reuses extract_candidates, transition_logp, route_distance_m, and viterbi_decode from the parent reference (import them or paste them alongside), and adds the osmnx download, projection, node-mapping, and edge-case handling. Every distance is computed in a metric projected CRS — the graph is projected with osmnx.project_graph and the trace shares the same EPSG code.

PYTHON 
from __future__ import annotations

import math

import geopandas as gpd
import networkx as nx
import numpy as np
import osmnx as ox
import pandas as pd
from shapely.geometry import Point
from shapely.strtree import STRtree

# Candidate scoring + Viterbi are defined in the parent reference:
from map_matching import (  # noqa: F401
    extract_candidates,
    route_distance_m,
    transition_logp,
    viterbi_decode,
)


def match_trajectory_to_roads(
    df: pd.DataFrame,
    x_col: str = "x_m",
    y_col: str = "y_m",
    metric_epsg: int = 32633,   # UTM zone 33N — must match the cleaned trace CRS
    sigma_m: float = 8.0,
    beta_m: float = 5.0,
    search_radius_m: float = 32.0,
    bbox_buffer_m: float = 500.0,
    route_cutoff_m: float = 2000.0,
) -> pd.DataFrame:
    """
    Match a cleaned, metric-projected trajectory to an OSM road network.

    Parameters
    ----------
    df : Cleaned trajectory. Must contain x_col, y_col (metric CRS) and be
         ordered in time. A single segment_id is assumed — call once per
         segment for multi-segment trips.
    metric_epsg : EPSG of the metric CRS shared by the trace AND the graph.
                  NEVER run scoring in EPSG:4326 — degrees are not metres.
    sigma_m : GPS perpendicular-error std dev for the emission Gaussian.
    beta_m : Transition decay for route-vs-straight-line mismatch.
    search_radius_m : Max snap distance per point (~4 * sigma_m).
    bbox_buffer_m : Padding added to the trajectory bbox before download so
                    boundary edges are not clipped.
    route_cutoff_m : Upper bound on Dijkstra route search between candidates.

    Returns
    -------
    DataFrame with the original rows plus:
      matched_edge_id : (u, v, key) tuple of the matched OSM edge, or None
      snap_x, snap_y  : snapped position in the metric CRS, or NaN
      matched         : bool, False for unmatched (off-road / no candidate)

    Raises
    ------
    ValueError : missing columns, empty frame, or fewer than 2 points.
    """
    required = {x_col, y_col}
    missing = required - set(df.columns)
    if missing:
        raise ValueError(f"Missing required columns: {missing}")
    if df.empty:
        raise ValueError("Input trajectory is empty.")

    df = df.reset_index(drop=True)
    n = len(df)
    if n < 2:
        raise ValueError("Need at least 2 points to run transition scoring.")

    # ── Step 1: buffered bounding box in WGS84 for the osmnx download ──────
    # Buffer in the METRIC CRS (bbox_buffer_m is metres) so boundary edges
    # are not clipped, then project the padded extent back to lat/lon to
    # request the graph.
    pts = gpd.GeoSeries(
        [Point(xy) for xy in zip(df[x_col], df[y_col])], crs=f"EPSG:{metric_epsg}"
    )
    minx, miny, maxx, maxy = pts.total_bounds
    padded = gpd.GeoSeries(
        [Point(minx - bbox_buffer_m, miny - bbox_buffer_m),
         Point(maxx + bbox_buffer_m, maxy + bbox_buffer_m)],
        crs=f"EPSG:{metric_epsg}",
    ).to_crs("EPSG:4326")
    west, south, east, north = padded.total_bounds

    # ── Step 2: download + project the drive graph ────────────────────────
    # osmnx returns EPSG:4326; project it into the SAME metric CRS as the
    # trace so candidate distances and route lengths are both in metres.
    graph_wgs = ox.graph_from_bbox(
        bbox=(north, south, east, west),
        network_type="drive",
        truncate_by_edge=True,
    )
    graph = ox.project_graph(graph_wgs, to_crs=f"EPSG:{metric_epsg}")

    edges_gdf = ox.graph_to_gdfs(graph, nodes=False).reset_index()
    if edges_gdf.empty:
        raise ValueError("Downloaded road graph has no edges; widen the bbox.")

    # STRtree over projected edge geometries for fast radius queries.
    edge_geoms = edges_gdf.geometry.values
    edge_index = STRtree(edge_geoms)
    # Map STRtree geometry position -> (u, v, key) once.
    edge_keys = list(zip(edges_gdf["u"], edges_gdf["v"], edges_gdf["key"]))
    edges_for_extract = gpd.GeoDataFrame(
        {"geometry": edge_geoms}, crs=f"EPSG:{metric_epsg}"
    )

    # ── Step 3a: candidate generation per observation ─────────────────────
    per_point: list[list] = []
    for i in range(n):
        cands = extract_candidates(
            float(df.loc[i, x_col]),
            float(df.loc[i, y_col]),
            edges_for_extract,
            edge_index,
            sigma_m=sigma_m,
            search_radius_m=search_radius_m,
        )
        # Adaptive widen once if empty before declaring the point unmatched.
        if not cands:
            cands = extract_candidates(
                float(df.loc[i, x_col]),
                float(df.loc[i, y_col]),
                edges_for_extract,
                edge_index,
                sigma_m=sigma_m,
                search_radius_m=search_radius_m * 2.0,
            )
        per_point.append(cands)

    # ── Step 3b: build lattice, splitting where a point has no candidate ──
    # A break restarts a fresh Viterbi sub-path (handles off-road gaps and
    # disconnected graph regions without forcing spurious transitions).
    result = df.copy()
    result["matched_edge_id"] = None
    result["snap_x"] = np.nan
    result["snap_y"] = np.nan
    result["matched"] = False

    def nearest_node(x: float, y: float) -> int:
        return ox.distance.nearest_nodes(graph, x, y)

    run_start = 0
    while run_start < n:
        if not per_point[run_start]:
            run_start += 1
            continue
        run_end = run_start
        while run_end + 1 < n and per_point[run_end + 1]:
            run_end += 1
        idxs = list(range(run_start, run_end + 1))

        if len(idxs) == 1:
            # Single isolated point: take the best emission candidate.
            best = min(per_point[idxs[0]], key=lambda c: c.perp_dist_m)
            _assign(result, idxs[0], best, edge_keys)
            run_start = run_end + 1
            continue

        emission = [[c.emission_logp for c in per_point[t]] for t in idxs]
        transition: list[list[list[float]]] = []
        for a, b in zip(idxs[:-1], idxs[1:]):
            gc = math.hypot(
                df.loc[b, x_col] - df.loc[a, x_col],
                df.loc[b, y_col] - df.loc[a, y_col],
            )
            step = []
            for ci in per_point[a]:
                row = []
                node_i = nearest_node(ci.snapped_x, ci.snapped_y)
                for cj in per_point[b]:
                    node_j = nearest_node(cj.snapped_x, cj.snapped_y)
                    rdist = route_distance_m(graph, node_i, node_j, route_cutoff_m)
                    row.append(transition_logp(gc, rdist, beta_m=beta_m))
                step.append(row)
            transition.append(step)

        path = viterbi_decode(emission, transition)
        for t, chosen in zip(idxs, path):
            _assign(result, t, per_point[t][chosen], edge_keys)

        run_start = run_end + 1

    return result


def _assign(result, i, cand, edge_keys) -> None:
    """Write one matched candidate back into the result frame."""
    # cand.edge_id is the positional index shared by edge_geoms / edge_keys,
    # so it maps straight back to the (u, v, key) tuple of the OSM edge.
    result.at[i, "matched_edge_id"] = edge_keys[cand.edge_id]
    result.at[i, "snap_x"] = cand.snapped_x
    result.at[i, "snap_y"] = cand.snapped_y
    result.at[i, "matched"] = True

The lattice-splitting loop is the part that makes this robust on real data. When a point has no candidate — an off-road excursion, or a hole where the OSM extract omitted a service road — the loop closes the current Viterbi run and starts a fresh one after the gap, so a single bad point never drags its neighbours onto the wrong edge. A disconnected graph is handled the same way: route_distance_m returns +inf, transition_logp returns -inf, and the decoder simply avoids that transition.

Validation block

Run these checks immediately after matching. Each targets a distinct failure — silent CRS mismatch, low match rate, and broken connectivity — so a passing set gives real confidence before the matched path feeds downstream aggregation.

PYTHON 
def validate_match(result: pd.DataFrame, metric_epsg: int = 32633) -> dict:
    """
    Sanity-check a matched trajectory. Raises AssertionError on critical
    failures; returns a metrics dict on success.
    """
    assert not result.empty, "Matched result is empty."
    assert {"matched", "snap_x", "snap_y", "matched_edge_id"} <= set(result.columns), \
        "Missing matcher output columns."

    match_rate = result["matched"].mean()
    # A healthy urban vehicle match rate is typically > 0.9. A very low rate
    # almost always means a CRS mismatch between trace and graph, or a
    # search_radius_m that is too tight for the data's sigma.
    assert match_rate > 0.5, f"Match rate {match_rate:.2f} too low — check CRS / radius."

    # Snapped points must sit within a few sigma of their observations.
    matched = result[result["matched"]]
    if not matched.empty and {"x_m", "y_m"}.issubset(matched.columns):
        resid = np.hypot(matched["snap_x"] - matched["x_m"],
                         matched["snap_y"] - matched["y_m"])
        p95 = float(np.nanpercentile(resid, 95))
        assert p95 < 60.0, f"95th-pct snap residual {p95:.1f} m — likely wrong CRS."
    else:
        p95 = float("nan")

    # Consecutive matched edges should reference real (u, v, key) tuples.
    edge_ids = matched["matched_edge_id"].dropna().tolist()
    assert all(isinstance(e, tuple) and len(e) == 3 for e in edge_ids), \
        "Matched edge IDs are not (u, v, key) tuples."

    return {
        "n_points": len(result),
        "match_rate": round(float(match_rate), 3),
        "p95_snap_residual_m": round(p95, 1),
        "n_unique_edges": len(set(edge_ids)),
        "crs_epsg": metric_epsg,
    }

A match rate above 0.9 with a 95th-percentile snap residual under about 20 m is typical for clean urban vehicle data at 1 Hz. A sudden collapse in either metric is the classic symptom of the trace and the graph living in different coordinate systems.

Common mistakes and gotchas

  • Leaving the graph in EPSG:4326. osmnx.graph_from_bbox returns WGS84 because OSM stores lat/lon. If you skip osmnx.project_graph, every perpendicular distance and route length is computed in degrees, so the emission Gaussian and the transition decay both distort with latitude. Always project the graph to the same metric CRS as the trace before scoring.
  • A tight download bounding box. Fetching only the exact trajectory extent clips edges at the boundary, so a point near the edge of the box has no candidate and appears “off-road.” Pad the bbox by at least a few hundred metres and use truncate_by_edge=True.
  • Forcing a match on every point. An observation in a parking lot, on a ferry, or on an unmapped drive has no correct edge. Snapping it to the nearest road injects a false transition that pulls neighbouring points off their true path. Break the Viterbi lattice at unmatched points instead.
  • Unbounded Dijkstra on sparse data. With 30–60 s taxi data, an uncapped shortest-path search between distant candidates is slow and occasionally spurious. Pass a route_cutoff_m and raise beta_m so the transition term tolerates the larger legitimate gaps between observations.
  • Feeding a raw, unsmoothed trace. Map-matching assumes clean input. Residual jitter forces a large σ, which flattens the emission term until the matcher can no longer separate parallel roads. Run the trajectory through drift correction first, and derive transport mode from speed and acceleration profiling to pick the right σ and β.

Back to Map-Matching to Road Networks