Shortest-Path Routing with NetworkX and OSMnx

Routing looks like a one-liner until the first answer comes back wrong: a route that drives the wrong way up a one-way, a distance that is 800 m short because the origin snapped to the far side of a river, or a "shortest" path that takes four minutes longer than the obvious one. This guide is for anyone turning a pair of coordinates into a defensible travel time and a drawable line. It sits under Street Network Analysis with OSMnx in Spatial Analysis & Advanced Query Techniques, and it is the routing half of the same graph the service areas in building isochrones from a street network are cut from.

Why This Approach / What Goes Wrong

A shortest path is only defined relative to a cost, and OSMnx hands you a graph with two plausible ones. length is present on every edge from the moment the graph is built — metres, computed geodesically between incident nodes. travel_time does not exist until you ask for it, and it is an imputation: add_edge_speeds fills the gaps in OpenStreetMap's patchy maxspeed tagging, and add_edge_travel_times divides metres by metres-per-second. The two costs disagree constantly, because the shortest way through a city is a zigzag of residential streets and the fastest way is a detour onto an arterial. Choosing weight="length" when the question was "how long does the ambulance take" is the single most common routing error, and it never raises — it just returns a slower route with a smaller number attached.

Three structural traps sit underneath that choice. The first is directionality: the graph is a MultiDiGraph, so a one-way street contributes exactly one directed edge, and the drive from A to B is genuinely a different problem from the drive from B to A. Any code that collapses the graph to an undirected Graph for speed silently legalises wrong-way driving. The second is connectivity. OSMnx's default keeps the largest weakly connected component, which ignores edge direction, so a block reachable only against a one-way survives in the graph and makes the router return None for every trip into it. The third is snapping: nearest_nodes matches your coordinate to the nearest intersection by straight-line distance, not to the nearest point on the nearest street, so both endpoints carry an error of up to half a block before the router has done any work.

Prerequisites

• osmnx>=2.1,<3 — graph construction plus the ox.routing and ox.distance helpers used here
• networkx>=3.3 — the algorithms themselves: Dijkstra, A*, single-source solves
• geopandas>=1.0 and shapely>=2.0 — projecting the endpoints and assembling the output geometry
• scipy>=1.11 — required by nearest_nodes for the KD-tree on a projected graph

python -m pip install "osmnx>=2.1,<3" "networkx>=3.3" "geopandas>=1.0" "shapely>=2.0" "scipy>=1.11"

Everything below assumes the OSMnx 2.x module layout, where routing helpers live under ox.routing rather than the 1.x ox.speed and ox.distance locations; the full set of 1.x-to-2.x breaks is catalogued in Street Network Analysis with OSMnx.

Step-by-Step Implementation

1. Build a graph that is actually routable, then project it.

Two preparation steps decide whether routing works at all. largest_component(..., strongly=True) guarantees every node can reach every other following edge direction, which is the precondition the router needs and the weak default does not give you. Projecting to the local UTM zone then makes node coordinates metric, which matters for the snap distances and the A* heuristic later — pick the zone with the graph's own centroid rather than hard-coding one, and never reach for Web Mercator, whose scale error grows with latitude.

import networkx as nx
import osmnx as ox

ox.settings.use_cache = True

streets = ox.graph_from_point((45.0703, 7.6869), dist=2500, network_type="drive")

# Routing needs *strong* connectivity: a block you can only leave by driving the
# wrong way up a one-way survives the default weak filter and breaks the solver.
streets = ox.truncate.largest_component(streets, strongly=True)

# Metric node coordinates — needed for snap distances and the A* heuristic
streets = ox.projection.project_graph(streets)

print(streets.graph["crs"])                              # EPSG:32632
print(len(streets), streets.number_of_edges())           # 3894 8721

2. Impute speeds, then derive travel times.

add_edge_speeds writes a speed_kph column: it parses each edge's maxspeed tag, and where the tag is missing it substitutes the mean of the observed values for that highway class in this graph. hwy_speeds overrides that per class, and fallback catches classes with no observed value anywhere. add_edge_travel_times then writes travel_time in seconds and rounds both attributes to one decimal.

streets = ox.routing.add_edge_speeds(
    streets,
    hwy_speeds={"residential": 30, "living_street": 20, "primary": 50},
    fallback=40,
)
streets = ox.routing.add_edge_travel_times(streets)

u, v, key, data = next(iter(streets.edges(keys=True, data=True)))
print(round(data["length"], 1), data["speed_kph"], data["travel_time"])
# 84.3 50.0 6.1

