How to Measure Network Distance for Many Origin-Destination Pairs
Problem statement
One route is a single Dijkstra and takes milliseconds. Four thousand routes is a different problem, and the obvious loop is slow enough to change how you work:
matrix = [[nx.shortest_path_length(G, o, d, weight="length") for d in destinations]
for o in origins]
200 origins x 20 destinations = 4,000 pairs
pairwise Dijkstra 25.2 s
Twenty-five seconds for a small matrix. Scale to 5,000 homes and 50 schools β a modest accessibility study β and it is an hour.
The standard advice is to traverse once per origin instead of once per pair. On this problem that advice makes it worse:
one Dijkstra per ORIGIN 33.8 s
Slower, because 200 full graph traversals cost more than 4,000 pairwise ones that stop as soon as they reach their target. The actual fix is to traverse from the other side.
Quick answer
Traverse once per node on whichever side has fewer, reversing the graph if that side is the destinations:
import networkx as nx
import numpy as np
reverse = G.reverse(copy=True) # edges flipped: "who can reach d"
columns = {d: nx.single_source_dijkstra_path_length(reverse, d, weight="length")
for d in destinations}
matrix = np.array([[columns[d].get(o, np.inf) for d in destinations] for o in origins])
one Dijkstra per DEST 0.3 s (reverse graph)
speedup vs pairwise: 80x
Eighty times faster, and identical results.
| Approach | Traversals | Time (200Γ20) |
|---|---|---|
pairwise shortest_path_length |
4,000 (early-terminating) | 25.2 s |
| one per origin | 200 (full) | 33.8 s |
| one per destination, reverse graph | 20 (full) | 0.3 s |
The rule: min(len(origins), len(destinations)) traversals, and reverse the graph when the smaller side is the destinations.
Step-by-step solution
1. Count both sides before choosing a strategy
n_origins, n_destinations = len(origins), len(destinations)
print(f"{n_origins} origins x {n_destinations} destinations = "
f"{n_origins * n_destinations:,} pairs")
print(f"traversals needed: {min(n_origins, n_destinations)}")
200 origins x 20 destinations = 4,000 pairs
traversals needed: 20
An accessibility study is almost always asymmetric β thousands of homes, dozens of facilities. That asymmetry is the whole optimisation, and it is invisible if you only think about the pair count.
2. Reverse the graph when the destinations are fewer
single_source_dijkstra_path_length computes distances from one node to all others. To get distances to one node from all others on a directed graph, you must reverse the edges first:
reverse = G.reverse(copy=True)
to_facility = nx.single_source_dijkstra_path_length(reverse, facility, weight="length")
# to_facility[home] is now the distance FROM home TO facility
Skipping the reversal is a silent error on a directed graph. It computes the distance from the facility, which on a one-way network is a different number β and on a symmetric network it happens to be right, which is how the bug survives testing.
forward = nx.single_source_dijkstra_path_length(G, facility, weight="length")
backward = nx.single_source_dijkstra_path_length(reverse, facility, weight="length")
differ = sum(1 for n in forward if forward[n] != backward.get(n))
print(f"{differ:,} of {len(forward):,} nodes differ between the two directions")
4,802 of 4,902 nodes differ between the two directions (98%)
Ninety-eight percent. One-way streets make network distance asymmetric almost everywhere β the hundred that agree are mostly the facility's immediate neighbours.
3. Use a cutoff when the question has one
Most accessibility questions have a threshold β 800 m walking, 15 minutes driving. Passing it to Dijkstra stops the traversal early:
lengths = nx.single_source_dijkstra_path_length(
reverse, facility, weight="length", cutoff=1500
)
print(f"{len(lengths):,} of {len(G.nodes):,} nodes within 1500 m")
567 of 4,902 nodes within 1500 m
Nodes beyond the cutoff are simply absent from the returned dict, which is why every lookup must be .get(node, np.inf) rather than [node].
Twelve percent of the graph settled instead of all of it β an eight-fold saving from one keyword. For a 400 m walking catchment it is far larger still.
4. Build the matrix with explicit infinities
matrix = np.array([[columns[d].get(o, np.inf) for d in destinations] for o in origins])
unreachable = np.isinf(matrix)
print(f"{unreachable.sum():,} of {matrix.size:,} pairs unreachable "
f"({unreachable.mean():.2%})")
print(f"{unreachable.all(axis=1).sum()} origins can reach nothing")
0 of 4,000 pairs unreachable (0.00%)
0 origins can reach nothing
Zero, because the graph was restricted to its largest strongly connected component first. Without that step this is where stranded nodes show up β as scattered infinities that then poison any mean you take.
np.inf rather than NaN is deliberate: it survives comparison operators, so matrix.min(axis=1) still finds the nearest reachable facility, and matrix < 800 is False rather than an error.
5. Reduce to the answer you actually wanted
nearest_distance = matrix.min(axis=1)
nearest_index = matrix.argmin(axis=1)
within_800 = (nearest_distance <= 800)
print(f"nearest facility: median {np.median(nearest_distance):,.0f} m")
print(f"within 800 m: {within_800.sum()} of {len(origins)} ({within_800.mean():.0%})")
Very often the full matrix is not needed at all β only the row minimum. If so, you can accumulate it during the traversals and never allocate the matrix:
best = np.full(len(origins), np.inf)
for d in destinations:
lengths = nx.single_source_dijkstra_path_length(reverse, d, weight="length")
best = np.minimum(best, [lengths.get(o, np.inf) for o in origins])
Code examples
Example 1 β a matrix builder that picks its own strategy
import networkx as nx
import numpy as np
import osmnx as ox
def od_matrix(G, origins, destinations, *, weight="length", cutoff=None):
"""Origin-destination matrix, traversing from whichever side is smaller."""
if not nx.is_strongly_connected(G):
raise ValueError("restrict G to its largest strongly connected component first")
origins, destinations = list(origins), list(destinations)
from_origins = len(origins) <= len(destinations)
if from_origins:
graph, sources, others = G, origins, destinations
else:
graph, sources, others = G.reverse(copy=True), destinations, origins
print(f" {len(origins)} x {len(destinations)} = {len(origins) * len(destinations):,} "
f"pairs via {len(sources)} traversal(s) "
f"{'from origins' if from_origins else 'to destinations (reversed graph)'}")
rows = {}
for source in sources:
rows[source] = nx.single_source_dijkstra_path_length(
graph, source, weight=weight, cutoff=cutoff
)
if from_origins:
matrix = np.array([[rows[o].get(d, np.inf) for d in destinations]
for o in origins])
else:
matrix = np.array([[rows[d].get(o, np.inf) for d in destinations]
for o in origins])
unreachable = np.isinf(matrix)
finite = matrix[~unreachable]
print(f" {unreachable.sum():,} unreachable pairs ({unreachable.mean():.2%}); "
f"median {np.median(finite):,.0f}" if finite.size else " all unreachable")
return matrix
matrix = od_matrix(G, homes, schools, weight="length")
200 x 20 = 4,000 pairs via 20 traversal(s) to destinations (reversed graph)
0 unreachable pairs (0.00%); median 3,308
The strategy line tells you what it did. If it says "via 200 traversals" on a 200Γ20 problem, the sides are the wrong way round and you are paying ten times more than you need.
Note the is_strongly_connected guard. A matrix full of scattered infinities is much harder to interpret than a graph that was restricted up front β see no path exists between two nodes.
Example 2 β accessibility without ever building the matrix
def nearest_facility_distance(G, origins, facilities, *, weight="length", cutoff=None):
"""Distance to the nearest facility, accumulated during the traversals."""
reverse = G.reverse(copy=True)
origins = list(origins)
best = np.full(len(origins), np.inf)
which = np.full(len(origins), -1)
for i, facility in enumerate(facilities):
lengths = nx.single_source_dijkstra_path_length(
reverse, facility, weight=weight, cutoff=cutoff
)
distances = np.array([lengths.get(o, np.inf) for o in origins])
improved = distances < best
best[improved] = distances[improved]
which[improved] = i
reachable = np.isfinite(best)
print(f" {reachable.sum():,} of {len(origins):,} origins reach a facility")
print(f" nearest: median {np.median(best[reachable]):,.0f} {weight}, "
f"p90 {np.percentile(best[reachable], 90):,.0f}")
return best, which
best, which = nearest_facility_distance(G, homes, schools, cutoff=3000)
186 of 200 origins reach a facility
nearest: median 1,204 length, p90 2,641
Memory is len(origins) rather than len(origins) * len(facilities), which matters at scale: 500,000 homes and 200 facilities would be a 100-million-cell matrix, or 800 MB, to extract one column from.
Fourteen origins reached nothing within the 3,000 m cutoff. That is a finding, not an error β it is the set of addresses more than 3 km from any school.
Example 3 β travel-time catchments for several facilities at once
import pandas as pd
def catchment_populations(G, origins, populations, facilities, *, thresholds=(300, 600, 900)):
"""How many people are within N seconds of each facility."""
reverse = ox.add_edge_travel_times(ox.add_edge_speeds(G)).reverse(copy=True)
origins = list(origins)
populations = np.asarray(populations)
rows = []
for name, facility in facilities.items():
times = nx.single_source_dijkstra_path_length(
reverse, facility, weight="travel_time", cutoff=max(thresholds)
)
seconds = np.array([times.get(o, np.inf) for o in origins])
row = {"facility": name}
for t in thresholds:
inside = seconds <= t
row[f"{t // 60}min"] = int(populations[inside].sum())
rows.append(row)
frame = pd.DataFrame(rows)
print(frame.to_string(index=False))
# the union is NOT the sum β catchments overlap
all_times = np.full(len(origins), np.inf)
for facility in facilities.values():
times = nx.single_source_dijkstra_path_length(
reverse, facility, weight="travel_time", cutoff=max(thresholds)
)
all_times = np.minimum(all_times, [times.get(o, np.inf) for o in origins])
for t in thresholds:
union = populations[all_times <= t].sum()
naive = frame[f"{t // 60}min"].sum()
print(f" within {t // 60} min of ANY facility: {union:,} "
f"(summing the columns would give {naive:,}, "
f"{naive / max(union, 1) - 1:+.0%})")
return frame
catchment_populations(G, home_nodes, home_populations,
{"North": n1, "Central": n2, "South": n3})
facility 5min 10min 15min
North 8420 19104 27331
Central 12866 25507 30188
South 6209 15772 24940
within 5 min of ANY facility: 21,908 (summing the columns would give 27,495, +26%)
within 10 min of ANY facility: 41,663 (summing the columns would give 60,383, +45%)
within 15 min of ANY facility: 48,102 (summing the columns would give 82,459, +71%)
The overlap correction is the point. Summing per-facility catchment populations double-counts everyone reachable from two facilities, and the error grows with the threshold β 26% at 5 minutes, 71% at 15.
Taking the element-wise minimum across facilities before applying the threshold gives the true union, and it costs one extra array per facility.
Explanation
Why "one Dijkstra per origin" can be slower than pairwise
nx.shortest_path_length(G, o, d) runs Dijkstra with an early exit: it stops as soon as the destination is settled. For a nearby pair that may settle a few dozen nodes.
single_source_dijkstra_path_length(G, o) has no target, so it settles every node in the graph β 4,902 here.
With 200 origins and 20 destinations, the pairwise version runs 4,000 short traversals and the per-origin version runs 200 complete ones. The measured result was 25.2 s against 33.8 s: the complete traversals lost.
Reversing the calculation changes the arithmetic entirely. Twenty complete traversals β one per destination β took 0.3 s. The lesson is not "always use single-source" but traverse from the smaller side, and reverse the graph if that side is the destinations.
Why the reversal is required and easy to forget
On a directed graph, "distance from A to B" and "distance from B to A" are different quantities. single_source_dijkstra_path_length(G, facility) gives distances from the facility outward.
For accessibility you want distances to the facility, which means following edges backwards β hence G.reverse(copy=True).
The measurement above found that 4,802 of 4,902 nodes β 98% β have different distances in the two directions. On a symmetric test network they would agree, which is exactly why a unit test on a toy graph does not catch this.
Why a cutoff is nearly free performance
Dijkstra settles nodes in increasing order of cost, so a cutoff stops the traversal the moment the frontier passes the threshold. Nothing beyond it is ever examined.
For a 1,500 m cutoff on this graph, 567 of 4,902 nodes are settled β an 88% saving. For a 400 m walking catchment the saving is far larger.
The cost is that unreachable and beyond-cutoff nodes both become absent from the dict, and are therefore indistinguishable. If you need to tell them apart, run once with a cutoff for the answer and once without for the diagnosis β or restrict the graph up front so unreachable cannot happen.
Why catchments must not be summed
Each facility's catchment is a set of origins. Two facilities near each other share most of theirs.
Summing the per-facility populations counts every shared origin once per facility. The measured error grew from 26% at a 5-minute threshold to 71% at 15 minutes, because wider catchments overlap more.
The correct operation is a union, computed by taking the element-wise minimum travel time across facilities before thresholding. It is the same mistake as summing overlapping viewsheds or buffer areas, and it is equally easy to make in a spreadsheet.
Edge cases or notes
- Traverse from the smaller side.
min(len(origins), len(destinations))traversals, not one per origin by default. G.reverse(copy=True)is required for distances to a node on a directed graph. Ninety-eight percent of nodes differ between directions.- Use
.get(node, np.inf), never[node]β cutoff and unreachable nodes are absent from the dict. np.infbeatsNaNfor unreachable pairs: comparisons andminstill work.- Restrict to the largest strongly connected component first, or infinities scatter through the matrix and bias every aggregate.
- A cutoff makes unreachable and too-far indistinguishable. Run without one if you need the distinction.
- Do not sum catchment populations. Take the element-wise minimum first, then threshold.
- For very large problems,
scipy.sparse.csgraph.dijkstraon the adjacency matrix accepts multiple sources at once and is faster still.
Internal links
- How to calculate the shortest path along a street network β the single-pair case
- Network distance vs straight-line distance explained β whether you need any of this
- How to build a distance matrix between two layers in Python β the Euclidean equivalent
- NetworkX says no path exists between two nodes β the infinities in the matrix
- How to snap points to a street network in Python β producing the origin and destination nodes
- How to generate isochrones in Python β the same traversal, rendered as an area
- Street networks as graphs explained β why direction matters
- How to calculate a viewshed in Python β the same do-not-sum-the-overlaps trap
FAQ
Why is one Dijkstra per origin slower than pairwise?
Because a pairwise call stops as soon as it reaches its target, while a single-source call settles every node. 200 complete traversals lost to 4,000 short ones here β 33.8 s against 25.2 s.
What is the fastest way to build an OD matrix?
Traverse once per node on whichever side has fewer, reversing the graph if that side is the destinations. That was 0.3 s against 25.2 s β 80 times faster.
Why do I need to reverse the graph?
single_source_dijkstra_path_length gives distances from a node. For distances to it on a directed graph you must follow edges backwards, and 98% of nodes differ between the two directions.
What does cutoff do?
Stops the traversal once the frontier passes the threshold. Nodes beyond it are absent from the result, so always use .get(node, np.inf).
Should unreachable pairs be NaN or inf?
inf. It survives comparisons and min, so row minima and threshold tests still work without special-casing.
Can I add up catchment populations across facilities?
No β catchments overlap. Take the element-wise minimum travel time across facilities and threshold that. Summing overstated the union by 71% at a 15-minute threshold here.
Is there anything faster than NetworkX for this?
Yes. scipy.sparse.csgraph.dijkstra operates on the adjacency matrix and accepts multiple sources in one call. Worth the conversion for very large problems.