Road networks: matching, routing and route flows

Question. How far did people actually travel on the road network rather than in a straight line, which roads did a GPS trace use, and how do I attach points of interest and events to trajectories?

You will get a road and a rail network for part of Beijing, a GPS trace matched to roads, road-distance versions of jump length and radius of gyration compared with the straight-line ones, flows assigned to road segments (route flows), and joins of a trajectory to nearby points of interest and events.

Prerequisites. Process trajectories. This notebook downloads road data from Amazon S3 (see Data).

Data. (1) GeoLife: the first 10 trips of user 000 (about 2 MB of the 313 MB archive) with fastmob.data.load_geolife; timestamps UTC. (2) Road and rail segments from the Overture Maps transportation theme, a bounding box of about 3 km by 3 km (roads) and 20 km by 20 km (rail) around Haidian, Beijing, fetched with fetch_road_network and fetch_rail_network from the public Overture bucket on Amazon S3. Overture removes old releases, so the notebook asks the bucket for the newest release and prints it; the amount downloaded is not reported by the library and was not measured. Overture data are open data with their own licence terms; see the Overture Maps documentation before reusing them. (3) The points of interest and events used in the last section are synthetic. The notebook stops if a download fails.

%pip install -q 'fastmob>=0.3.4' pandas matplotlib duckdb

Trajectory and road network

We take user 000's trips inside a small box, and fetch the road network of the same box. fetch_road_network returns two Arrow tables: graph nodes (road connectors with coordinates) and edges (road segments with length in metres, speed and a weight in deciseconds for time-optimal routing). We save them to disk and read them back with build_road_graph, which loads a saved graph, or fetches and saves one if the files do not exist; this is how you would cache a network you analyse repeatedly. RoadNetwork.build then prepares the routing structure once; reuse it for every query.

import re
import tempfile
import urllib.request
from pathlib import Path

import duckdb
import fastmob
import numpy as np
import pandas as pd
import pyarrow as pa
import pyarrow.parquet as pq
from fastmob.data import load_geolife
from fastmob.network import (
    RoadNetwork,
    build_rail_graph,
    build_road_graph,
    fetch_rail_network,
    fetch_road_network,
    match_trajectory,
    od_desire_lines,
    snap_locations_to_graph,
)

# Overture keeps only recent releases: ask the bucket for the newest one.
listing = (
    urllib.request.urlopen(
        "https://overturemaps-us-west-2.s3.amazonaws.com/?prefix=release/&delimiter=/"
    )
    .read()
    .decode()
)
release = max(re.findall(r"release/([0-9][0-9.\-]*)/", listing))
print("Overture release used:", release)
duckdb.sql("SET enable_progress_bar = false")

trips = load_geolife(users=["000"], max_trips_per_user=10).to_pandas()
nodes, edges = fetch_road_network(
    116.30, 39.98, 116.34, 40.01, overture_release=release
)

folder = Path(tempfile.mkdtemp())
pq.write_table(nodes, folder / "nodes.parquet")
pq.write_table(edges, folder / "edges.parquet")
nodes, edges = build_road_graph(
    116.30,
    39.98,
    116.34,
    40.01,
    overture_release=release,
    nodes_output=str(folder / "nodes.parquet"),
    edges_output=str(folder / "edges.parquet"),
)
network = RoadNetwork.build(edges, nodes)

area = trips[
    trips["lat"].between(39.985, 40.005) & trips["lng"].between(116.305, 116.335)
].reset_index(drop=True)
print(
    f"{nodes.num_rows} road nodes, {edges.num_rows} road segments; {len(area)} of {len(trips)} fixes inside the box"
)

Overture release used: 2026-10-21.0


2955 road nodes, 6365 road segments; 1184 of 6137 fixes inside the box

The box holds 2,955 road nodes and 6,365 road segments, and 1,184 of the 6,137 fixes of the first ten trips fall inside it. The notebook prints the Overture release it used, which is part of the data provenance.

Snap points to the network