Run this after truncating the graph, not before: any edge added or restored later will be missing speed_kph, and add_edge_travel_times raises rather than guessing. Be honest with yourself about what the number means. OSMnx parses plain integers and mph suffixes and nothing else, so implicit country codes such as RO:urban or a walk value are dropped and the class mean silently takes over. Coverage varies enormously between cities — a well-mapped European core may tag maxspeed on 60% of drivable edges, a sparsely mapped one on under 5%, and in the second case almost every travel_time in your output is your own hwy_speeds dictionary reflected back at you. The values are free-flow limits regardless: no congestion, no signals, no turn delay, so treat the result as a lower bound on real drive time and calibrate against observed trips before publishing minutes to anyone.

3. Snap the endpoints, and look at how far they moved.

Addresses arrive in longitude/latitude. Project them with the graph's own CRS so the KD-tree query happens in the same metric space as the node coordinates. Shapely and GeoPandas are unconditionally x-then-y, so Point(lon, lat) is correct for EPSG:4326 regardless of that authority's declared lat-lon axis order — the same always_xy=True convention described in Coordinate Systems with PyProj.

import geopandas as gpd
from shapely.geometry import Point

stops = gpd.GeoSeries(
    [Point(7.6771, 45.0723), Point(7.7035, 45.0629)],   # (lon, lat) — depot, client
    crs="EPSG:4326",
).to_crs(streets.graph["crs"])

nodes, snap_dist = ox.distance.nearest_nodes(
    streets, X=stops.x.to_numpy(), Y=stops.y.to_numpy(), return_dist=True
)
orig, dest = nodes
print(nodes, [round(d, 1) for d in snap_dist])
# [1745200343, 258161573] [18.4, 41.2]

return_dist=True is not optional in production code. On a projected graph the distances come back in metres, and a value in the hundreds means the address landed across a motorway, a river or a rail cutting — the KD-tree measures straight lines, exactly like the cKDTree path benchmarked in sjoin_nearest vs cKDTree, and it has no idea whether the two points are connected.

4. Route once per weight, and total the real edge attributes.

route_len = ox.routing.shortest_path(streets, orig, dest, weight="length")
route_time = ox.routing.shortest_path(streets, orig, dest, weight="travel_time")


def route_totals(G, route, weight):
    """Sum the true edge attributes along a route.

    Pass the same weight you routed with: route_to_gdf resolves each parallel
    edge bundle by picking its minimum-weight member, so a mismatched weight
    can total a slightly different set of edges than the router traversed.
    """
    edges = ox.routing.route_to_gdf(G, route, weight=weight)
    return edges["length"].sum(), edges["travel_time"].sum()


for label, route in [("length", route_len), ("travel_time", route_time)]:
    metres, seconds = route_totals(streets, route, label)
    print(f"{label:>12}: {len(route):3d} nodes {metres:7.0f} m {seconds:6.0f} s")
# Expected output:
#       length:  61 nodes    2361 m    281 s
#  travel_time:  48 nodes    2604 m    243 s

The time-optimal route is 243 m longer and 38 seconds quicker, and it uses thirteen fewer intersections — the signature of an arterial detour. ox.routing.shortest_path is a thin wrapper over nx.shortest_path that swallows NetworkXNoPath and returns None instead, which is convenient in a batch and dangerous in a script that then indexes into the result.

5. Rebuild the node list as a LineString.

A route is a list of integers; nothing downstream can draw that. route_to_gdf gives you the ordered edges with their geometries, and the one-liner linemerge(edges.geometry.tolist()) works when the edges are contiguous — but it discards travel direction and degrades to a MultiLineString at the first gap. When the direction of travel matters (arrows on a map, animating a vehicle, splitting the trip at a midpoint), assemble the coordinates yourself and orient each segment to start at its u node.

from shapely.geometry import LineString


def route_linestring(G, route, weight="travel_time"):
    """Concatenate the route's edge geometries into one directed LineString."""
    coords = []
    for u, v in zip(route[:-1], route[1:]):
        key = min(G[u][v], key=lambda k: G[u][v][k][weight])
        data = G[u][v][key]
        ux, uy = G.nodes[u]["x"], G.nodes[u]["y"]
        if "geometry" in data:
            seg = list(data["geometry"].coords)
        else:                                   # straight edge: no stored geometry
            seg = [(ux, uy), (G.nodes[v]["x"], G.nodes[v]["y"])]
        # Edge geometry is not guaranteed to run u -> v; flip it if it does not
        if (seg[0][0] - ux) ** 2 + (seg[0][1] - uy) ** 2 > (seg[-1][0] - ux) ** 2 + (seg[-1][1] - uy) ** 2:
            seg.reverse()
        coords.extend(seg if not coords else seg[1:])
    return LineString(coords)


line = route_linestring(streets, route_time)
trip = gpd.GeoDataFrame({"trip_id": ["depot-to-client"]}, geometry=[line],
                        crs=streets.graph["crs"])
trip.to_crs("EPSG:4326").to_file("route.geojson", driver="GeoJSON")

