Routing and Shortest Paths

5.4. Routing and Shortest Paths#

5.4.1. Learning objectives#

By the end of this section, you should be able to:

  • explain a shortest path as a minimum-cost path and choose a cost that matches the question;

  • snap coordinates to graph nodes and check the quality of the snap;

  • compare a distance route with a time route and report both correctly on a multigraph;

  • handle origins and destinations that cannot be connected;

  • interpret a detour factor and say why a nearest-neighbour stop order is not optimal; and

  • validate and document a route.

Section goal

Turn a pair of coordinates into a route that can be defended. A small synthetic city with known lengths and speeds shows how cost, snapping, parallel edges, and no-path cases behave, and the same code then runs on the Columbia graph.

5.4.2. Define the route cost#

A shortest route is only meaningful after choosing a cost. length answers a distance question; travel_time answers a time question. Other costs can represent tolls, slope, risk, or a policy penalty, provided the edge values and units are documented. Table 5.14 compares them.

Table 5.14 Route cost fields and what each one answers.#

Cost field

Question

Unit

Caution

length

Which route is shortest in distance?

Metres

Ignores speed and delay

travel_time

Which route is fastest?

Seconds

Estimated from speed limits and imputed speeds, not from observed traffic

A custom cost

Which route avoids a penalty such as slope, risk, or tolls?

Documented by the analyst

Values must be recorded and must not be negative

A shortest path is the path whose edge costs add up to the smallest total. OSMnx and NetworkX find it with Dijkstra’s algorithm [Dijkstra, 1959], which grows outward from the origin and always extends the cheapest unexplored path first. That strategy is only correct when no cost is negative, so a custom penalty must never subtract.

Travel time is derived. OSMnx [Boeing, 2025] reads each edge’s maxspeed tag where one exists and fills the others with an estimate based on the road type, stores the result as speed_kph, and computes travel_time in seconds as the edge length divided by the speed.

The toy city below has a short, slow route and a long, fast one between an origin 1 and a destination 3. A 1,000 m residential street runs straight between them at 30 km/h. A 1,600 m detour of primary roads runs at 80 km/h, and its long section has a parallel, shorter but slower tertiary road. Nodes 6 and 7 form a separate pair of streets that are not connected to the rest.

Hide code cell source

 
import math

import networkx as nx
import osmnx as ox

ORIGIN, MID, DESTINATION, NORTH_WEST, NORTH_EAST, ISLAND_A, ISLAND_B = 1, 2, 3, 4, 5, 6, 7
positions = {
    ORIGIN: (0, 0), MID: (500, 0), DESTINATION: (1000, 0),
    NORTH_WEST: (0, 300), NORTH_EAST: (1000, 300), ISLAND_A: (5000, 0), ISLAND_B: (5100, 0),
}

city = nx.MultiDiGraph(crs="EPSG:32615")
for node, (x, y) in positions.items():
    city.add_node(node, x=500_000 + x, y=4_300_000 + y)


def add_street(a, b, highway, maxspeed, length=None):
    if length is None:
        length = math.dist((city.nodes[a]["x"], city.nodes[a]["y"]), (city.nodes[b]["x"], city.nodes[b]["y"]))
    for u, v in ((a, b), (b, a)):
        city.add_edge(u, v, length=length, highway=highway, oneway=False, maxspeed=maxspeed)


add_street(ORIGIN, MID, "residential", "30")
add_street(MID, DESTINATION, "residential", "30")
add_street(ORIGIN, NORTH_WEST, "primary", "80")
add_street(NORTH_EAST, DESTINATION, "primary", "80")
add_street(NORTH_WEST, NORTH_EAST, "primary", "80")                 # the fast long section
add_street(NORTH_WEST, NORTH_EAST, "tertiary", "30", length=950)    # a shorter, slower parallel road
add_street(ISLAND_A, ISLAND_B, "residential", "30")

city = ox.add_edge_speeds(city)
city = ox.add_edge_travel_times(city)

print(f"{'edge':>8} {'type':>12} {'length m':>9} {'km/h':>5} {'seconds':>8}")
for u, v, key, data in city.edges(keys=True, data=True):
    if u < v:
        print(f"{u}->{v} #{key:<1} {data['highway']:>12} {data['length']:>9.0f} {data['speed_kph']:>5.0f} {data['travel_time']:>8.1f}")

    edge         type  length m  km/h  seconds