snap_locations_to_graph finds the nearest graph node of each point (rows with lat and lng), up to max_distance_m metres; farther points get -1. Use it to attach zones, stops or facilities to the network.

places = pd.DataFrame(
    {"lat": [39.99, 40.00, 39.995, 39.9], "lng": [116.31, 116.32, 116.33, 116.0]}
)
snapped = snap_locations_to_graph(places, nodes, max_distance_m=500).to_pylist()
snapped

[1488, 1742, 1419, -1]

Three of the four points snap to a node within 500 m (node indexes 1,488, 1,742 and 1,419); the fourth, far outside the road box, returns -1.

Distances and routes on the network

RoadNetwork.batch_distances returns the length in metres of the fastest route between pairs of nodes (the route with the least travel time, which is not always the shortest one) (and a flag saying whether a route exists); batch_routes returns the route itself as a table of points with the accumulated travel time in deciseconds. Both take arrays of node indexes.

node_ids = nodes.column("node_idx").to_numpy()
origins, destinations = pa.array(node_ids[:5]), pa.array(node_ids[10:15])

distance_m, reachable = network.batch_distances(origins, destinations)
routes = network.batch_routes(origins, destinations)

print(np.round(distance_m.to_numpy(zero_copy_only=False)), reachable.to_pylist())
print(
    routes.num_rows,
    "route points for",
    len(set(routes.column("query_id").to_pylist())),
    "routes",
)

[ 942. 3945. 1248.    0. 3041.] [True, True, True, False, True]
143 route points for 4 routes

Between the five pairs of nodes the road distances along the fastest routes are 942, 3,945, 1,248 and 3,041 m; the fourth pair (False) has no route, so its distance is 0 and must be ignored. The routes table has 143 points for four routes: the pair without a route is absent.

Match a GPS trace to roads

match_trajectory assigns each GPS fix to the most likely road segment (a hidden-Markov map matcher). The result has the matched edge, the projected position on the road, the distance from the fix to the road, a status (matched, gap_start, or unmatched when no road was within candidate_radius_m) and a score.

matched = match_trajectory(area.head(400), network=network).to_pandas()

print(matched["match_status"].value_counts().to_dict())
print(
    f"median distance from fix to matched road: {matched['distance_to_road_m'].median():.0f} m"
)
matched[
    ["lat", "lng", "matched_lat", "matched_lng", "distance_to_road_m", "match_status"]
].head()

{'matched': 369, 'gap_start': 23, 'unmatched': 8}
median distance from fix to matched road: 20 m

Of 400 fixes, 369 are matched, 23 start a new segment after a gap (gap_start) and 8 could not be matched. The median fix lies 20 m from its matched road. That distance mixes GPS error with the offset between the road centreline and where people actually walk or drive, the width of the road and roads missing from the map, so read it as a quality signal, not as the GPS error. The projected coordinates (matched_lat, matched_lng) are on the road centreline.

Road distance versus straight line

jump_lengths_road and radius_of_gyration_road are the network versions of two earlier measures: the distance between consecutive fixes is measured along the road network instead of the great-circle line. To see a difference we keep every 20th fix (GeoLife records every 1 to 5 seconds, so about 20 to 100 seconds apart), so that consecutive fixes are far enough apart.

import fastmob.measures.individual as measures

sparse = area.iloc[::20].reset_index(drop=True)
road_jumps = measures.jump_lengths_road(sparse, network=network)
line_jumps = fastmob.jump_lengths(sparse, merge=True)
road_radius = measures.radius_of_gyration_road(sparse, network=network)
line_radius = fastmob.radius_of_gyration(sparse)

moving = np.asarray(line_jumps) > 0.05
same = np.isclose(np.asarray(road_jumps), np.asarray(line_jumps), rtol=0, atol=1e-6)
print(
    f"median road/straight-line ratio for jumps over 50 m: {np.median(np.asarray(road_jumps)[moving] / np.asarray(line_jumps)[moving]):.2f}; {same[moving].mean():.0%} of these jumps have identical road and straight-line length"
)
pd.concat(
    [
        road_radius.rename(columns={"radius_of_gyration": "road_km"}),
        line_radius.rename(columns={"radius_of_gyration": "straight_line_km"})[
            "straight_line_km"
        ],
    ],
    axis=1,
)