print(f"{line.length:.0f} m planar vs {route_totals(streets, route_time, 'travel_time')[0]:.0f} m geodesic")
# 2603 m planar vs 2604 m geodesic

Those two numbers differ by the UTM scale factor — about one part in two and a half thousand. Report the summed length attribute, which OSMnx computed geodesically at build time; treat the projected geometry as a drawing, not a measurement.

6. Swap Dijkstra for A* when the graph gets large.

nx.shortest_path runs Dijkstra, which expands outward from the origin in every direction until the destination is settled. A* adds a heuristic estimate of the remaining cost and steers the frontier toward the target, but only stays correct if that heuristic never over-estimates. On a projected graph the straight-line distance is admissible for weight="length" by definition. For weight="travel_time" you must convert those metres to seconds at the graph's top speed, otherwise the estimate can exceed the true remaining time and A* will happily return a suboptimal route.

import math
import timeit

top_speed_m_s = max(d["speed_kph"] for _, _, d in streets.edges(data=True)) / 3.6


def seconds_as_the_crow_flies(node, target):
    """Admissible heuristic: no route can beat a straight line at the top speed."""
    dx = streets.nodes[node]["x"] - streets.nodes[target]["x"]
    dy = streets.nodes[node]["y"] - streets.nodes[target]["y"]
    return math.hypot(dx, dy) / top_speed_m_s


route_astar = nx.astar_path(streets, orig, dest,
                            heuristic=seconds_as_the_crow_flies, weight="travel_time")

cost = lambda r: sum(min(streets[u][v].values(), key=lambda d: d["travel_time"])["travel_time"]
                     for u, v in zip(r[:-1], r[1:]))
assert abs(cost(route_astar) - cost(route_time)) < 1e-6   # same cost, maybe different ties

runs = 20
dij = timeit.timeit(lambda: nx.shortest_path(streets, orig, dest, weight="travel_time"),
                    number=runs) / runs
ast = timeit.timeit(lambda: nx.astar_path(streets, orig, dest,
                                          heuristic=seconds_as_the_crow_flies,
                                          weight="travel_time"), number=runs) / runs
print(f"dijkstra {dij * 1000:.1f} ms | a* {ast * 1000:.1f} ms")
# dijkstra 41.3 ms | a* 18.7 ms

The heuristic is weakest exactly where the speed spread is widest: dividing by a 130 km/h motorway limit makes the estimate for a 30 km/h residential grid so loose that A* degenerates toward Dijkstra. A* earns its keep on regional graphs where the origin and destination are far apart relative to the extract; on a few thousand nodes the difference is real but rarely decisive, and it is never worth the extra code if you are about to batch the calls anyway.

7. Batch many pairs without paying per pair.

ox.routing.shortest_path accepts iterables for orig and dest and solves them in parallel, with cpus=None meaning every core. That is the right tool for unrelated pairs, but each worker process receives a full pickled copy of the graph, so a large graph times sixteen cores is a memory problem, not a speed-up.

When the pairs share an origin — one depot to every customer, one clinic to every census tract — a single-source solve answers all of them in one Dijkstra instead of n. Because the graph is directed, "depot to everyone" and "everyone to the depot" are different problems: run the second one on a reversed view, which NetworkX gives you for free without copying the graph.

customers = list(streets.nodes)[:250]          # stand-in for snapped delivery points

# Unrelated pairs: parallel, one process per core
pairs_out = ox.routing.shortest_path(streets, [orig] * 5, customers[:5],
                                     weight="travel_time", cpus=1)

# Shared origin: one solve answers every destination, with a 30-minute horizon
out_cost, out_path = nx.single_source_dijkstra(streets, orig, cutoff=1800,
                                               weight="travel_time")

# Shared destination: reverse the edge directions, not the graph in memory
inbound = streets.reverse(copy=False)          # read-only view
in_cost, in_path = nx.single_source_dijkstra(inbound, orig, cutoff=1800,
                                              weight="travel_time")

reachable = [c for c in customers if c in out_cost]
asymmetric = [c for c in reachable if c in in_cost and abs(out_cost[c] - in_cost[c]) > 30]
print(len(reachable), "reachable within 30 min;", len(asymmetric), "differ by >30 s each way")
# 231 reachable within 30 min; 47 differ by >30 s each way

Those 47 asymmetric pairs are the one-way system made visible: the same physical trip costs measurably more in one direction. Note that in_path lists run from the origin in the reversed graph, so reverse each list to recover the real drive order.

The same idea scales to a full cost matrix. A naive n×m matrix costs n×m separate solves; one solve per distinct origin costs n, and single_source_dijkstra_path_length skips building the path dictionaries you are going to discard. Add a cutoff in the units of the weight — seconds for travel_time — so the frontier stops instead of settling the whole graph for destinations you would reject anyway.

