How to Measure Network Distance for Many Origin-Destination Pairs

On this page

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:

t = time.perf_counter()
pairwise = [[nx.shortest_path_length(G, o, d, weight="length") for d in destinations]
            for o in origins]
pairwise_s = time.perf_counter() - t
print(f"{len(origins)} origins x {len(destinations)} destinations = "
      f"{len(origins) * len(destinations):,} pairs")
print(f"  pairwise Dijkstra   {pairwise_s:.1f} s")

200 origins x 20 destinations = 4,000 pairs
  pairwise Dijkstra   11.2 s

Twelve seconds for a small matrix. Scale to 5,000 homes and 50 schools — a modest accessibility study, 62 times the pairs — and it is over twelve minutes.

The standard advice is to traverse once per origin instead of once per pair. It helps — about nine-fold:

t = time.perf_counter()
rows = [nx.single_source_dijkstra_path_length(G, o, weight="length") for o in origins]
print(f"  one Dijkstra per ORIGIN   {time.perf_counter() - t:.1f} s")

  one Dijkstra per ORIGIN   1.2 s

Two hundred full traversals beat 4,000 that stop at their target, because stopping at the target is not as early as it sounds: a pairwise call here had settled 46% of the graph, on average, by the time it got there. But the origins are still the wrong side to traverse from.

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

t = time.perf_counter()
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])
elapsed = time.perf_counter() - t
print(f"  one Dijkstra per DEST    {elapsed:.1f} s  (reverse graph)")
print(f"  speedup vs pairwise: {pairwise_s / elapsed:.0f}x; "
      f"identical: {np.allclose(matrix, pairwise)}")

  one Dijkstra per DEST    0.3 s  (reverse graph)
  speedup vs pairwise: 36x; identical: True

Forty times faster than pairwise, five times faster than one per origin, and identical results.

Approach	Traversals	Time (200×20)
pairwise `shortest_path_length`	4,000 (early-terminating)	11.2 s
one per origin	200 (full)	1.2 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.

One per origin already beats pairwise. One per destination, on the reversed graph, beats both.

Count both sides, reverse, cut off, use infinities

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,227 of 4,985 nodes differ between the two directions

Eighty-eight percent. One-way streets make network distance asymmetric almost everywhere in a city centre.

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")

385 of 4,985 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].

Seven percent of the graph settled instead of all of it — nearly a fourteen-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%})")

nearest facility: median 809 m
within 800 m: 99 of 200 (50%)

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])

Eighty-eight percent of nodes have a different distance in the two directions. Forgetting the reversal is silent.

Accessibility, catchments and a SciPy matrix

The data used on this page

The drive network within 2.5 km of Manchester Town Hall from OpenStreetMap, restricted to its largest strongly connected component. The origins are the 2021 Census output areas whose population-weighted centroids (ONS) fall within 2 km, with their usual-resident counts from Nomis table TS001; 200 of them are sampled for the matrix. The destinations are OpenStreetMap amenity=school features within 2 km, 20 sampled. Everything is snapped to its nearest graph node. Timings were taken at a load average of 0.2–0.5 on a 24-thread machine; they vary with load, the ratios much less. Overpass via the maps.mail.ru mirror while overpass-api.de was refusing connections. The OpenStreetMap answers are read as the Overpass mirror gave them on 10 October 2026.

import io
import os
import time
import zipfile

import geopandas as gpd
import networkx as nx
import numpy as np
import osmnx as ox
import pandas as pd
import requests
from shapely.geometry import Point

# OpenStreetMap as Overpass returned it on 10 October 2026, saved as osmnx's own cache so the counts
# below do not drift as OSM is edited; delete the cache folder to query live instead
if not os.path.exists("cache"):
    reply = requests.get("https://data.spatialworkflow.io/osm-cache-network-distance-matrix-python-2026-10-10.zip", timeout=300)
    reply.raise_for_status()
    with zipfile.ZipFile(io.BytesIO(reply.content)) as archive:
        archive.extractall("cache")

