Multi-Modal Route Matching for Mixed Fleets

Modern logistics operations rarely rely on a single vehicle class. Mixed fleets combine heavy trucks, light commercial vans, electric cargo bikes, and passenger shuttles, each operating under distinct physical constraints, regulatory restrictions, and routing preferences. Applying a single road graph to all of them causes systematic snapping errors: an HGV gets pinned to a residential street it cannot legally enter, a cargo bike gets routed onto a motorway, and the resulting matched traces corrupt every downstream metric — from stop detection timestamps to fuel-consumption KPIs.

The core problem is that naive approaches assume homogeneous mobility. A production-grade pipeline must dynamically adapt graph topology, transition probabilities, and spatial tolerances based on real-time vehicle metadata and inferred travel modes. This page walks through exactly that, from raw telemetry through to a validated GeoJSON output, building on the Hidden Markov Model map matching pattern and the trajectory analysis framework it sits within.

Prerequisites

Before implementing the matching pipeline, ensure your environment meets the following requirements.

Python & libraries:

• Python 3.10 or later
• geopandas >= 0.14, osmnx >= 1.9, networkx >= 3.2, numpy >= 1.26, shapely >= 2.0, pandas >= 2.1
• pyarrow >= 15 for columnar telemetry storage

Data requirements:

• GPS records joined with a vehicle metadata schema that carries at minimum: vehicle_id, vehicle_class (hgv | lcv | bicycle), max_weight_kg, propulsion (ice | ev | human), and regulatory_tags (an array of OSM access-tag values).
• OSM road network extracts for your operational geography, pre-fetched via osmnx or GeoFabrik and refreshed at least monthly.
• All spatial data projected to a metric CRS — EPSG:32633 (UTM zone 33N) or EPSG:27700 (British National Grid) — from the point of ingestion. See coordinate reference system mapping for projection patterns.

Background knowledge:

• Familiarity with graph traversal and spatial indexing.
• A working understanding of the HMM emission/transition model. The HMM map matching page covers the probabilistic foundations in depth; this page extends that model with per-mode probability matrices and graph-level constraints.

Step-by-Step Workflow

Step 1 — Telemetry Ingestion & Signal Conditioning

Raw GPS streams from mixed-fleet telematics hubs contain multipath noise, clock drift, and stationary pings that degrade matching accuracy regardless of graph quality. Parse logs into a structured DataFrame, discard points where accuracy_radius > 15 metres or satellite_count < 4, then apply a Savitzky-Golay filter or lightweight Kalman smoother to suppress high-frequency jitter without distorting legitimate sharp turns.

Project coordinates to your working metric CRS immediately after cleaning. Compute instantaneous speed and bearing using vectorized operations so they are available for mode classification in the next step:

import numpy as np
import pandas as pd
from geopandas import GeoDataFrame


def compute_kinematics(gdf: GeoDataFrame) -> GeoDataFrame:
    """
    Derive speed (m/s) and heading (degrees, 0=East) from projected geometry.

    Parameters
    ----------
    gdf : GeoDataFrame
        Must carry a 'timestamp' column (datetime64[ns, UTC]) and geometry
        projected to a metric CRS so that .x / .y are in metres.

    Returns
    -------
    GeoDataFrame with added columns: dt, speed, heading.
    Rows where dt <= 0 or speed is NaN are dropped.
    """
    gdf = gdf.sort_values("timestamp").copy()
    gdf["dt"] = gdf["timestamp"].diff().dt.total_seconds()
    gdf["dx"] = gdf.geometry.x.diff()
    gdf["dy"] = gdf.geometry.y.diff()

    # Suppress division-by-zero for zero-second intervals
    valid_dt = gdf["dt"].replace(0, np.nan)
    gdf["speed"] = np.sqrt(gdf["dx"] ** 2 + gdf["dy"] ** 2) / valid_dt

    # arctan2(dy, dx): East = 0°, North = 90°
    gdf["heading"] = np.degrees(np.arctan2(gdf["dy"], gdf["dx"])) % 360

    return gdf.dropna(subset=["dt", "speed"])

Expected output shape: same row count as input minus leading NaN rows and filtered outliers; new columns dt (float64, seconds), speed (float64, m/s), heading (float64, degrees).

For deriving reliable velocity metrics from noisy pings, the speed profiling methodology covers haversine-based distance computation and outlier suppression in detail.