median road/straight-line ratio for jumps over 50 m: 1.06; 16% of these jumps have identical road and straight-line length

For jumps longer than 50 m the road distance is a median 1.06 times the straight-line distance, a small detour. But 16% of those jumps have exactly the same road and straight-line length: that is how a fallback looks, because the measure uses the straight-line distance when an end is not snapped to the network or the two ends are not connected (as for 4 of the 7 pairs in the route-flow section below). The true ratio is therefore somewhat above 1.06. The radius of gyration of the sparse trace is 1.48 km by road and 0.97 km in a straight line: the road version is larger because distances to the centre are travel distances, not geometric ones.

Accumulate flows on roads

od_desire_lines takes pairs of origin and destination nodes with a flow each, walks the fastest route of every pair and adds the flow to every road segment it crosses, so many overlapping trips pile up on shared segments. (In transport planning a desire line is the straight line between an origin and a destination; this function assigns the flows to the road network instead, which planners call an overline or route-network assignment.) It returns the segments with their accumulated flow and the dropped flow: the total flow of pairs that could not be routed because an end was not snapped or the two ends lie in disconnected parts of the graph. We use the snapped first and last fix of each recording day of user 000 as flows of one trip each.

ends = trips.groupby(trips["datetime"].dt.date).agg(
    first_lat=("lat", "first"),
    first_lng=("lng", "first"),
    last_lat=("lat", "last"),
    last_lng=("lng", "last"),
)
starts = snap_locations_to_graph(
    ends.rename(columns={"first_lat": "lat", "first_lng": "lng"}),
    nodes,
    max_distance_m=2000,
).to_numpy(zero_copy_only=False)
stops = snap_locations_to_graph(
    ends.rename(columns={"last_lat": "lat", "last_lng": "lng"}),
    nodes,
    max_distance_m=2000,
).to_numpy(zero_copy_only=False)
ok = (starts >= 0) & (stops >= 0) & (starts != stops)

lines, dropped = od_desire_lines(
    network, pa.array(starts[ok]), pa.array(stops[ok]), pa.array(np.ones(ok.sum()))
)

print(
    len(ends),
    "recording days;",
    ok.sum(),
    "with two different snapped ends;",
    lines.num_rows,
    "road segments carry flow; dropped flow",
    dropped,
    "of",
    int(ok.sum()),
)
lines.to_pandas().head(3)

9 recording days; 7 with two different snapped ends; 76 road segments carry flow; dropped flow 4.0 of 7

The nine recording days give seven pairs with two different snapped ends. Three of them could be routed, and 76 road segments carry flow (up to 2 each); the dropped flow of 4 is made of pairs whose ends could not be routed to each other, because they lie in disconnected parts of this small road box. Always check the dropped flow: a small box cuts the network, and flow that cannot be routed vanishes from the map.

Rail network

The same functions exist for rail. fetch_rail_network returns a bidirectional rail graph (nodes and edges) from Overture rail segments, with a speed for each rail class, and build_rail_graph loads a saved rail graph or fetches and saves one. We use a larger box (about 20 km by 20 km), as there are few rail lines.

rail_nodes, rail_edges = fetch_rail_network(
    116.20, 39.85, 116.45, 40.05, overture_release=release
)

pq.write_table(rail_nodes, folder / "rail_nodes.parquet")
pq.write_table(rail_edges, folder / "rail_edges.parquet")
rail_nodes, rail_edges = build_rail_graph(
    116.20,
    39.85,
    116.45,
    40.05,
    overture_release=release,
    nodes_output=str(folder / "rail_nodes.parquet"),
    edges_output=str(folder / "rail_edges.parquet"),
)

