How to Count the Population Within Reach of Each Facility
On this page
Problem statement
"How many people live within reach of each pharmacy?" sounds like one number per pharmacy. It is several decisions: straight line or road, which threshold, how a census block that straddles the edge is counted, whether people over the study-area boundary count, and what happens when the per-facility numbers are added up.
Measured for the 22 pharmacies in Chittenden County, Vermont, with the 209,041 residents of every populated Vermont census block in the county and 10 km around it:
- Within a 10-minute drive, the median pharmacy reached 72,314 people.
- Adding up the county residents reached by each pharmacy gave 1,371,533 β 9.16 times the 149,702 different residents actually within 10 minutes of one.
- For 1 km buffers, counting block points instead of block areas changed a pharmacy's count by a median of 10.6%, and by up to 30.6%.
- Ranking pharmacies by people within 1 km and by people within a 10-minute drive barely agreed: Spearman 0.531.
Quick answer
Count per facility from a travel-time matrix, and count unique people separately:
import numpy as np
within = minutes <= 10 # homes x facilities, travel time in minutes
per_facility = within.T @ population # people within reach of each facility
reached = population[within.any(axis=1)].sum() # people within reach of at least onepopulation must include the blocks beyond the study area that a facility near the boundary can reach. For straight-line buffers, share each block's population by area rather than testing a single point.
What you get β from The data used on this page
3,452 populated blocks (2,241 in the county), 209,041 people; 22 pharmacies; matrix in 13 ms
Define reach, include the margin, count by area
1. Decide what "within reach" means
Pick the mode and threshold from the service: a walk of 1 km, a drive of 10 minutes. Straight-line buffers are quick and reproducible; travel time follows the roads. They are not interchangeable here: the pharmacy with 15,970 people within 1 km ranked fifth by 10-minute drive, and the one with the fewest, 831 people within 1 km, reached 70,197 within 10 minutes.
2. Include the population beyond the boundary
A pharmacy near the county line serves people on both sides. Counting only county residents undercounted one pharmacy's 10-minute population by 8.6%. At 5 minutes no pharmacy was affected, because none reached across the line that quickly. Load demand for the study area plus a margin at least as wide as the reach.
3. Count buffers by area, not by point
A census block is one point in the simplest count, even when half of it lies outside the buffer. With 1 km buffers, a block point count and an area-weighted count differed by a median of 10.6% per pharmacy, and by 30.6% for one: 866 people by point, 663 by area. At 5 km the median difference fell to 0.9% and the largest to 4.3%, because few blocks are large relative to a 5 km circle.
4. Count travel-time reach from a matrix
Build the homes Γ facilities travel-time matrix once β 3,452 blocks Γ 22 pharmacies took 16 ms β and compare it with the threshold. The per-facility count is a matrix product. At 10 minutes the median pharmacy reached 72,314 people; at 5 minutes, 18,330.
5. Count unique people separately
Per-facility counts overlap wherever facilities are close together. Within 5 minutes, the sum of county residents reached by each pharmacy was 437,119, but only 114,853 different residents β 68.2% of the county β were within 5 minutes of any. At 10 minutes the sum was 1,371,533 against 149,702 residents (88.9%), and 75.2% of residents could reach two or more pharmacies. Never report the sum as a number of people.
6. Use the finest population units you have
Block groups are coarser than blocks. With 5 km buffers and county residents only, counting block-group centroids instead of block points changed a pharmacy's count by a median of 2.0% and by up to 10.0% (22,148 against 20,138). The error grows as the reach shrinks relative to the units.
7. Compare the rankings, not only the counts
The ranking of pharmacies changed with the definition. Spearman rank correlations between definitions:
10-min drive vs 1 km buffer 0.531
10-min drive vs 5 km buffer 0.912
5-min drive vs 1 km buffer 0.726
10-min drive vs 5-min drive 0.833
If a decision rests on which facility serves the most people, show that it survives a second definition.
Pharmacy reach by buffer and by drive time
The data used on this page
Chittenden County, Vermont, and 10 km around it. Population is the 2020 Census: every populated Vermont block polygon from the TIGER/Line 2020 blocks within 10 km of the county, with its POP20; New York blocks across Lake Champlain have no road link inside the margin and are left out. The facilities are the 22 OpenStreetMap amenity=pharmacy features inside the county. Drive times come from OpenStreetMap's drive network for the county plus 10.5 km, with OSMnx's imputed free-flow speeds, road-only between nearest nodes. 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 time
import urllib.request
import zipfile
import geopandas as gpd
import numpy as np
import osmnx as ox
import pandas as pd
import scipy.sparse as sp
import shapely
from scipy.sparse.csgraph import dijkstra
# 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-population-within-reach-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.POP20 > 0]
blocks = blocks[shapely.distance(blocks.representative_point().values, county) <= 10_000].reset_index(drop=True)
in_county = (blocks.COUNTYFP20 == "007").to_numpy()
population = blocks["POP20"].to_numpy(float)
reach = gpd.GeoSeries([county.buffer(10_500)], crs=CRS).to_crs(4326).iloc[0]
pharmacies = ox.features_from_polygon(reach, tags={"amenity": "pharmacy"}).to_crs(CRS)
pharmacies["geometry"] = pharmacies.geometry.representative_point()
pharmacies = pharmacies[shapely.distance(pharmacies.geometry.values, county) == 0].reset_index()
G = ox.add_edge_travel_times(ox.add_edge_speeds(ox.graph_from_polygon(reach, network_type="drive")))
G = ox.project_graph(G, to_crs=CRS)
index = {n: i for i, n in enumerate(G.nodes)}
secs = {}
for u, v, d in G.edges(data=True): # shortest of any parallel edges
secs[(index[u], index[v])] = min(secs.get((index[u], index[v]), np.inf), d["travel_time"])
rows, cols = zip(*secs)
R = sp.csr_array((list(secs.values()), (rows, cols)), shape=(len(index), len(index))).T.tocsr()
p_idx = [index[n] for n in ox.distance.nearest_nodes(G, pharmacies.geometry.x, pharmacies.geometry.y)]
pts = blocks.representative_point()
b_idx = [index[n] for n in ox.distance.nearest_nodes(G, pts.x, pts.y)]
start = time.perf_counter()
minutes = dijkstra(R, directed=True, indices=p_idx)[:, b_idx].T / 60 # blocks x pharmacies
print(f"{len(blocks):,} populated blocks ({in_county.sum():,} in the county), {population.sum():,.0f} people; "
f"{len(pharmacies)} pharmacies; matrix in {(time.perf_counter() - start) * 1000:.0f} ms")3,452 populated blocks (2,241 in the county), 209,041 people; 22 pharmacies; matrix in 13 ms
Example 1 β buffer counts by point and by area
import geopandas as gpd
import numpy as np
import pandas as pd
def buffer_population(facilities, blocks, radius, population="POP20"):
"""People within radius of each facility, counted by block point and by block area."""
rings = gpd.GeoDataFrame({"facility": np.arange(len(facilities))},
geometry=facilities.buffer(radius).to_numpy(), crs=facilities.crs)
points = gpd.GeoDataFrame(blocks[[population]], geometry=blocks.representative_point(), crs=blocks.crs)
by_point = gpd.sjoin(points, rings, predicate="within").groupby("facility")[population].sum()
parts = gpd.overlay(blocks[[population, "geometry"]].assign(block_m2=blocks.area), rings,
how="intersection", keep_geom_type=True)
by_area = (parts[population] * parts.area / parts.block_m2).groupby(parts.facility).sum()
return pd.DataFrame({"by_point": by_point, "by_area": by_area}).reindex(rings.facility, fill_value=0)
counts = buffer_population(pharmacies, blocks, 1000)
gap = (counts.by_point - counts.by_area).abs() / counts.by_area.clip(lower=1)
print(f"median by point {counts.by_point.median():,.0f}, by area {counts.by_area.median():,.0f}; "
f"difference median {gap.median():.1%}, max {gap.max():.1%}")median by point 2,528, by area 2,750; difference median 10.6%, max 30.6%
blocks holds every populated block polygon in the study area and its margin, in a projected CRS in metres.
Example 2 β drive-time reach, per facility and unique
def reach_counts(minutes, population, threshold, in_study_area):
within = minutes <= threshold
per_facility = within.T @ population
summed = (within.T @ (population * in_study_area)).sum()
reached = population[in_study_area & within.any(axis=1)].sum()
share = reached / population[in_study_area].sum()
print(f"{threshold} min: per facility median {np.median(per_facility):,.0f}; "
f"summed over facilities {summed:,.0f}; different residents {reached:,.0f} ({share:.1%})")
return per_facility
drive_5 = reach_counts(minutes, population, 5, in_county)
drive_10 = reach_counts(minutes, population, 10, in_county)5 min: per facility median 18,330; summed over facilities 437,119; different residents 114,853 (68.2%)
10 min: per facility median 72,314; summed over facilities 1,371,533; different residents 149,702 (88.9%)
minutes has one row per populated block, including the margin, and one column per pharmacy; in_county flags the rows inside the study area.
Example 3 β what the boundary and the definition cost
from scipy.stats import spearmanr
def boundary_shortfall(minutes, population, threshold, in_study_area):
"""Share of each facility's reach lost by ignoring people outside the study area."""
within = minutes <= threshold
everyone = within.T @ population
local = within.T @ (population * in_study_area)
return 1 - local / np.maximum(everyone, 1)
for threshold in (5, 10):
lost = boundary_shortfall(minutes, population, threshold, in_county)
print(f"{threshold} min: shortfall median {np.median(lost):.1%}, max {lost.max():.1%}")
print(f"rank agreement, 10-min drive vs 1 km buffer: {spearmanr(drive_10, counts.by_area)[0]:.3f}")5 min: shortfall median 0.0%, max 0.0%
10 min: shortfall median 0.0%, max 8.6%
rank agreement, 10-min drive vs 1 km buffer: 0.531
Example 4 β the other comparisons quoted above
w_c = population * in_county
within10 = minutes <= 10
print(f"county residents who can reach two or more pharmacies in 10 min: "
f"{w_c[within10.sum(axis=1) >= 2].sum() / w_c.sum():.1%}")
rank10 = pd.Series(drive_10).rank(ascending=False)
top_1km = int(np.argmax(counts.by_point))
fewest_1km = int(np.argmin(counts.by_point))
print(f"most people within 1 km: {counts.by_point.iloc[top_1km]:,.0f}, ranked {int(rank10[top_1km])} by 10-min drive; "
f"fewest: {counts.by_point.iloc[fewest_1km]:,.0f} within 1 km, {drive_10[fewest_1km]:,.0f} within 10 min")
c5 = buffer_population(pharmacies, blocks, 5000)
gap5 = (c5.by_point - c5.by_area).abs() / c5.by_area.clip(lower=1)
worst = int(np.argmax(gap))
print(f"1 km worst: {counts.by_point.iloc[worst]:,.0f} by point, {counts.by_area.iloc[worst]:,.0f} by area; "
f"5 km difference median {gap5.median():.1%}, max {gap5.max():.1%}")
# block groups against blocks, 5 km buffers, county residents only
county_blocks = blocks[in_county]
pts = county_blocks.representative_point()
bg = pd.DataFrame({"bg": county_blocks.GEOID20.str[:12].to_numpy(), "p": county_blocks.POP20.to_numpy(),
"x": pts.x.to_numpy() * county_blocks.POP20.to_numpy(),
"y": pts.y.to_numpy() * county_blocks.POP20.to_numpy()}).groupby("bg").sum()
bg_pts = gpd.points_from_xy(bg.x / bg.p, bg.y / bg.p)
by_block = np.array([county_blocks.POP20[pts.distance(g) <= 5000].sum() for g in pharmacies.geometry])
by_bg = np.array([bg.p[bg_pts.distance(g) <= 5000].sum() for g in pharmacies.geometry])
d = np.abs(by_bg - by_block) / by_block
i = int(np.argmax(d))
print(f"{len(bg)} block groups: difference median {np.median(d):.1%}, max {d.max():.1%} "
f"({by_bg[i]:,.0f} against {by_block[i]:,.0f})")
print("Spearman rank correlations")
for label, a, b in (("10-min drive vs 1 km buffer", drive_10, counts.by_area),
("10-min drive vs 5 km buffer", drive_10, c5.by_area),
(" 5-min drive vs 1 km buffer", drive_5, counts.by_area),
("10-min drive vs 5-min drive", drive_10, drive_5)):
print(f" {label} {spearmanr(a, b)[0]:.3f}")county residents who can reach two or more pharmacies in 10 min: 75.2%
most people within 1 km: 15,970, ranked 5 by 10-min drive; fewest: 831 within 1 km, 70,197 within 10 min
1 km worst: 866 by point, 663 by area; 5 km difference median 0.9%, max 4.3%
118 block groups: difference median 2.0%, max 10.0% (22,148 against 20,138)
Spearman rank correlations
10-min drive vs 1 km buffer 0.531
10-min drive vs 5 km buffer 0.912
5-min drive vs 1 km buffer 0.726
10-min drive vs 5-min drive 0.833
Why the sums exceed the population
Why the sums exceed the population
Every person within reach of three pharmacies is counted three times. The overlap is the normal state of a town with several pharmacies on the same streets, so the sum grows with facility density, not with population: at 10 minutes it was more than nine times the number of people.
Why area weighting matters for small buffers
A 1 km buffer covers 3.14 kmΒ². Rural census blocks can be larger than that, so whether a block's point falls inside decides the fate of everyone in it. Area weighting assumes people are spread evenly across the block, which is also wrong, but its errors are smaller and do not jump.
Why buffers and drive times rank facilities differently
A buffer rewards a pharmacy in a dense centre; a drive-time reach rewards one near a fast road that connects many places. A central pharmacy on slow streets can have the most people within 1 km and fewer within 10 minutes than an edge-of-town store by a main road.
What "within reach" leaves out
Reach is binary and ignores competition: a pharmacy reaching 80,000 people shares most of them with a dozen others. When supply relative to demand is the question, use the two-step floating catchment method.
Units, snapping and unreachable homes
- Threshold units must match the matrix. A matrix in seconds compared with 10 counts almost nobody.
- Snap distances are not in the matrix. For walking thresholds, add the leg from the block point to the network.
- Unreachable homes are infinite. They fail every threshold; count them separately.
- Areal weighting assumes even spread. Mask uninhabited land first if you have a land cover or building layer.
- Direction matters on one-way streets. Measure home to facility, not the reverse, for trips to a service.
- Population vintage must match the facilities. 2020 blocks against 2026 pharmacies is a stated compromise.
- Buffers across water overstate reach. A buffer over a lake counts the far shore as close.
Internal links
- Catchment areas explained: buffers, isochrones and Voronoi compared β the three catchment definitions
- How to measure access to services with the two-step floating catchment method β reach with competition
- How to measure network distance for many originβdestination pairs β the travel-time matrix
- How to generate isochrones in Python β reach as a polygon
- How to create buffers in GeoPandas for spatial analysis β buffer basics
- Fixing accessibility scores that are wrong near the study area edge β demand beyond the boundary
- How to build Voronoi service areas around facilities in Python β non-overlapping areas instead
- Census geographies explained: blocks, tracts, output areas and why they nest β choosing population units
- How to measure the population served by frequent transit β the same count for transit stops
FAQ
How do I count the population within a distance of each facility?
Buffer each facility and share every census block's population by the fraction of its area inside the buffer. For drive times, compare a homes Γ facilities travel-time matrix with the threshold and multiply by population.
Can I add up the population within reach of each facility?
Not as a count of people. For 22 pharmacies the per-pharmacy 10-minute counts summed to 1,371,533, while 149,702 different residents were within 10 minutes of any pharmacy.
Is counting block centroids inside a buffer accurate enough?
For large buffers, usually. For 1 km buffers the point count differed from an area-weighted count by a median of 10.6% per pharmacy; for 5 km buffers by 0.9%.
Do I need population outside my study area?
Yes, for facilities near the boundary. Ignoring it cut one pharmacy's 10-minute population by 8.6%.
Should I use a buffer or a drive time?
Use the one that matches how people reach the service. They ranked pharmacies differently: a 1 km buffer and a 10-minute drive had a rank correlation of only 0.531.
Does the size of the census unit matter?
Yes, and more as the reach gets smaller. Block-group centroids changed 5 km buffer counts by a median of 2.0% and up to 9.4% compared with blocks.