Step 2 — Mode Classification & Metadata Fusion

Vehicle metadata provides a strong prior, but observed kinematics frequently diverge. A delivery van crawling through a pedestrianized zone exhibits bicycle-like acceleration profiles. A cargo e-bike descending a steep grade may temporarily exceed the speed threshold associated with motorways. Relying on metadata alone produces misclassified segments; relying on kinematics alone cannot distinguish an LCV from a passenger car.

Implement a rule-based classifier that fuses both signals. Store the inferred mode as a categorical column; every downstream component reads from this label rather than raw metadata:

def classify_mode(row: pd.Series) -> str:
    """
    Fuse vehicle metadata with observed kinematics to assign a travel mode.

    Decision logic (in priority order):
      1. Metadata declares hgv → confirm with low turning rate (heading_delta < 15°).
      2. Metadata declares bicycle OR speed < 8 m/s with high turn rate → bicycle.
      3. Everything else → lcv (default motor-vehicle fallback).

    Parameters
    ----------
    row : pd.Series
        Must contain 'vehicle_class' (str), 'speed' (m/s), 'heading_delta' (degrees).

    Returns
    -------
    str: one of 'hgv', 'lcv', 'bicycle'.
    """
    declared = str(row.get("vehicle_class", "")).lower()
    speed = float(row.get("speed", 0.0))
    turn = float(row.get("heading_delta", 0.0))

    if declared == "hgv" and turn < 15.0:
        return "hgv"
    if declared == "bicycle" or (speed < 8.0 and turn > 30.0):
        return "bicycle"
    return "lcv"

Apply the classifier row-wise after computing a heading_delta column (absolute difference between consecutive headings, capped at 180°). For segment-level mode stability, apply a rolling mode (window = 5 pings) to prevent single-point anomalies from triggering graph switches.

Step 3 — Dynamic Graph Filtering & Topology Adaptation

A static road network cannot serve mixed fleets efficiently. You must generate mode-specific subgraphs by applying OSM tag filters and physical constraints against the base MultiDiGraph returned by osmnx. In osmnx, each edge carries OSM attributes as dictionary keys — the highway value may be a plain string or a list when OSM encodes multiple tags.

import osmnx as ox
import networkx as nx


def build_mode_graph(G: nx.MultiDiGraph, mode: str) -> nx.MultiDiGraph:
    """
    Return a copy of G restricted to edges accessible to `mode`.

    Parameters
    ----------
    G    : nx.MultiDiGraph returned by osmnx.graph_from_bbox / graph_from_place.
    mode : str — 'hgv', 'lcv', or 'bicycle'.

    Returns
    -------
    nx.MultiDiGraph — a new graph containing only permitted edges.
    Raises ValueError for unknown modes.

    Notes
    -----
    - Edges where access='private' are always excluded.
    - highway can be a list (OSM multi-value); set-intersection handles both cases.
    - Call .copy() so callers can mutate the result independently.
    """
    allowed: dict[str, set[str]] = {
        "hgv": {
            "motorway", "trunk", "primary", "secondary",
            "motorway_link", "trunk_link", "primary_link", "secondary_link",
        },
        "bicycle": {
            "cycleway", "residential", "living_street", "tertiary",
            "unclassified", "path", "footway",
        },
        "lcv": {
            "motorway", "trunk", "primary", "secondary", "tertiary",
            "residential", "unclassified", "living_street",
            "motorway_link", "trunk_link", "primary_link",
        },
    }

    if mode not in allowed:
        raise ValueError(f"Unknown mode '{mode}'. Expected one of {set(allowed)}")

    permitted = allowed[mode]
    keep: list[tuple] = []

    for u, v, k, data in G.edges(keys=True, data=True):
        if data.get("access") == "private":
            continue
        hw = data.get("highway", "")
        hw_set = set(hw) if isinstance(hw, list) else {hw}
        if hw_set & permitted:
            keep.append((u, v, k))

    return G.edge_subgraph(keep).copy()

Caching is mandatory in production. Serialize filtered graphs to GraphML or pickle after the first build; load them at service startup. Runtime OSM queries during high-throughput matching degrade latency from milliseconds to seconds.

For routing engine configuration notes, see the OSM Wiki routing guidelines for tag semantics, and GraphHopper’s custom model documentation for declarative per-profile constraints that mirror this filtering logic at the engine level.

Step 4 — Probabilistic Spatial Alignment

