How to Find the Areas Nobody Can Reach Within a Travel Time

On this page

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,347 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 — the search itself in 3 milliseconds on a county of 2,749 blocks and a 10,019-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,929 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.

A gap map coloured by area shows mostly empty hills; weight by people before you rank anything.

What you get — from The data used on this page

2,749 blocks, 168,323 people; 48 fire stations; 10,019 nodes and 23,077 arcs after dropping 36 outside the main component

One graph, neighbours included, one multi-source search

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 10,019 nodes and 23,077 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.9% of people and 12.6% of land were beyond ten minutes. Adding the 19 stations within 10 km outside the county brought that to 1.5% and 10.9%. 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/h

Measured, adding the last leg raised the population beyond eight minutes from 4.1% to 6.9%. That is a large change for a detail, and it is not random: the blocks in the gaps had a median snapping distance of 443 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,047      30.9%       72.9%
  6 min        19,128      11.4%       51.3%
  8 min         7,039       4.2%       26.6%
 10 min         2,619       1.6%       11.9%
 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 > 8

The cross-tabulation of blocks:

                    buffer: covered   buffer: uncovered
network: ≤ 8 min         2,406               117
network: > 8 min            72               154

Of the 343 blocks at least one method calls uncovered, only 149 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.7% 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.

Similar totals, different people: the buffer and the network disagree about 189 blocks.

Travel time to the nearest facility for every block

The data used on this page

Chittenden County, Vermont. The blocks are every TIGER/Line 2020 block in the county with land (ALAND20 > 0), with POP20 for people and ALAND20 for land area, each as a representative point. The facilities are OpenStreetMap amenity=fire_station features in the county and within 10 km of it. The road network is OpenStreetMap's drive network for the county plus 10.5 km, restricted to its largest strongly connected component, with OSMnx's imputed free-flow travel times, held as a sparse matrix of seconds from node to node. 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 requests
import urllib.request
import zipfile

import geopandas as gpd
import networkx as nx
import numpy as np
import osmnx as ox
import scipy.sparse as sp
import shapely

# 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-find-service-coverage-gaps-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/)"
CRS = 32145                                               # NAD83 / Vermont, metres

if not os.path.exists("tl_2020_50_tabblock20.zip"):
    urllib.request.urlretrieve("https://www2.census.gov/geo/tiger/TIGER2020/TABBLOCK20/"
                               "tl_2020_50_tabblock20.zip", "tl_2020_50_tabblock20.zip")
vt_blocks = gpd.read_file("tl_2020_50_tabblock20.zip").to_crs(CRS)
county = vt_blocks[vt_blocks.COUNTYFP20 == "007"].union_all()           # Chittenden
blocks = vt_blocks[(vt_blocks.COUNTYFP20 == "007") & (vt_blocks.ALAND20 > 0)].reset_index(drop=True)
pts = blocks.representative_point()
block_points_xy = np.column_stack([pts.x, pts.y])
pop, land = blocks["POP20"].to_numpy(float), blocks["ALAND20"].to_numpy(float)

reach = gpd.GeoSeries([county.buffer(10_500)], crs=CRS).to_crs(4326).iloc[0]
stations = ox.features_from_polygon(reach, tags={"amenity": "fire_station"}).to_crs(CRS)
stations["geometry"] = stations.geometry.representative_point()
stations = stations[shapely.distance(stations.geometry.values, county) <= 10_000]
stations = stations[~stations.geometry.duplicated()].reset_index()
stations_xy = np.column_stack([stations.geometry.x, stations.geometry.y])
stations_in_county = (shapely.distance(stations.geometry.values, county) == 0)

G = ox.graph_from_polygon(reach, network_type="drive")
G = ox.project_graph(ox.add_edge_travel_times(ox.add_edge_speeds(G)), to_crs=CRS)
n_all = len(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)}
node_xy = np.array([(G.nodes[n]["x"], G.nodes[n]["y"]) for n in nodes])
best = {}
for u, v, d in G.edges(data=True):                      # the minimum of any parallel edges
    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)))
print(f"{len(blocks):,} blocks, {pop.sum():,.0f} people; {len(stations)} fire stations; "
      f"{len(nodes):,} nodes and {len(best):,} arcs after dropping {n_all - len(nodes)} outside the main component")