ox.settings.overpass_url = "https://maps.mail.ru/osm/tools/overpass/api"
ox.settings.overpass_rate_limit = False
ox.settings.http_user_agent = "spatialworkflow.io guide examples (https://www.spatialworkflow.io/)"
centre = (53.4808, -2.2426)
G = ox.graph_from_point(centre, dist=2500, network_type="drive")
G = ox.truncate.largest_component(G, strongly=True)
Gp = ox.project_graph(G, to_crs=27700)
town_hall = gpd.GeoSeries([Point(centre[1], centre[0])], crs=4326).to_crs(27700).iloc[0]

ONS = ("https://services1.arcgis.com/ESMARspQHYMw9BZ9/arcgis/rest/services/"
       "OA_December_2021_EW_PWC_V4/FeatureServer/0/query")
reply = requests.get(ONS, params={"geometry": f"{centre[1]},{centre[0]}", "geometryType": "esriGeometryPoint",
                                  "inSR": 4326, "distance": 2000, "units": "esriSRUnit_Meter",
                                  "outFields": "OA21CD", "outSR": 27700, "f": "geojson"}, timeout=120).json()
oas = gpd.GeoDataFrame.from_features(reply["features"], crs=27700).set_index("OA21CD").sort_index()
population = pd.read_csv("https://www.nomisweb.co.uk/api/v01/dataset/NM_2021_1.data.csv?date=latest"
                         f"&geography={','.join(oas.index)}&c2021_restype_3=0&measures=20100"
                         "&select=geography_code,obs_value").set_index("GEOGRAPHY_CODE")["OBS_VALUE"]
oas["population"] = population
oas["node"] = ox.nearest_nodes(Gp, oas.geometry.x, oas.geometry.y)

school_points = ox.features_from_point(centre, tags={"amenity": "school"}, dist=2000)
school_points = school_points[school_points["name"].notna()].reset_index().to_crs(27700)
school_points["geometry"] = school_points.geometry.representative_point()
school_points = school_points.sort_values(["element", "id"]).reset_index(drop=True)
school_points["node"] = ox.nearest_nodes(Gp, school_points.geometry.x, school_points.geometry.y)

homes = list(oas.sample(200, random_state=0)["node"])
schools = list(school_points.sample(min(20, len(school_points)), random_state=0)["node"])
origins, destinations = homes, schools
facility = schools[0]
home_nodes, home_populations = list(oas["node"]), oas["population"].to_numpy()
print(f"graph {len(G.nodes):,} nodes; {len(oas)} output areas, {oas['population'].sum():,} residents; "
      f"{len(school_points)} named schools")

graph 4,985 nodes; 283 output areas, 81,423 residents; 28 named schools

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 2,791

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=1500)

  191 of 200 origins reach a facility
  nearest: median 759 length, p90 1,272

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.

Nine origins reached nothing within the 1,500 m cutoff. That is a finding, not an error — it is the set of output areas more than 1.5 km by road from any of the 20 sampled schools.

Example 3 — travel-time catchments for several facilities at once

import pandas as pd