With a mode-filtered graph and cleaned trajectory, execute the HMM matcher. States represent candidate road segments; observations represent GPS coordinates. The key departure from a homogeneous HMM is that both the emission model and the transition model are parameterized per mode.

Emission probability

The standard Gaussian emission uses perpendicular distance from a GPS observation to each candidate edge. Add a cosine-similarity penalty between GPS bearing and road segment orientation — this substantially reduces false matches at intersections where several segments are equidistant:

p_emission(obs | edge) = N(d_perp; 0, sigma_z) * cos_similarity(bearing_gps, bearing_edge)

Mode-specific sigma_z values (one standard deviation of the perpendicular error distribution):

Transition probability

Transition probability encodes how likely it is to move between two candidate edges given the distance the vehicle actually travelled between observations. Use the great-circle distance as the observation, and the shortest network distance in the mode-filtered graph as the expected travel distance:

p_transition(e_t | e_{t-1}) ∝ exp(-|d_network - d_gc| / beta)

Mode-specific beta (smoothness factor for the distance-difference distribution):

Run all probability calculations in log-space to avoid underflow across long traces. A trajectory with 2,000 pings will produce products of probabilities that collapse to zero in float64 if you use linear space. See the HMM map matching page for the full Viterbi implementation in log-space.

Step 5 — Post-Processing & Trajectory Reconstruction

After Viterbi decoding, reconstruct the matched path by concatenating selected edges. Validate topological connectivity: if the algorithm jumps between disconnected components due to GPS dropouts, insert interpolated waypoints using the shortest feasible path in the mode-filtered graph rather than a straight-line segment.

Output the result as a GeoJSON FeatureCollection:

• GPS point features: original coordinates with match_status property (matched, interpolated, dropped).
• Snapped trajectory linestring: the reconstructed path geometry projected back to WGS84 (EPSG:4326) for GeoJSON compliance.
• Edge metadata per segment: osm_id, highway class, maxspeed, travel_time_s, mode.

Before stop detection can be applied, validate that the matched output is in WGS84 and that timestamps are UTC — mixed-timezone traces cause systematic dwell-time errors of hours, not seconds. The timestamp synchronization guide covers OBD-II and mobile device clock alignment in detail.

Probability & Math Model Notes

The full model is a first-order HMM. Let $O = {o_1, \ldots, o_T}$ be GPS observations and $S = {s_1, \ldots, s_T}$ be the hidden state sequence (road segment assignments). The Viterbi path maximises:

$$\hat{S} = \arg\max_S \prod_{t=1}^{T} p(o_t | s_t) \cdot p(s_t | s_{t-1})$$

Log-space rewrite (required for production):

$$\log \hat{S} = \arg\max_S \sum_{t=1}^{T} \left[\log p(o_t | s_t) + \log p(s_t | s_{t-1})\right]$$

Numerical stability notes:

• Pre-compute log-emission and log-transition matrices once per trajectory segment; avoid calling np.log inside the inner Viterbi loop.
• Clip d_perp values to a minimum of 1e-6 metres before feeding into the Gaussian to avoid log(0).
• Normalize candidate-edge emission scores per observation step so that the magnitudes across modes remain comparable when you switch graphs mid-trajectory (e.g., a driver handing off from an HGV leg to a bicycle leg).

Routing Engine Integration Notes

Mixed-fleet matching is best run against a local routing engine rather than directly against a raw networkx graph for production-scale deployments.

Valhalla exposes a costing_options object per vehicle profile that maps cleanly onto the HGV/LCV/bicycle modes above. Set use_highways, use_trails, and exclude_unpaved per mode. Valhalla’s trace_route API accepts a GPS trace and returns a matched path, offloading Viterbi decoding to the engine itself — useful if you want to benchmark against your own HMM implementation.

GraphHopper supports custom models via a JSON DSL where you declare edge-level penalties and exclusions. A custom model for HGV might penalize road_class=residential with a large multiplier rather than a hard exclusion, which produces softer transitions at the boundary of permitted networks.

Coordinate order gotcha: both Valhalla and GraphHopper accept [longitude, latitude] in their JSON payloads, consistent with GeoJSON. osmnx returns geometries in the same order. However, many Python spatial libraries (including pyproj.Transformer) default to (latitude, longitude) unless always_xy=True is passed. Verify coordinate order at every API boundary — transposed coordinates produce matching errors of hundreds of kilometres that are extremely difficult to debug from match-rate metrics alone.

