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.

Three strategies for a 200 by 20 origin-destination matrix, taking 25.2, 33.8 and 0.3 seconds.
The naive "one per origin" advice is the slowest of the three here. Traverse from the smaller side.

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])
A forward Dijkstra giving distances from a facility and a reverse-graph Dijkstra giving distances to it, differing on a one-way network.
Ninety percent of nodes have a different distance in the two directions. Forgetting the reversal is silent.

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.

A pairwise Dijkstra settling a small neighbourhood before stopping, against a single-source traversal settling every node in the graph.
The pairwise call stops early. The single-source call cannot, which is why 200 of them lose to 4,000 short ones.

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.inf beats NaN for unreachable pairs: comparisons and min still 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.dijkstra on the adjacency matrix accepts multiple sources at once and is faster still.

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.