⚙️ How-To Location & accessibility Networks & routing Joins & relationships Location & Accessibility · step 6 of 20 Updated

How to Measure Distance to the Nearest Facility for Every Home

Problem statement

The distance from each home to its nearest pharmacy, school or fire station is the input to most access studies. It is also easy to compute in a way that is quick, plausible and wrong: a straight-line nearest-neighbour join.

Measured for the 168,323 residents of Chittenden County, Vermont, in 2,241 populated census blocks, against 27 pharmacies in and within 10 km of the county:

  • Straight line: the population-weighted median distance to the nearest pharmacy was 1,381 m.
  • Road network: it was 2,342 m. For homes more than 500 m from a pharmacy, the road distance was a median 1.46 times the straight line, and up to 5.07 times.
  • 13.4% of residents had a different nearest pharmacy by road than in a straight line.

The network version is not slow. One multi-source shortest-path search from all 27 pharmacies at once took 1.1 ms with SciPy — against 54 ms for the equivalent NetworkX call and 1,253 ms for the common loop of one search per pharmacy.

Quick answer 

import numpy as np
from scipy.sparse.csgraph import dijkstra

# csr: sparse matrix of edge lengths (or seconds) between node indices
dist, _, source = dijkstra(csr.T, directed=True, indices=pharmacy_nodes,
                           min_only=True, return_predecessors=True)
metres = dist[home_nodes] + home_snap_m            # add the walk from the home to the network
nearest_node = source[home_nodes]                  # which pharmacy node won

min_only=True seeds all pharmacies at once and returns, for every node, the distance to the nearest one and which one it was. Searching on csr.T, the reversed graph, gives the distance from each home to a pharmacy, which is the direction people travel on a network with one-way streets.

Bar chart of the time to compute nearest-pharmacy network distances for 2,241 homes with SciPy multi-source Dijkstra, NetworkX multi-source Dijkstra and one NetworkX search per pharmacy.
The same distances three ways; the per-facility loop is over a thousand times slower than it needs to be.

Step-by-step solution

1. Represent homes as points inside populated blocks

Use a representative point for each populated census block, not its centroid, so every point lies inside its block. In the county that gave 2,241 points carrying 168,323 people. Weight every summary by population; a median over blocks treats a farm and a housing estate as equals.

2. Include facilities beyond the boundary

People near a county line use pharmacies across it. With the 10 km ring included, 3.4% of the county's residents turned out to have their nearest pharmacy by road outside the county. Leaving those pharmacies out overstates their distances.

3. Get the straight-line answer first, as a baseline 

import geopandas as gpd

nearest = gpd.sjoin_nearest(homes, pharmacies[["pid", "geometry"]], how="left", distance_col="d_euc")

It took 3.7 ms for 2,241 homes and found no ties. It is worth keeping, because the ratio between network and straight-line distance is the most useful single diagnostic of a network result.

4. Build a sparse graph with the minimum over parallel edges

OSMnx returns a multigraph, where two nodes can be joined by several edges — a dual carriageway, a road and a slip road. When converting to a SciPy matrix, keep the shortest of them. NetworkX's to_scipy_sparse_array sums parallel edge weights, which silently lengthens those links; building the matrix from a dictionary of minimum weights avoids it (Example 2). The county's drive graph had 9,835 nodes and 23,264 edges, and the matrix took 94 ms to build.

5. Snap homes and facilities to nodes, and keep the snap distances 

home snap median 135 m, p90 569 m, max 2,562 m; pharmacy snap median 73 m, max 540 m

The snap leg is not in the network distance, so add it. For the median resident it was 5.6% of the total distance, and for one in ten more than 15.4%.

scipy dijkstra, min_only, reversed graph (length)          1.1 ms
networkx multi_source_dijkstra_path_length (length)       54.1 ms
one networkx Dijkstra per pharmacy, then minimum        1,253.2 ms

All three gave identical distances to the metre. The loop is common in tutorials because it also yields a full origin–destination matrix; if only the nearest facility is needed, it does 27 searches to answer one question.

7. Count homes that cannot reach anything

