Geometric Nearest-Edge Snapping in Python

This page extends the snapper sketched in choosing a map-matching algorithm into a self-contained, copy-paste class you can drop into a preprocessing job. The specific edge case it addresses is the one that makes naive snapping usable in practice: a single nearest-edge query flickers between parallel roads and jumps across junctions, so this implementation retrieves several candidates per ping and disciplines the choice with a perpendicular-distance gate, an optional heading gate, and a one-step continuity check. It is still a stateless geometric method — it has no memory of connectivity and no probabilities — so it remains the right tool only where roads are well-separated or sampling is dense enough to keep matches unambiguous.

Compatibility and Configuration Requirements

The class assumes the road network is already loaded as a projected GeoDataFrame of LineString edges — typically fetched with osmnx and reprojected. Both the pings and the edges must be in the same metric CRS.

The Snapper Class

import numpy as np
from shapely import STRtree
from shapely.geometry import Point, LineString
import geopandas as gpd
from typing import Optional


class NearestEdgeSnapper:
    """
    Stateless geometric map matcher: snap GPS pings to the nearest road edge.

    Retrieves the k nearest candidate edges per ping and selects the best one
    that passes a perpendicular-distance gate and an optional heading gate,
    then flags matches that jump implausibly far from the previous match.

    Parameters
    ----------
    edges : geopandas.GeoDataFrame
        Road segments with a LineString geometry column, in a metric CRS.
    max_snap_distance_m : float
        Reject any candidate whose perpendicular distance to the ping exceeds
        this many metres. Set to roughly 3x the GPS error std dev: ~25-40 m for
        consumer trackers, tighter for RTK. Too small drops valid matches near
        wide roads; too large lets pings snap to unrelated nearby streets.
    candidate_k : int
        Number of nearest edges retrieved per ping. k=1 is pure nearest-edge
        (fast, but no heading disambiguation). k=4-6 lets the heading gate pick
        the correctly-oriented edge on parallel roads. Larger k costs query time.
    heading_gate_deg : float or None
        If a ping heading is supplied, reject candidates whose edge bearing
        differs from it by more than this many degrees (mod 180, so direction of
        digitisation does not matter). None disables the gate. 45 deg is a safe
        default; tighten to 30 for well-mapped urban grids.
    max_jump_m : float
        Continuity check: flag (do not drop) a match whose snapped point is more
        than this far from the previous ping's snapped point. Surfaces junction
        jumps for review. Scale with sampling interval and speed limit.
    """

    def __init__(
        self,
        edges: gpd.GeoDataFrame,
        max_snap_distance_m: float = 30.0,
        candidate_k: int = 5,
        heading_gate_deg: Optional[float] = 45.0,
        max_jump_m: float = 400.0,
    ):
        if edges.crs is None or edges.crs.is_geographic:
            raise ValueError(
                "edges must be in a projected metric CRS (e.g. a UTM zone); "
                "snapping on geographic degrees corrupts distances."
            )
        self.edges = edges.reset_index(drop=True)
        self.geoms = self.edges.geometry.values
        self.tree = STRtree(self.geoms)          # built once; reuse per trace
        self.max_snap_distance_m = max_snap_distance_m
        self.candidate_k = candidate_k
        self.heading_gate_deg = heading_gate_deg
        self.max_jump_m = max_jump_m

    @staticmethod
    def _edge_bearing_deg(edge: LineString, along: Point) -> float:
        """Local bearing of the edge near the projected point, in [0, 360)."""
        s = edge.project(along)                  # arc-length of nearest point
        a = edge.interpolate(max(s - 1.0, 0.0))
        b = edge.interpolate(min(s + 1.0, edge.length))
        ang = np.degrees(np.arctan2(b.x - a.x, b.y - a.y))
        return ang % 360.0

    @staticmethod
    def _heading_diff_mod180(a: float, b: float) -> float:
        """Smallest angular difference ignoring travel direction, in [0, 90]."""
        d = abs(a - b) % 180.0
        return min(d, 180.0 - d)

    def snap_point(
        self,
        x: float,
        y: float,
        heading_deg: Optional[float] = None,
        prev_xy: Optional[tuple] = None,
    ) -> dict:
        """
        Snap a single ping. Returns a dict with the chosen edge index, snapped
        coordinates, perpendicular distance, jump distance, and a suspect flag.
        edge_i is None when no candidate passes the gates.
        """
        pt = Point(x, y)
        # query_nearest returns candidate edge indices, nearest first
        idx = self.tree.query_nearest(pt, all_matches=True)
        idx = np.atleast_1d(idx)[: self.candidate_k]

        best = None
        for edge_i in idx:
            edge = self.geoms[int(edge_i)]
            perp_m = pt.distance(edge)
            if perp_m > self.max_snap_distance_m:
                continue
            proj = edge.interpolate(edge.project(pt))
            if heading_deg is not None and self.heading_gate_deg is not None:
                bearing = self._edge_bearing_deg(edge, proj)
                if self._heading_diff_mod180(bearing, heading_deg) > self.heading_gate_deg:
                    continue
            if best is None or perp_m < best["perp_m"]:
                best = {"edge_i": int(edge_i), "snap_x": proj.x,
                        "snap_y": proj.y, "perp_m": float(perp_m)}

        if best is None:
            return {"edge_i": None, "snap_x": None, "snap_y": None,
                    "perp_m": None, "jump_m": None, "suspect": True}

        jump_m = 0.0
        if prev_xy is not None:
            jump_m = Point(best["snap_x"], best["snap_y"]).distance(Point(prev_xy))
        best["jump_m"] = float(jump_m)
        best["suspect"] = bool(jump_m > self.max_jump_m)
        return best

    def snap_trace(self, pings: np.ndarray, headings: Optional[np.ndarray] = None) -> list:
        """
        Snap an ordered trace. pings is shape (N, 2) of [x, y] in the edge CRS;
        headings is an optional (N,) array of bearings in degrees. Returns a list
        of per-ping result dicts in input order.
        """
        results = []
        prev_xy = None
        for i in range(len(pings)):
            hd = None if headings is None else float(headings[i])
            res = self.snap_point(pings[i, 0], pings[i, 1], heading_deg=hd, prev_xy=prev_xy)
            results.append(res)
            if res["edge_i"] is not None:
                prev_xy = (res["snap_x"], res["snap_y"])
        return results