1->2 #0  residential       500    30     60.0
1->4 #0      primary       300    80     13.5
2->3 #0  residential       500    30     60.0
3->5 #0      primary       300    80     13.5
4->5 #0      primary      1000    80     45.0
4->5 #1     tertiary       950    30    114.0
6->7 #0  residential       100    30     12.0

Each travel_time is the length divided by the speed: the 500 m residential edge at 30 km/h takes 60 s, and the 1,000 m primary edge at 80 km/h takes 45 s. The two parallel edges between nodes 4 and 5 differ in length and in time, which matters below. On the Columbia graph the same two calls add the speeds and travel times. OSMnx estimates missing speeds, so treat the resulting travel times as estimates, not observations, and inspect the imputed values before making a decision.

The cell below loads the drivable network of Columbia from S01 so that this section runs on its own. The first download is cached.

Hide code cell source

 
import osmnx as ox

ox.settings.use_cache = True

G = ox.graph_from_place("Columbia, Missouri, USA", network_type="drive", simplify=True, retain_all=False)
print(G.number_of_nodes(), G.number_of_edges(), G.graph.get("crs"))

5766 13888 epsg:4326

Hide code cell source

 
import osmnx as ox

G_routing = ox.add_edge_speeds(G)
G_routing = ox.add_edge_travel_times(G_routing)

5.4.3. Snap locations to the graph#

Coordinates describe arbitrary locations, while routing requires graph nodes. Table 5.15 lists the two ways to snap.

Table 5.15 Ways to snap a coordinate to a graph.#

Method

Returns

Strength

Caution

ox.distance.nearest_nodes

The closest node

Fast, and routing needs a node

Error up to about half a block, and the node can lie across a barrier

ox.distance.nearest_edges

The closest edge and the point on it

Locates the point on the street itself

A new node must be inserted into the graph before routing

Whichever method is used, ask for the snap distance with return_dist=True. A snap of a few metres is routine, whereas a snap of hundreds of metres means the point is outside the network or the CRS is wrong. The graph and the coordinates must use compatible CRS conventions: a projected graph needs projected coordinates, and an unprojected graph needs longitude and latitude. The cell below supplies the same location both ways to a projected graph.

Hide code cell source

 
import geopandas as gpd
from shapely.geometry import Point

# A location a user might supply: a point near node 4, given in longitude and latitude.
x, y = 500_010, 4_300_290
as_lonlat = gpd.GeoSeries([Point(x, y)], crs=city.graph["crs"]).to_crs("EPSG:4326").iloc[0]

wrong_node, wrong_dist = ox.distance.nearest_nodes(city, X=as_lonlat.x, Y=as_lonlat.y, return_dist=True)
print(f"lon/lat passed to a projected graph: node {wrong_node}, snap distance {wrong_dist / 1000:,.0f} km")

projected = gpd.GeoSeries([as_lonlat], crs="EPSG:4326").to_crs(city.graph["crs"]).iloc[0]
node, dist = ox.distance.nearest_nodes(city, X=projected.x, Y=projected.y, return_dist=True)
print(f"projected first:                    node {node}, snap distance {dist:,.1f} m")

lon/lat passed to a projected graph: node 1, snap distance 4,329 km
projected first:                    node 4, snap distance 14.1 m

Passing longitude and latitude to a projected graph still returns a node, but one thousands of kilometres away, and only the snap distance reveals it. Projecting the point first gives a 14 m snap to node 4. Check the same things on the Columbia graph, which is unprojected and therefore takes longitude and latitude directly (OSMnx needs the scikit-learn package to snap on an unprojected graph). Check that the nearest nodes fall inside the intended study area and that the selected network type matches the traveler:

Hide code cell source

 
origin = (38.9517, -92.3341)
destination = (38.9466, -92.3277)

origin_node = ox.distance.nearest_nodes(
    G_routing, X=origin[1], Y=origin[0]
)
destination_node = ox.distance.nearest_nodes(
    G_routing, X=destination[1], Y=destination[0]
)

5.4.4. Compare distance and time routes#

ox.shortest_path takes the graph, the two nodes, and the weight to minimize. Different weights can choose different routes between the same two nodes, which is the point of comparing them. The toy city has such a pair:

Hide code cell source

 
distance_route = ox.shortest_path(city, ORIGIN, DESTINATION, weight="length")
time_route = ox.shortest_path(city, ORIGIN, DESTINATION, weight="travel_time")

print("distance route nodes:", distance_route)
print("time route nodes:    ", time_route)


