How to Find the Areas Nobody Can Reach Within a Travel Time
Problem statement
A service promises a response within a travel time — a fire engine in eight minutes, a clinic within twenty, a depot delivery in an hour — and you need to know who is outside that promise, how many of them there are, and where they live.
The usual first attempt is a buffer around each facility. It produces a number that looks reasonable and a map that is wrong. Measured on Chittenden County, Vermont, a 5 km straight-line buffer around 48 fire stations said 4.4% of the population was uncovered; driving times on the road network said 4.1%. The totals almost agree. The people do not: 2,815 residents the buffer marked as covered are more than eight minutes away by road, and 3,335 it marked as uncovered are inside eight minutes.
This guide computes network travel time from every facility to every census block, turns thresholds into population and area shares, and extracts the gaps as ranked patches — all of it in 0.3 seconds on a county of 2,749 blocks and a 9,799-node drive network.
Quick answer
Run one multi-source shortest-path search from all facilities at once, read off the time at each block's nearest node, and threshold:
import numpy as np
from scipy.sparse.csgraph import dijkstra
from scipy.spatial import cKDTree
tree = cKDTree(node_xy) # network nodes, projected metres
_, station_nodes = tree.query(stations_xy)
snap_m, block_nodes = tree.query(block_points_xy)
seconds = dijkstra(travel_time_csr, directed=True, # facility -> everywhere
indices=np.unique(station_nodes), min_only=True)
minutes = seconds[block_nodes] / 60
beyond = minutes > 8
print(f"{pop[beyond].sum() / pop.sum():.1%} of people, "
f"{land[beyond].sum() / land.sum():.1%} of land beyond 8 minutes")Measured: 4.1% of the county's 168,323 people — 6,941 residents — and 25.5% of its land were more than eight minutes' drive from the nearest station. The search itself took 0.003 s, because min_only=True solves one problem from all 48 stations instead of 48 separate ones.
Step-by-step solution
1. Build a travel-time graph once
Load a drive network, add speeds and travel times, keep the largest strongly connected component, and store it as a sparse matrix of seconds between nodes:
import osmnx as ox
import networkx as nx
import numpy as np
import scipy.sparse as sp
G = ox.project_graph(ox.load_graphml("drive.graphml"), to_crs=32145)
G = ox.add_edge_travel_times(ox.add_edge_speeds(G))
G = G.subgraph(max(nx.strongly_connected_components(G), key=len)).copy()
nodes = list(G.nodes)
index = {n: i for i, n in enumerate(nodes)}
best = {}
for u, v, d in G.edges(data=True):
key = (index[u], index[v])
best[key] = min(best.get(key, np.inf), d["travel_time"])
rows, cols = zip(*best)
csr = sp.csr_matrix((list(best.values()), (rows, cols)), shape=(len(nodes), len(nodes)))The county plus a 10 km margin gave 9,799 nodes and 22,563 one-way arcs after dropping 36 nodes outside the main component. Parallel edges must be reduced to their minimum, not summed — csr_matrix adds duplicates by default.
Check the speeds you have been given. OSMnx imputes missing maxspeed from the mean for each road type, and here the median edge speed came out at exactly 40.6 km/h, the imputed value for the most common class. The travel times are a model, and the thresholds inherit its assumptions.
2. Include the facilities just outside your boundary
People near a county line are served across it. Measured with only the 29 stations inside the county, 1.8% of people and 12.3% of land were beyond ten minutes. Adding the 19 stations within 10 km outside the county brought that to 1.5% and 10.6%. Buffer the facility search, then report on the study area.
3. Search from all facilities at once
The question is "how far is the nearest facility", so one multi-source search answers it:
seconds = dijkstra(csr, directed=True, indices=np.unique(station_nodes), min_only=True)Direction matters on a directed graph. For an emergency response the trip runs from the station to the home, which is csr as it is. For people travelling to a service, search on the transpose, csr.T.
4. Snap blocks to the network and decide about the last leg
A block's representative point is rarely on a road. Snapping it to the nearest node gives a network time, but the walk or drive to that node is not in it:
minutes_net = seconds[block_nodes] / 60
minutes_leg = (seconds[block_nodes] + snap_m / (30 / 3.6)) / 60 # last leg at 30 km/hMeasured, adding the last leg raised the population beyond eight minutes from 4.1% to 7.0%. That is a large change for a detail, and it is not random: the blocks in the gaps had a median snapping distance of 436 m, against 143 m for all blocks. The places that are far from a station are also the places far from a road.
5. Report population and area at several thresholds
threshold people beyond share land beyond
4 min 52,734 31.3% 73.2%
6 min 19,118 11.4% 51.3%
8 min 6,941 4.1% 25.5%
10 min 2,514 1.5% 10.6%
12 min 1,050 0.6% 3.3%
15 min 89 0.1% 0.0%
Two things stand out. Area and population tell different stories — at four minutes, 73% of the land but 31% of the people. And the curve is steep: going from eight to ten minutes cuts the uncovered population by almost two-thirds. A single threshold hides both.
6. Compare with a straight-line buffer, once, to see what it gets wrong
from scipy.spatial import cKDTree
straight_m, _ = cKDTree(stations_xy).query(block_points_xy)
buffer_says_uncovered = straight_m > 5000
network_says_uncovered = minutes_net > 8The cross-tabulation of blocks:
buffer: covered buffer: uncovered
network: ≤ 8 min 2,406 121
network: > 8 min 72 150
Of the 343 blocks at least one method calls uncovered, only 150 are uncovered by both. The buffer is not a conservative approximation; it is a different map.
7. Turn gap blocks into patches and rank them
import geopandas as gpd
gaps = blocks[network_says_uncovered]
merged = gaps.geometry.buffer(1).union_all()
patches = gpd.GeoDataFrame(geometry=list(getattr(merged, "geoms", [merged])), crs=blocks.crs)
patches["people"] = [gaps.loc[gaps.intersects(p), "POP20"].sum() for p in patches.geometry]
patches = patches.sort_values("people", ascending=False)At eight minutes the county's gaps formed 34 patches. The five most populated held 50.6% of everyone beyond eight minutes — the largest had 1,083 people over 48 km² — and 5 patches had nobody living in them at all. That ranking, not the county-wide percentage, is what a new station would be sited against.
Code examples
Example 1 — nearest-facility travel time for every block
import numpy as np
from scipy.sparse.csgraph import dijkstra
from scipy.spatial import cKDTree
def travel_time_to_nearest(csr, node_xy, facility_xy, point_xy,
direction="from_facility", last_leg_kmh=None):
"""Minutes from the nearest facility to each point (or to it, with direction='to_facility')."""
tree = cKDTree(node_xy)
_, fac_nodes = tree.query(facility_xy)
snap_m, pt_nodes = tree.query(point_xy)
graph = csr if direction == "from_facility" else csr.T.tocsr()
seconds = dijkstra(graph, directed=True, indices=np.unique(fac_nodes), min_only=True)
t = seconds[pt_nodes]
if last_leg_kmh:
t = t + snap_m / (last_leg_kmh / 3.6)
unreachable = np.isinf(t).sum()
if unreachable:
print(f"warning: {unreachable} points cannot be reached — check graph connectivity")
return t / 60, snap_mExample 2 — a coverage table at several thresholds
import pandas as pd
def coverage_table(minutes, population, land_area, thresholds=(4, 6, 8, 10, 12, 15)):
"""People and land beyond each travel-time threshold."""
rows = []
for thr in thresholds:
beyond = minutes > thr
rows.append({
"threshold_min": thr,
"people_beyond": int(population[beyond].sum()),
"share_people": population[beyond].sum() / population.sum(),
"share_land": land_area[beyond].sum() / land_area.sum(),
})
table = pd.DataFrame(rows)
print(table.to_string(index=False, formatters={
"share_people": "{:.1%}".format, "share_land": "{:.1%}".format}))
return tableRun on the county's blocks with network minutes:
threshold_min people_beyond share_people share_land
4 52734 31.3% 73.2%
6 19118 11.4% 51.3%
8 6941 4.1% 25.5%
10 2514 1.5% 10.6%
12 1050 0.6% 3.3%
15 89 0.1% 0.0%
Example 3 — gap patches, ranked by the people in them
import geopandas as gpd
def gap_patches(blocks, minutes, threshold, pop_col="POP20"):
"""Dissolve adjacent blocks beyond the threshold into patches and rank them."""
gaps = blocks.loc[minutes > threshold].copy()
if gaps.empty:
return gpd.GeoDataFrame(columns=["people", "km2", "geometry"], crs=blocks.crs)
merged = gaps.geometry.buffer(1).union_all() # 1 m closes slivers between blocks
parts = list(getattr(merged, "geoms", [merged]))
patches = gpd.GeoDataFrame(geometry=parts, crs=blocks.crs)
joined = gpd.sjoin(gaps[[pop_col, "geometry"]].set_geometry(gaps.representative_point()),
patches, predicate="within")
patches["people"] = joined.groupby("index_right")[pop_col].sum().reindex(patches.index, fill_value=0)
patches["km2"] = patches.area / 1e6
patches = patches.sort_values("people", ascending=False).reset_index(drop=True)
top5 = patches.head(5)["people"].sum() / max(patches["people"].sum(), 1)
print(f"{len(patches)} patches beyond {threshold} min; the 5 largest hold {top5:.1%} of those people; "
f"{(patches['people'] == 0).sum()} are empty")
return patchesAssigning people through representative points rather than intersects means a block touching two patches is counted once.
Explanation
Why one search is enough
Dijkstra's algorithm from a single source finds the shortest time to every node. With several sources and min_only=True, SciPy seeds all of them at distance zero and runs once, which gives each node its time from the nearest source — exactly the question a coverage analysis asks. On this network that was 0.003 s, against a full origin–destination matrix that would repeat the search for every station and then take a minimum.
Why the buffer gets the total right and the map wrong
A straight-line radius assumes the road network is equally dense in every direction. It is not: valleys channel roads, lakes and ridges force detours, and fast roads stretch reach along their length. Those errors run both ways — a village on a highway 6 km out is covered, a hamlet 3 km away up a dead-end valley road is not — so they partly cancel in the county total. At the level of the block, where the decision is made, 193 blocks changed side.
Why the last leg matters most in the gaps
The network time starts at a node, and the distance from a home to that node is invisible to it. In a town that distance is a few tens of metres. In the rural blocks that make up the gaps, the representative point can be hundreds of metres from the nearest mapped road — median 436 m here — because blocks are large and roads few. So including or excluding the last leg changes the answer most exactly where the answer matters. Decide explicitly, and state which you used.
Why area and population diverge
Census blocks are drawn around people, so rural blocks are large and nearly empty while town blocks are small and dense. The land beyond four minutes was 73% of the county; the population beyond four minutes, 31%. A map of uncovered area is dominated by forest and farmland. A ranking by uncovered people is dominated by the edges of villages. Both are legitimate; mixing them is not.
Edge cases or notes
- Travel time is not response time. A fire service adds call handling and turnout; subtract them from the standard before thresholding.
- Imputed speeds are averages. Where
maxspeedis missing, every road of a type gets the same speed; calibrate against observed times if the thresholds are contractual. - Directed graphs are asymmetric. One-way systems make from-facility and to-facility times differ; choose the direction the service actually travels.
- Unreachable nodes return
inf. Keep the largest strongly connected component, or points on isolated fragments silently count as beyond every threshold. - OSM facility coverage varies. A missing station creates a gap that does not exist; check the facility list against an official one before publishing.
- Blocks with no land are lake or river. Drop them (
ALAND20 > 0) before computing land shares. - Time of day is not modelled. Free-flow speeds understate peak travel times, so a gap at 8 minutes free-flow is wider at rush hour.
- Snapping to the nearest node can cross a river. For large rural blocks, check that the snapped node is on the same side as the homes.
Internal links
- Isochrones explained: travel time areas and what they assume — the polygon view of the same calculation
- How to generate isochrones in Python — drawing the covered area rather than scoring blocks
- How to measure network distance for many origin–destination pairs — when you need every pair, not the nearest
- How to maximise coverage with a limited number of sites — where to put the next station
- Fixing accessibility scores that are wrong near the study area edge — why the stations outside the county matter
- How to count the population within reach of each facility — the facility-side view
- How to snap points to a street network in Python — the last-leg problem in detail
- Network distance vs straight-line distance explained — why the buffer map differs
FAQ
How do I find areas outside a travel time from any facility?
Run one multi-source Dijkstra search from all facilities on a travel-time graph, read the time at each block's nearest node, and select blocks above the threshold. Measured on a county network, the search took 0.003 s.
Is a buffer good enough for a coverage analysis?
For a county total it can be close — 4.4% against 4.1% here — but it misclassifies individual areas. A 5 km buffer called 2,815 people covered who were more than eight minutes away by road.
Should I include facilities outside my study area?
Yes. With only in-county stations, 1.8% of people were beyond ten minutes; adding the 19 stations just outside the boundary brought it to 1.5%.
Does the distance from a home to the nearest road matter?
In the gaps, a lot. Adding that last leg raised the share beyond eight minutes from 4.1% to 7.0%, because uncovered blocks were a median 436 m from a road node.
Should I report area or population?
Both, separately. At four minutes 73.2% of the land but 31.3% of the people were uncovered, and ranking gaps by area sends attention to empty forest.
How do I prioritise the gaps?
Dissolve adjacent uncovered blocks into patches and rank them by population. Here the five largest of 34 patches held half of everyone beyond eight minutes.