One block, holding 229 people, returned an infinite distance: its nearest node sat on a fragment of the graph with no route to any pharmacy. Report such homes and fix the snapping or the graph, rather than letting infinity into a mean or dropping them silently.

8. Compare network and straight line 

                 straight line   road network
median distance     1,381 m        2,342 m
90th percentile     7,502 m       10,040 m
within 1 km          33.0%          14.6%
within 2 km          59.9%          44.1%
within 5 km          82.1%          69.9%

The network distance is longer everywhere, and proportionally longest at short range, where a detour around a block or a river doubles a trip. The drive time to the nearest pharmacy, from OSMnx free-flow speeds, was a median 2.8 minutes and at most 22.9.

Bar chart of the share of residents within 1, 2, 5 and 10 km of the nearest pharmacy by straight line and by road.
Straight-line distance puts more than twice as many people within a kilometre of a pharmacy as the roads do.

Code examples

Example 1 — straight-line nearest facility, with ties checked 

import geopandas as gpd


def nearest_straight_line(homes, facilities, id_col="pid"):
    """Straight-line nearest facility per home; sjoin_nearest returns one row per tie."""
    joined = gpd.sjoin_nearest(homes, facilities[[id_col, "geometry"]], how="left", distance_col="d_euc")
    ties = len(joined) - len(homes)
    if ties:
        print(f"{ties} extra rows from ties; keeping the first facility for each home")
    return joined[~joined.index.duplicated()][[id_col, "d_euc"]]

Both layers must be in the same projected CRS, in metres. Without the de-duplication, a home equidistant from two facilities appears twice and its population is counted twice in every summary.

Example 2 — network nearest facility with SciPy 

import numpy as np
import osmnx as ox
from scipy.sparse import csr_array
from scipy.sparse.csgraph import dijkstra


def graph_to_csr(G, weight="length"):
    """Sparse adjacency with the minimum weight over parallel edges."""
    nodes = list(G.nodes)
    index = {n: i for i, n in enumerate(nodes)}
    best = {}
    for u, v, data in G.edges(data=True):
        key = (index[u], index[v])
        best[key] = min(best.get(key, np.inf), data[weight])
    rows, cols = zip(*best)
    return csr_array((list(best.values()), (rows, cols)), shape=(len(nodes), len(nodes))), index


def nearest_on_network(G, homes, facilities, weight="length"):
    """Network distance from every home to its nearest facility, including both snap legs."""
    A, index = graph_to_csr(G, weight)
    h_nodes, h_snap = ox.distance.nearest_nodes(G, homes.geometry.x, homes.geometry.y, return_dist=True)
    f_nodes, f_snap = ox.distance.nearest_nodes(G, facilities.geometry.x, facilities.geometry.y, return_dist=True)
    f_idx = np.array([index[n] for n in f_nodes])
    h_idx = np.array([index[n] for n in h_nodes])
    dist, _, source = dijkstra(A.T, directed=True, indices=f_idx, min_only=True, return_predecessors=True)
    facility_snap = {}
    for i, s in zip(f_idx, f_snap):
        facility_snap[i] = min(s, facility_snap.get(i, np.inf))
    total = dist[h_idx] + h_snap + np.array([facility_snap.get(s, 0.0) for s in source[h_idx]])
    unreachable = np.isinf(total)
    print(f"{unreachable.sum()} homes cannot reach any facility ({homes.loc[unreachable, 'POP20'].sum():,} people)")
    return total, source[h_idx]
1 homes cannot reach any facility (229 people)

source holds the node index of the winning facility; map it back to facility identifiers with the same index. Two facilities snapped to one node share a source, so give them the smaller snap distance, as here.

Example 3 — compare the two, by population 

def wquantile(values, weights, q):
    order = np.argsort(values)
    cum = np.cumsum(weights[order]) / weights.sum()
    return values[order][np.searchsorted(cum, q)]