print(
    rail_nodes.num_rows,
    "rail nodes,",
    rail_edges.num_rows,
    "rail segments; classes:",
    sorted(set(rail_edges.column("class").to_pylist())),
)

2157 rail nodes, 4904 rail segments; classes: ['light_rail', 'standard_gauge', 'subway', 'tram']

The rail graph has 2,157 nodes and 4,904 segments of four classes: light rail, standard gauge, subway and tram, with a default speed per class. The rail and road graphs are separate networks.

Join trajectories to points of interest and events

join_with_pois attaches to every fix its single nearest point of interest and the distance in metres (dist_poi); join_with_pois_by_category does so for the nearest one of each category; join_with_events attaches the nearest event that happened within time_window_s seconds of the fix. The points of interest and events below are synthetic (random places on the trace, random categories, two invented events), to show the output.

from fastmob.integration import (
    join_with_events,
    join_with_pois,
    join_with_pois_by_category,
)

sample = area.sample(5, random_state=1)
pois = pd.DataFrame(
    {
        "id": range(5),
        "name_poi": list("abcde"),
        "lat": sample["lat"].to_numpy(),
        "lng": sample["lng"].to_numpy(),
        "type_poi": ["cafe", "shop", "cafe", "park", "shop"],
    }
)
events = pd.DataFrame(
    {
        "event_id": [1, 2],
        "event_type": ["concert", "match"],
        "lat": area["lat"].head(2).to_numpy(),
        "lng": area["lng"].head(2).to_numpy(),
        "datetime": area["datetime"].head(2).to_numpy(),
    }
)

nearest = join_with_pois(area.head(3), pois)
by_category = join_with_pois_by_category(area.head(3), pois)
with_events = join_with_events(area.head(5), events, time_window_s=900)

print(nearest[["id_poi", "name_poi", "dist_poi"]].round(1).to_string(index=False))
print(by_category.columns.tolist())
print(with_events[["event_id", "event_type", "dist_event"]].to_string(index=False))

 id_poi name_poi  dist_poi
      4        e    1858.2
      4        e    1844.8
      4        e    1831.8
['uid', 'lat', 'lng', 'datetime', 'id_cafe', 'dist_cafe', 'id_park', 'dist_park', 'id_shop', 'dist_shop']
event_id event_type  dist_event
       1    concert    0.000000
       2      match    0.000000
       2      match   13.210709
       2      match   26.613978
       2      match   41.540837

The nearest synthetic point of interest to the first fixes is e, about 1.8 km away. join_with_pois_by_category returns one id and distance per category (cafe, park, shop). Event 2 is attached to four of the five fixes, at distances up to 42 m, because the fixes are within 15 minutes of it; event 1 matches only the first fix.

Interpretation and limits

• Network distances are only as good as the network. Overture roads are a snapshot of the latest release, while the GeoLife data are from 2008: roads built since then appear as a straight-through route that did not exist.
• Map matching is probabilistic; check the status and score columns and inspect the result on curved roads, where the projection onto the segment is approximate.
• distance_to_road_m is a quality signal: a median of tens of metres means the fixes and the roads agree well enough for segment-level analysis, but not for lane-level analysis.
• Routing uses the fastest car route. Walkers, cyclists and trains need other networks (the rail graph is separate and not connected to the road graph).
• The synthetic points of interest and events say nothing about Beijing; replace them with real data.
• Fetching a network downloads data from a public bucket and takes time. Save the tables and reuse them.

Next steps

• Use the stops and trips of the same trajectories: Build the trajectory hierarchy.
• Compare a model's flows with observed ones: Flows between places builds origin-destination flows, and Fit and compare distributions scores them.

Related API

fetch_road_network · build_road_graph · fetch_rail_network · build_rail_graph · snap_locations_to_graph · build · batch_distances · batch_routes · match_trajectory · od_desire_lines · jump_lengths_road · radius_of_gyration_road · join_with_pois · join_with_pois_by_category · join_with_events