def route_cost(graph, route, attribute):
    total = 0.0
    for u, v in zip(route[:-1], route[1:]):
        choices = graph[u][v].values()
        total += min(data[attribute] for data in choices)
    return total


print()
print("distance route: length", route_cost(city, distance_route, "length"), "m, time",
      route_cost(city, distance_route, "travel_time"), "s")
print("time route:     length", route_cost(city, time_route, "length"), "m (wrong), time",
      route_cost(city, time_route, "travel_time"), "s")

distance route nodes: [1, 2, 3]
time route nodes:     [1, 4, 5, 3]

distance route: length 1000.0 m, time 120.0 s
time route:     length 1550.0 m (wrong), time 72.0 s

The distance route takes the 1,000 m street in 120 s, and the time route takes the 1,600 m primary detour in 72 s, 60% longer in distance and 40% shorter in time. Neither is wrong; they answer different questions.

The last line exposes a trap on multigraphs. For a multigraph, calculate a route’s edge attributes using the edge selected by the algorithm rather than assuming that key 0 is always the correct parallel edge. route_cost selects the least-cost parallel edge, which is correct when the attribute is the weight the route was found with. For another attribute it can pick a different edge: the time route used the 1,000 m primary road between nodes 4 and 5, but the least length is the slower 950 m tertiary road, so route_cost reports 1,550 m instead of 1,600 m. The function below selects the edge by the routing weight and reports both attributes of that same edge.

Hide code cell source

 
def route_summary(graph, route, weight):
    length = seconds = 0.0
    for u, v in zip(route[:-1], route[1:]):
        edge = min(graph[u][v].values(), key=lambda data: data[weight])
        length += edge["length"]
        seconds += edge["travel_time"]
    return length, seconds


print("distance route (length m, time s):", route_summary(city, distance_route, "length"))
print("time route     (length m, time s):", route_summary(city, time_route, "travel_time"))

distance route (length m, time s): (1000.0, 120.0)
time route     (length m, time s): (1600.0, 72.0)

The time route now reports 1,600 m and 72 s, using one edge for both numbers. OSMnx also offers ox.routing.route_to_gdf(G, route, weight=...), which returns the chosen edges as a table. On the Columbia graph the comparison and the cost function are:

Hide code cell source

 
import networkx as nx

distance_route = ox.shortest_path(
    G_routing,
    origin_node,
    destination_node,
    weight="length",
)
time_route = ox.shortest_path(
    G_routing,
    origin_node,
    destination_node,
    weight="travel_time",
)

print("distance route nodes:", len(distance_route))
print("time route nodes:", len(time_route))

distance route nodes: 13
time route nodes: 13

Hide code cell source

 
def route_cost(graph, route, attribute):
    total = 0.0
    for u, v in zip(route[:-1], route[1:]):
        choices = graph[u][v].values()
        total += min(data[attribute] for data in choices)
    return total

print("distance (m):", route_cost(G_routing, distance_route, "length"))
print("time (s):", route_cost(G_routing, time_route, "travel_time"))

distance (m): 1134.1625754239046
time (s): 98.9447276988518

Hide code cell source

 
print("distance route (length m, time s):", route_summary(G_routing, distance_route, "length"))
print("time route     (length m, time s):", route_summary(G_routing, time_route, "travel_time"))

distance route (length m, time s): (np.float64(1134.1625754239046), 107.34951487773185)
time route     (length m, time s): (np.float64(1135.4146831285439), 98.9447276988518)

On Columbia the two routes are nearly the same: the time route is 1.2 m longer (0.1%) and 8.4 s faster (8%), and both use 13 nodes, so the choice of cost matters little for this short trip. In the toy city the same choice changed the route completely.

The shortest-time route may be longer because a faster road can reduce total travel time. Not every pair of nodes can be connected. nx.shortest_path raises networkx.NetworkXNoPath, whereas ox.shortest_path returns None, and either way the origin should be reported as disconnected instead of silently substituting a straight-line distance. The toy city has an island of two streets that no route can reach:

Hide code cell source

 
route = ox.shortest_path(city, ORIGIN, ISLAND_A, weight="length")
print("ox.shortest_path returns:", route)

try:
    nx.shortest_path(city, ORIGIN, ISLAND_A, weight="length")
except nx.NetworkXNoPath as error:
    print("nx.shortest_path raises: NetworkXNoPath -", error)

reachable = nx.single_source_dijkstra_path_length(city, ORIGIN, weight="length")
unreachable = sorted(set(city.nodes) - set(reachable))
print("nodes unreachable from the origin:", unreachable)