def catchment_populations(G, origins, populations, facilities, *, thresholds=(120, 240, 360)):
    """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


# three real schools: the northernmost, the one nearest the Town Hall, the southernmost
ends = school_points.iloc[[school_points.geometry.y.argmax(),
                           school_points.distance(town_hall).argmin(),
                           school_points.geometry.y.argmin()]]
catchment_populations(G, home_nodes, home_populations, dict(zip(ends["name"], ends["node"])))

                                        facility  2min  4min  6min
                     Saviour CofE Primary School  1330 20911 49736
                        Abbey College Manchester  3600 51854 81423
St Chrysostom's Church of England Primary School   773 27705 64020
  within 2 min of ANY facility: 5,703 (summing the columns would give 5,703, +0%)
  within 4 min of ANY facility: 70,114 (summing the columns would give 100,470, +43%)
  within 6 min of ANY facility: 81,423 (summing the columns would give 195,179, +140%)

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 — nothing at 2 minutes, where the three catchments do not yet touch; 43% at 4; 140% at 6, where every one of the 81,423 residents inside 2 km is within reach of at least one school and the columns add up to 195,179. The times are OSMnx's imputed free-flow speeds, so they are shorter than a real weekday drive.

Taking the element-wise minimum across facilities before applying the threshold gives the true union, and it costs one extra array per facility.

Catchments that do not touch add up; catchments that overlap do not.

Example 4 — all destinations in one call with SciPy

import scipy.sparse as sp
from scipy.sparse.csgraph import dijkstra

index = {n: i for i, n in enumerate(G.nodes)}
weights = {}
for u, v, data in G.edges(data=True):             # keep the SHORTEST of any parallel edges
    key = (index[u], index[v])
    weights[key] = min(weights.get(key, np.inf), data["length"])
rows_, cols_ = zip(*weights)
A = sp.csr_array((list(weights.values()), (rows_, cols_)), shape=(len(G), len(G)))

t = time.perf_counter()
D = dijkstra(A.T, indices=[index[d] for d in destinations])   # A.T is the reversed graph
fast = D.T[[index[o] for o in origins]]
print(f"  scipy, 20 sources in one call: {time.perf_counter() - t:.3f} s; "
      f"identical: {np.allclose(fast, matrix)}")

summed = nx.to_scipy_sparse_array(G, weight="length", format="csr")
print(f"  to_scipy_sparse_array differs in {(summed != A).nnz} cells: "
      f"it ADDS the lengths of parallel edges")

  scipy, 20 sources in one call: 0.006 s; identical: True
  to_scipy_sparse_array differs in 28 cells: it ADDS the lengths of parallel edges

About fifty times faster than the reversed NetworkX version, because the traversals run in compiled code. The trap is the conversion. nx.to_scipy_sparse_array collapses a MultiDiGraph by summing parallel edges, so a pair of nodes joined by two carriageways gets their combined length — a distance no vehicle drives. Build the matrix yourself and keep the minimum.

Early exits, reversal and cutoffs

Why an early exit saves less than it sounds

nx.shortest_path_length(G, o, d) runs Dijkstra with an early exit: it stops as soon as the destination is settled. But Dijkstra settles nodes in order of distance, so by then it has settled every node nearer than the destination — and a destination picked at random is rarely near. Across the 4,000 pairs here, a pairwise call settled 46% of the graph on average:

# rows[i] holds every distance from origins[i]; a pairwise call settles all nodes nearer than its target
settled = [sum(1 for v in rows[i].values() if v <= rows[i][d]) / len(G)
           for i in range(len(origins)) for d in destinations]
print(f"a pairwise call settles {np.mean(settled):.0%} of the graph on average")

a pairwise call settles 46% of the graph on average

single_source_dijkstra_path_length(G, o) has no target, so it settles every node — 4,985 here — and answers all 20 destinations at once.

So the pairwise version does 4,000 × 46% ≈ 1,840 graphs' worth of work, and the per-origin version 200. The measured result was 11.2 s against 1.2 s.

Reversing the calculation cuts it again. Twenty complete traversals — one per destination — took 0.3 s. The lesson is 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,227 of 4,985 nodes — 85% — 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.

Without a cutoff, one-to-many pays for the whole graph on every call — and still wins, because an early exit settles nearly half of it anyway.

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, 385 of 4,985 nodes are settled — a 92% 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 nothing at a 2-minute threshold to 43% at 4 minutes and 140% at 6, 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.

Smaller-side traversal, missing keys and summed catchments

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. Eighty-eight percent of nodes differed between directions here.
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 accepts multiple sources at once and ran about fifty times faster than NetworkX here — but build the sparse matrix yourself. nx.to_scipy_sparse_array sums parallel edges, and 28 node pairs in this graph have them.

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
How to choose facility locations with a p-median model — the matrix as an optimisation input
How to compute transit travel times from GTFS in Python — a time-dependent matrix from a timetable

FAQ

Is one Dijkstra per origin faster than pairwise?

Here, yes — 1.2 s against 11.2 s. A pairwise call stops at its target, but it has settled every nearer node by then, 46% of the graph on average. Traversing from the smaller side was faster again.

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 11.2 s — 36 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 88% of nodes differed between the two directions in central Manchester.

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 140% at a 6-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; it took 0.005 s for the 20 destinations here. Build the matrix yourself, keeping the shortest of any parallel edges — nx.to_scipy_sparse_array adds them together.