import numpy as np

depots, clients = customers[:8], customers[8:60]
matrix = np.full((len(depots), len(clients)), np.inf)

for i, source in enumerate(depots):
    reach = nx.single_source_dijkstra_path_length(streets, source, cutoff=900,
                                                  weight="travel_time")
    for j, target in enumerate(clients):
        if target in reach:
            matrix[i, j] = reach[target]

print(np.isfinite(matrix).sum(), "of", matrix.size, "pairs inside 15 minutes")
# 341 of 416 pairs inside 15 minutes

Rows that come back entirely infinite are the tell for a depot stranded behind a one-way loop or left over from a weakly connected component — a cheaper diagnostic than inspecting the graph.

Verification

Three assertions catch nearly every routing bug: that the route is a legal walk in the directed graph, that the cost you report is the cost the solver actually minimised, and that each weight really did win on its own metric.

# 1. Every consecutive pair is a real directed edge — catches wrong-way routes
assert route_time[0] == orig and route_time[-1] == dest
assert all(streets.has_edge(u, v) for u, v in zip(route_time[:-1], route_time[1:])), \
    "Route traverses an edge that does not exist in this direction"

# 2. The optimum is a lower bound on any route's cost under the same weight
best_m = nx.shortest_path_length(streets, orig, dest, weight="length")
best_s = nx.shortest_path_length(streets, orig, dest, weight="travel_time")
m_by_len, s_by_len = route_totals(streets, route_len, "length")
m_by_time, s_by_time = route_totals(streets, route_time, "travel_time")
assert best_m <= m_by_time + 1e-6 and best_s <= s_by_len + 1e-6

# 3. The endpoints did not snap into the next neighbourhood
assert max(snap_dist) < 150, f"Endpoint moved {max(snap_dist):.0f} m during snapping"

print(f"detour of {m_by_time - m_by_len:.0f} m buys {s_by_len - s_by_time:.0f} s")
# detour of 243 m buys 38 s

Edge Cases & Debugging

• shortest_path returns None, or nx.NetworkXNoPath is raised. The destination is unreachable following edge direction. Re-run ox.truncate.largest_component(streets, strongly=True) and confirm with nx.is_strongly_connected(streets).
• KeyError: 'travel_time'. add_edge_speeds and add_edge_travel_times were never run, or were run before a step that added edges. Impute after every structural change to the graph, and check with sum("travel_time" not in d for _, _, d in streets.edges(data=True)).
• The route drives the wrong way up a one-way. Something converted the graph: ox.convert.to_undirected and nx.Graph(streets) both discard direction. Route on the MultiDiGraph and only collapse to a DiGraph (via ox.convert.to_digraph) for algorithms that reject parallel edges.
• A* returns a slower route than Dijkstra. The heuristic over-estimates. On an unprojected graph the node coordinates are degrees and math.hypot is meaningless; on travel_time you must divide by the graph's maximum speed, not an average.
• route_to_gdf raises on a one-node route. orig == dest because both coordinates snapped to the same intersection. Guard with if orig == dest and return a zero-cost trip before calling it.
• Parallel batching exhausts RAM. Each cpus worker pickles the whole graph. Drop to cpus=1 and restructure as single-source solves, which reuse one in-process graph and are usually faster anyway.

Frequently Asked Questions

Should I route on the projected graph or the lat-lon one? Distance routing gives the identical answer either way, because length is a stored geodesic attribute rather than something measured off the coordinates. Everything else needs metres: nearest_nodes returns degrees as distances on an unprojected graph (and needs scikit-learn rather than SciPy to do it), the A* heuristic is nonsense, and any buffer or area you derive from the result is wrong. Project once, immediately after truncating.

Why does my route ignore turn restrictions and traffic lights? Because the graph has no concept of them. Costs live on edges, so a solver cannot charge for the manoeuvre between two edges — no left-turn penalty, no banned turn, no signal delay. Approximating them means expanding to a line graph where each turn becomes its own edge, which multiplies the graph size; past that point a purpose-built engine is the honest answer, and the trade-off is laid out in OSMnx vs pgRouting for network analysis.

How do I get a second-best or alternative route? ox.routing.k_shortest_paths(streets, orig, dest, k=3, weight="travel_time") yields routes in ascending cost order using Yen's algorithm. Expect near-duplicates: the second-best path often differs from the best by a single block, so filter the results on geometric dissimilarity rather than presenting them as genuine alternatives.

Can I route to a point that is not an intersection? Not directly — every node in the graph is a junction. ox.distance.nearest_edges gives you the closest street segment instead, and you can interpolate the position along it and add the partial edge cost to the route total by hand. The pragmatic alternative is to accept the snap error and report it: on a dense urban network it is typically 10–50 m, which is smaller than the uncertainty in the imputed speeds anyway.