def compare_distances(straight, network, population, thresholds=(1000, 2000, 5000, 10000)):
    """Population-weighted comparison of straight-line and network nearest distances."""
    ok = np.isfinite(network)
    s, n, w = straight[ok], network[ok], population[ok]
    print(f"median: straight {wquantile(s, w, .5):,.0f} m, network {wquantile(n, w, .5):,.0f} m")
    far = s > 500
    print(f"network / straight for homes beyond 500 m: median {wquantile(n[far] / s[far], w[far], .5):.2f}")
    for t in thresholds:
        print(f"within {t / 1000:>4.0f} km: straight {w[s <= t].sum() / w.sum():.1%}, network {w[n <= t].sum() / w.sum():.1%}")
median: straight 1,381 m, network 2,342 m
network / straight for homes beyond 500 m: median 1.46
within    1 km: straight 33.0%, network 14.6%
within    2 km: straight 59.9%, network 44.1%
within    5 km: straight 82.1%, network 69.9%
within   10 km: straight 95.3%, network 90.0%

The ratio is the check: a median near 1.0 suggests the network is missing roads or the snap is too generous, and a median above about 2 suggests missing connections.

Explanation

Why one search answers the nearest-facility question

Dijkstra's algorithm grows outwards from its sources in order of distance. Starting from every facility at once, each node is first reached from whichever facility is closest, so a single run labels every node with its nearest facility and the distance to it. The cost is one search, whatever the number of facilities; the loop multiplies it by 27.

Why the reversed graph

On a directed graph, the distance from A to B is not the distance from B to A when one-way streets are involved. A search from facilities over csr gives distances from facilities to homes. Transposing the matrix reverses every edge, so the same search gives distances from homes to facilities — the trip residents make. For an emergency service driving out to homes, use csr as it is.

Why the nearest facility changes

A straight line picks the facility with the least displacement; a road picks the one with the shortest route. They differ wherever a river, rail line or disconnected street pattern makes the geometrically nearest facility a long way round. For 299 of the 2,240 reachable blocks, holding 13.4% of residents, the two answers named different pharmacies — which matters for any analysis that assigns people to facilities, not just for distances.

Why the snap leg matters

Network distance starts at a node, and the nearest node can be hundreds of metres from a rural home. In towns the leg is negligible; in the countryside it is a large share of the trip. Adding it keeps rural homes from appearing closer than they are.

Flow from homes and facilities to snapped nodes, a sparse graph, one reversed multi-source Dijkstra, and distances with snap legs added.
The snap legs go on at both ends; the search is the cheap part in the middle.

Edge cases or notes

  • Remove duplicate facilities mapped as both a point and a building before snapping, or one pharmacy counts twice.
  • Keep the largest strongly connected component if many homes return infinity, then report what was dropped.
  • Travel time needs speeds. ox.add_edge_speeds imputes missing speeds from averages by road type; state that.
  • Walking uses an undirected network. Use network_type="walk" and the reversal no longer matters.
  • Distances in degrees are meaningless. Project the graph and the points before snapping.
  • Nearest-node snapping can cross a river. Check the longest snap legs in rural areas.
  • Straight-line ties are rare with real coordinates — none here — but must be handled when facilities share an address.

FAQ

How do I find the nearest facility to each home by road in Python?

Build a sparse graph of the road network, snap homes and facilities to nodes, and run SciPy's dijkstra with all facility nodes as indices and min_only=True on the transposed matrix. For 2,241 homes and 27 pharmacies it took 1.1 ms.

Is straight-line distance close enough?

Usually not for short distances. The median resident was 1,381 m from a pharmacy in a straight line and 2,342 m by road, and straight lines put 33.0% of people within 1 km against 14.6% by road.

Why is a loop over facilities slow?

It runs one full shortest-path search per facility. The per-pharmacy NetworkX loop took 1,253 ms; one multi-source search gave identical distances in 1.1 ms.

Why do some homes get an infinite distance?

Their snapped node is on a part of the network with no route to any facility. One block with 229 residents did here; fix the snap or the graph and report it.

Should I add the distance from the home to the road?

Yes. It was 5.6% of the total distance for the median resident and over 15% for one in ten, concentrated in rural areas.

Which direction should the search run?

From homes to facilities for trips residents make, which means searching from the facilities on the reversed graph. For services that travel out to homes, search the graph as it is.