Rate-limit mitigations for local engines: spin up a dedicated Valhalla or GraphHopper instance per region. For bursty batch workloads, use a semaphore-limited thread pool (max 4 concurrent requests per engine instance) to prevent request queue saturation and OOM kills.

Operational Troubleshooting

Deployment Checklist

Project all spatial inputs to a metric CRS (EPSG:32633 or equivalent) before the first distance calculation.

Verify that the OSM graph is projected to the same CRS as the trajectory (G.graph['crs'] matches gdf.crs).

Run Viterbi scoring in log-space; confirm no -inf values propagate through the backpointer matrix.

Cache mode-filtered subgraphs to disk; validate on startup by checking node/edge counts against expected ranges.

Implement GPS dropout detection (accuracy threshold + satellite count) and dead-reckoning fallback.

Add a mode-mismatch detector that auto-reclassifies and emits mode_override warnings.

Segment traces longer than 15 minutes into overlapping windows before running the HMM.

Validate coordinate order ([lon, lat]) at every routing engine API boundary.

Confirm output GeoJSON is in WGS84 before handing off to stop detection and downstream analytics.

Set latency SLAs per mode: HGV batch (< 200 ms per trajectory-km), real-time LCV streaming (< 50 ms per ping window), bicycle (< 30 ms per ping window).

Scaling & Production Deployment

Spatial partitioning: Divide operational geography into H3 hexagonal tiles (resolution 6–7) or quadtree cells. Process trajectories within each tile independently and merge at tile boundaries. This prevents a single large trajectory from monopolizing a worker and enables straightforward horizontal scaling.

Graph caching: Serialize mode-filtered graphs to GraphML or pickle after the first build. Load them into memory at service startup and never issue runtime OSM queries during matching. Invalidate caches on a weekly schedule when OSM extracts are refreshed.

Batch vs. stream: Use Apache Kafka for real-time telemetry ingestion and run matching in micro-batches (5-minute windows) to balance latency with computational efficiency. For overnight batch jobs over historical traces, multiprocessing.Pool with chunksize=50 trajectories per worker is straightforward and avoids the serialization overhead of distributed frameworks for moderate fleet sizes.

Outlier removal at ingestion: Filter impossible-speed pings (> 55 m/s for any ground vehicle) before they reach the mode classifier. A single bad ping can force the classifier into bicycle mode for an HGV and corrupt the entire trajectory match.

Directionality & heading synchronization: For vehicles equipped with gyroscopes or accelerometers, fuse IMU heading into the emission model to improve accuracy in GPS-degraded corridors. Even a 1 Hz IMU sample reduces heading ambiguity errors by roughly 40% in dense urban grids.

Validation & Continuous Improvement

Deploy a shadow-matching pipeline alongside production. Compare algorithmic outputs against manually verified ground-truth routes for a weekly sample (200+ trajectories per mode). Track:

• Match rate: percentage of GPS points successfully snapped to a valid edge in the mode-filtered graph.
• Topological error rate: frequency of illegal turns or disconnected jumps, detected by running nx.is_weakly_connected on the matched path subgraph.
• Mode misclassification rate: instances where kinematic behaviour contradicts the assigned vehicle class (measured against manually labelled ground truth).
• Interpolation ratio: percentage of output path length that is interpolated rather than directly matched; values above 15% indicate systematic signal quality or graph connectivity problems.

Retrain sigma_z and beta parameters quarterly using maximum-likelihood estimation over the matched ground-truth sample. As OSM coverage improves and vehicle telematics hardware upgrades (4G → 5G, GPS → GNSS multi-constellation), your emission noise floor will decrease — tightening sigma_z will improve precision without sacrificing recall.

Related

• Hidden Markov Model Map Matching in Python — the probabilistic foundation this page extends with mode-specific probability matrices
• Speed Profiling from Raw GPS Coordinates — kinematic signal preparation that feeds the mode classifier
• Directionality & Heading Synchronization — IMU fusion techniques for improving heading reliability in the emission model
• Stop Detection & Dwell Time Analytics — the downstream consumer of matched trajectories; requires WGS84 output and UTC timestamps
• Coordinate Reference System Mapping for Fleet Data — CRS projection patterns referenced in the prerequisites

Parent: Trajectory Analysis & Map Matching Techniques