Building an OSM road graph with OSMnx for map matching

import networkx as nx
import osmnx as ox
import geopandas as gpd


def build_match_graph(
    bounds: tuple,
    target_crs: str,
    network_type: str = "drive_service",
    buffer_m: float = 500.0,
    simplify: bool = True,
) -> tuple:
    """
    Download, project and prune an OSM graph ready for map matching.

    Parameters
    ----------
    bounds : tuple
        (minx, miny, maxx, maxy) of the trajectory, in target_crs units.
    target_crs : str
        The trace's PROJECTED metric CRS. osmnx downloads in EPSG:4326;
        every distance in an emission or transition model must be metres,
        so the graph is projected before it is used, not after.
    network_type : str
        'drive_service' includes service roads and car parks, which delivery
        fleets use constantly and 'drive' omits. 'all' adds footways, which
        a matcher will happily snap a van onto.
    buffer_m : float
        Expansion of the download box. 500 m is a safe urban default; raise
        it for motorway or sparse data where consecutive fixes are far apart.

    Returns
    -------
    (graph, edges_gdf)

    Raises
    ------
    ValueError
        On a degenerate box, a geographic target CRS, or an empty download.
    """
    minx, miny, maxx, maxy = bounds
    if maxx <= minx or maxy <= miny:
        raise ValueError(f"Degenerate bounds {bounds}.")

    crs = gpd.GeoSeries([], crs=target_crs).crs
    if crs.is_geographic:
        raise ValueError(
            f"target_crs {target_crs} is geographic. Matching in degrees makes "
            "the emission Gaussian latitude-dependent and the transition "
            "distances meaningless."
        )

    # Buffer in metres, then convert the box to WGS84 for the download.
    box_m = gpd.GeoSeries.from_wkt(
        [f"POLYGON(({minx-buffer_m} {miny-buffer_m}, {maxx+buffer_m} {miny-buffer_m}, "
         f"{maxx+buffer_m} {maxy+buffer_m}, {minx-buffer_m} {maxy+buffer_m}, "
         f"{minx-buffer_m} {miny-buffer_m}))"],
        crs=target_crs,
    )
    west, south, east, north = box_m.to_crs(4326).total_bounds

    g = ox.graph_from_bbox(
        bbox=(west, south, east, north),
        network_type=network_type,
        simplify=simplify,           # keeps geometry on the edge, not the nodes
        retain_all=False,            # drop islands unreachable from the main graph
    )
    if g.number_of_edges() == 0:
        raise ValueError("Download returned an empty graph — check the bbox order.")

    g = ox.project_graph(g, to_crs=target_crs)
    edges = ox.graph_to_gdfs(g, nodes=False, edges=True)
    return g, edges


def graph_health(g, edges: gpd.GeoDataFrame, trace: gpd.GeoDataFrame,
                 search_radius_m: float = 60.0) -> dict:
    """
    Coverage and connectivity checks. Run these BEFORE matching, not after —
    a matcher on a broken graph produces confident nonsense rather than errors.
    """
    from shapely import STRtree

    # 1. Connectivity. A matcher's transition term needs a route between
    #    consecutive candidates; a fragmented graph makes that impossible.
    comps = list(nx.weakly_connected_components(g))
    largest = max(len(c) for c in comps)
    frac_connected = largest / g.number_of_nodes()

    # 2. Coverage. What share of fixes have at least one candidate edge?
    tree = STRtree(edges.geometry.values)
    have = 0
    for pt in trace.geometry.values:
        if len(tree.query(pt.buffer(search_radius_m))) > 0:
            have += 1
    coverage = have / max(len(trace), 1)

    return {
        "nodes": g.number_of_nodes(), "edges": g.number_of_edges(),
        "components": len(comps), "largest_component_frac": frac_connected,
        "fix_coverage": coverage,
    }

def validate_graph(health: dict, min_coverage: float = 0.98) -> None:
    """Three assertions that separate a usable graph from an expensive mistake."""
    # 1. Coverage: a fix with no candidate cannot be matched at all.
    assert health["fix_coverage"] >= min_coverage, (
        f"only {health['fix_coverage']:.1%} of fixes have a candidate edge — "
        "increase the buffer, or the network type excludes roads in use"
    )

    # 2. Connectivity: the transition term needs routes between candidates.
    assert health["largest_component_frac"] > 0.95, (
        f"largest component holds only {health['largest_component_frac']:.0%} "
        "of nodes — the graph is fragmented and transitions will fail"
    )

    # 3. Size sanity. A graph with very few edges for its area usually means
    #    the network type filtered out almost everything.
    assert health["edges"] > 100, (
        f"only {health['edges']} edges — check network_type and the bbox order "
        "(osmnx takes west, south, east, north)"
    )
    print(f"OK — {health['nodes']} nodes, {health['edges']} edges, "
          f"coverage {health['fix_coverage']:.1%}")