2,749 blocks, 168,323 people; 48 fire stations; 10,019 nodes and 23,077 arcs after dropping 36 outside the main component

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_m


minutes, snap_m = travel_time_to_nearest(csr, node_xy, stations_xy, block_points_xy)

Example 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 table


table = coverage_table(minutes, pop, land)

 threshold_min  people_beyond share_people share_land
             4          52734        31.3%      73.2%
             6          19092        11.3%      51.1%
             8           6929         4.1%      25.5%
            10           2592         1.5%      10.9%
            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 patches


patches = gap_patches(blocks, minutes, 8)
print(f"largest: {patches.people.iloc[0]:,} people over {patches.km2.iloc[0]:.0f} km²")

34 patches beyond 8 min; the 5 largest hold 50.7% of those people; 5 are empty
largest: 1,083 people over 48 km²

Assigning people through representative points rather than intersects means a block touching two patches is counted once.

Example 4 — the comparisons quoted above

import time

from scipy.spatial import cKDTree

start = time.perf_counter()
m, _ = travel_time_to_nearest(csr, node_xy, stations_xy, block_points_xy)
print(f"search and lookup: {time.perf_counter() - start:.3f} s")
beyond = m > 8
print(f"beyond 8 min: {pop[beyond].sum():,.0f} people ({pop[beyond].sum() / pop.sum():.1%}), "
      f"{land[beyond].sum() / land.sum():.1%} of land; median edge speed "
      f"{np.median([d['speed_kph'] for *_, d in G.edges(data=True)]):.1f} km/h")

inside = stations_in_county
m_in, _ = travel_time_to_nearest(csr, node_xy, stations_xy[inside], block_points_xy)
for label, mm in (("county stations only", m_in), ("with those within 10 km", m)):
    b = mm > 10
    print(f"{label}: beyond 10 min {pop[b].sum() / pop.sum():.1%} of people, {land[b].sum() / land.sum():.1%} of land")
print(f"stations: {inside.sum()} in the county, {(~inside).sum()} within 10 km outside")

leg, _ = travel_time_to_nearest(csr, node_xy, stations_xy, block_points_xy, last_leg_kmh=30)
print(f"with the last leg at 30 km/h: {pop[leg > 8].sum() / pop.sum():.1%} beyond 8 min; median snap "
      f"{np.median(snap_m[beyond]):.0f} m in the gaps, {np.median(snap_m):.0f} m overall")

straight_m, _ = cKDTree(stations_xy).query(block_points_xy)
buf_out, net_out = straight_m > 5000, m > 8
print(f"5 km buffer: {pop[buf_out].sum() / pop.sum():.1%} uncovered; network {pop[net_out].sum() / pop.sum():.1%}")
print(f"covered by buffer but over 8 min: {pop[~buf_out & net_out].sum():,.0f} people; "
      f"uncovered by buffer but within 8 min: {pop[buf_out & ~net_out].sum():,.0f}")
print(f"blocks: net<=8 & buffer covered {np.sum(~net_out & ~buf_out):,}, net<=8 & buffer uncovered "
      f"{np.sum(~net_out & buf_out)}, net>8 & buffer covered {np.sum(net_out & ~buf_out)}, both uncovered "
      f"{np.sum(net_out & buf_out)}")

search and lookup: 0.003 s
beyond 8 min: 6,929 people (4.1%), 25.5% of land; median edge speed 40.6 km/h
county stations only: beyond 10 min 1.9% of people, 12.6% of land
with those within 10 km: beyond 10 min 1.5% of people, 10.9% of land
stations: 29 in the county, 19 within 10 km outside
with the last leg at 30 km/h: 6.9% beyond 8 min; median snap 443 m in the gaps, 143 m overall
5 km buffer: 4.4% uncovered; network 4.1%
covered by buffer but over 8 min: 2,815 people; uncovered by buffer but within 8 min: 3,347
blocks: net<=8 & buffer covered 2,406, net<=8 & buffer uncovered 122, net>8 & buffer covered 72, both uncovered 149

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, 189 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.

The only slow part of the workflow is building the graph; every threshold after that is free.

Response times, imputed speeds and directed graphs

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 maxspeed is 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.9% 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 6.9%, because uncovered blocks were a median 443 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.