Execution and Tuning

Instantiate once against the projected road network, then call snap_trace per trip. Build the STRtree a single time and reuse the instance — reconstructing the index per trace costs roughly 200 ms per 100k edges and dwarfs the query time.

snapper = NearestEdgeSnapper(
    edges_utm,                 # GeoDataFrame in EPSG:326xx
    max_snap_distance_m=30.0,
    candidate_k=5,
    heading_gate_deg=45.0,
)
matched = snapper.snap_trace(pings_xy, headings=ping_bearings)
suspect_rate = sum(m["suspect"] for m in matched) / len(matched)

The knobs, and the effect of moving each:

• max_snap_distance_m — the hard reject threshold. Raise it when clean pings near wide multi-lane roads are being dropped (edge_i returns None); lower it when pings are snapping to unrelated streets one block over. A value near 3 · σ_z (three times the GPS error standard deviation) is the principled starting point.
• candidate_k — how many edges the heading gate gets to choose among. At k=1 the heading gate is inert because there is no alternative to fall back on. Raise to 5–6 on networks with parallel carriageways and service roads; each extra candidate adds query and projection cost.
• heading_gate_deg — how much bearing disagreement is tolerated before an edge is rejected. Tighten toward 30° on well-mapped urban grids to aggressively kill parallel-road flicker; loosen toward 60° where the heading signal is noisy or the road geometry is coarsely digitised. Supply headings from the ping stream or derive them during speed profiling; without a heading the gate is skipped entirely.
• max_jump_m — the continuity flag. It never drops a match, only marks suspect, so use it as a monitoring signal. Scale it with the sampling interval and the local speed limit: a 30-second gap at 90 km/h legitimately covers 750 m, so a flat 400 m threshold will over-flag sparse traces.

Track suspect_rate per region as your escalation trigger. When it stays high on dense-urban traces despite tuning, the network is too ambiguous for a stateless method and those traces belong on a hidden Markov model matcher.

Common Pitfalls

Up: Choosing a Map-Matching Algorithm: HMM vs Geometric Snapping vs Viterbi | Trajectory Analysis & Map Matching Techniques

Related

• Choosing a Map-Matching Algorithm: HMM vs Geometric Snapping vs Viterbi — the decision guide that places this snapper against the HMM and explains when to escalate
• Hidden Markov Model Map Matching in Python — where to go when parallel-road flicker and junction jumps prove structural
• Speed Profiling from Raw GPS Coordinates — derives the per-ping heading the snapper’s heading gate consumes
• Trajectory Analysis & Map Matching Techniques — the parent guide covering the full matching pipeline and its position in the telematics data flow