ox.shortest_path returns: None
nx.shortest_path raises: NetworkXNoPath - No path between 1 and 6.
nodes unreachable from the origin: [6, 7]

A None that is passed on to later code, for example as len(route), fails far from its cause, so test for it where the route is computed. Listing the unreachable nodes, as above, is also the cheapest way to report disconnected origins in bulk.

5.4.5. Detour and multi-stop routing#

The detour factor compares a network route to a straight-line reference: the route length divided by the straight-line distance. A factor of 1 means the network offers a direct route, and larger values mean barriers or sparse links. The reference must join the same two points as the route, which are the snapped nodes, not the original coordinates, or the snap error leaks into the factor.

Hide code cell source

 
toy_distance_route = ox.shortest_path(city, ORIGIN, DESTINATION, weight="length")
toy_time_route = ox.shortest_path(city, ORIGIN, DESTINATION, weight="travel_time")

straight_line = math.dist(
    (city.nodes[ORIGIN]["x"], city.nodes[ORIGIN]["y"]),
    (city.nodes[DESTINATION]["x"], city.nodes[DESTINATION]["y"]),
)

for name, route_nodes, weight in [("distance route", toy_distance_route, "length"), ("time route", toy_time_route, "travel_time")]:
    length, _ = route_summary(city, route_nodes, weight)
    print(f"{name}: {length:,.0f} m / {straight_line:,.0f} m straight line = detour factor {length / straight_line:.2f}")

distance route: 1,000 m / 1,000 m straight line = detour factor 1.00
time route: 1,600 m / 1,000 m straight line = detour factor 1.60

The distance route has a factor of 1.00 because the street runs straight, and the time route has 1.60. A high factor is not an error by itself, but it is a prompt to ask whether the network, the cost, or the snapped points explain it. On Columbia the factor is computed from the original coordinates:

Hide code cell source

 
straight_line_m = ox.distance.great_circle(
    origin[0], origin[1], destination[0], destination[1]
)
network_m = route_cost(G_routing, distance_route, "length")
detour_factor = network_m / straight_line_m
print("detour factor:", detour_factor)

detour factor: 1.4312913941734

A factor of 1.43 means about 1,134 m of street for a straight line of 792 m. As the text above advises, the reference should join the snapped nodes, whereas this one joins the original coordinates, so any snap error is included in the factor.

For multiple destinations, create a distance matrix and then apply a route-ordering method. A nearest-neighbor order is fast but is not guaranteed to be optimal [Rosenkrantz et al., 1977]; larger problems need a documented optimization method. The toy below places a depot and four stops on one long street, so the exact answer can be found by trying every order.

Hide code cell source

 
import itertools

stop_positions = [0, 500, -600, 2000, -2400]          # metres along one street; 0 is the depot
sorted_positions = sorted(stop_positions)
street = nx.Graph()
for a, b in zip(sorted_positions[:-1], sorted_positions[1:]):
    street.add_edge(a, b, weight=b - a)
matrix = dict(nx.all_pairs_dijkstra_path_length(street))      # the distance matrix

# Nearest-neighbour order: always visit the closest unvisited stop next.
current, remaining, nn_order, nn_total = 0, set(stop_positions[1:]), [], 0
while remaining:
    nearest = min(remaining, key=lambda stop: matrix[current][stop])
    nn_total += matrix[current][nearest]
    nn_order.append(nearest)
    remaining.remove(nearest)
    current = nearest

# Exact order: try every permutation, which is only feasible for a handful of stops.
def total(order):
    return sum(matrix[a][b] for a, b in zip((0,) + order, order))

best_order = min(itertools.permutations(stop_positions[1:]), key=total)
print("nearest-neighbour:", nn_order, f"{nn_total:,} m")
print("optimal:          ", list(best_order), f"{total(best_order):,} m")

nearest-neighbour: [500, -600, -2400, 2000] 7,800 m
optimal:           [500, 2000, -600, -2400] 6,400 m

The nearest-neighbour rule is lured by the closest stop each time. It visits the stop 500 m east, then turns west for two stops, and finally has to cross the whole street to reach the far east stop, for 7,800 m in total. The optimal order finishes the east side first and then the west, for 6,400 m, so the nearest-neighbour route is about 22% longer. Trying every order grows factorially, so beyond roughly nine stops a heuristic or a solver is needed, and its method and settings belong in the record. On the Columbia graph, a distance matrix is built with one Dijkstra search from each stop:

Hide code cell source

 
stops = [origin_node, destination_node]  # add more snapped stop nodes as needed
distance_matrix = {
    stop: nx.single_source_dijkstra_path_length(
        G_routing,
        stop,
        weight="travel_time",
    )
    for stop in stops
}

5.4.6. Visualize and validate#

A route that looks plausible on a map is not yet validated. Check that it starts and ends at the snapped nodes, that every consecutive pair of nodes is joined by an edge, that it does not revisit a node, that it is not None, and that the cost recomputed from the edges matches the cost the algorithm reports. The toy city shows both routes and runs those checks.

../_images/routing_toy_routes.png

Fig. 5.3 The distance route (blue) and the travel-time route (red) between the origin and destination of the toy city; the unreachable island is left out.#

In Fig. 5.3, the blue distance route follows the straight 1,000 m residential street along the bottom. The red time route climbs the 300 m primary road, crosses the 1,000 m primary road at the top, and descends the second 300 m primary road, 1,600 m in all. The two parallel roads between the top corners overlap on the map, so the figure cannot show which one was chosen. The validation below reads the chosen edges directly, which is the multigraph trap that the cost functions above are written to handle.

Hide code cell source

 
def validate_route(graph, route, origin_node, destination_node, weight):
    assert route is not None, "no route was found"
    assert route[0] == origin_node and route[-1] == destination_node, "route does not join the snapped nodes"
    assert all(graph.has_edge(u, v) for u, v in zip(route[:-1], route[1:])), "route uses a missing edge"
    assert len(set(route)) == len(route), "route revisits a node"
    recomputed = route_cost(graph, route, weight)
    reported = nx.shortest_path_length(graph, origin_node, destination_node, weight=weight)
    assert math.isclose(recomputed, reported), f"recomputed {recomputed} != reported {reported}"
    return recomputed


for name, route_nodes, weight in [("distance route", toy_distance_route, "length"),
                                  ("time route", toy_time_route, "travel_time")]:
    cost = validate_route(city, route_nodes, ORIGIN, DESTINATION, weight)
    print(f"{name}: passed all checks, {weight} = {cost:,.1f}")

distance route: passed all checks, length = 1,000.0
time route: passed all checks, travel_time = 72.0

route_cost is valid here because each route is checked against the weight it was found with. The function raises an AssertionError on the first failed check, so a route that prints has passed all of them. To see which parallel edge the time route used, list the chosen edges:

Hide code cell source

 
edges = ox.routing.route_to_gdf(city, toy_time_route, weight="travel_time")
print(edges[["highway", "length", "travel_time"]])

         highway  length  travel_time
u v key                              
1 4 0    primary   300.0         13.5
4 5 0    primary  1000.0         45.0
5 3 0    primary   300.0         13.5

The table shows that the segment between nodes 4 and 5 is the 1,000 m primary road at 45 s, not the 950 m tertiary road at about 114 s. The figure could not show this.

On the Columbia graph the routes are drawn and checked the same way. The two routes differ by about a metre, so the distance route is drawn wider and underneath to keep both visible:

Hide code cell source

 
fig, ax = ox.plot_graph_routes(
    G_routing, [distance_route, time_route],
    route_colors=["tab:blue", "tab:red"], route_linewidths=[6, 3],
    node_size=0, show=False, close=False,
)
ax.set_title("Columbia: distance route (blue) and travel-time route (red)")

Text(0.5, 1.0, 'Columbia: distance route (blue) and travel-time route (red)')
../_images/2ecab514cf9f1cee26846c3b5048e4da205b2de88ae75a7b40abd89c82429c42.png

Hide code cell source

 
for name, route_nodes, weight in [("distance route", distance_route, "length"),
                                  ("time route", time_route, "travel_time")]:
    cost = validate_route(G_routing, route_nodes, origin_node, destination_node, weight)
    print(f"{name}: passed all checks, {weight} = {cost:,.1f}")

distance route: passed all checks, length = 1,134.2
time route: passed all checks, travel_time = 98.9

Record the origin and destination, graph mode, cost field, route length, travel time, and any missing or imputed edge values. A route that passes these checks is internally consistent. It is only as accurate as the speeds behind it, which here are speed-limit tags and road-type estimates rather than observed traffic.

5.4.7. Summary#

Routing is a weighted graph problem, not a nearest-point calculation. Snap locations carefully and check the snap distance, choose a cost that matches the question, preserve directionality, read parallel edges by the routing weight, handle no-path cases, and report the assumptions behind every route, including the stop-ordering method when there is more than one destination.