⚙️ How-To Location & accessibility Census & demographics Joins & relationships Location & Accessibility · step 4 of 20 Updated

How to Count the Population Within Reach of Each Facility

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 212,225 residents of every populated census block in the county and 10 km around it:

  • Within a 10-minute drive, the median pharmacy reached 72,589 people.
  • Adding up the county residents reached by each pharmacy gave 1,376,616 — 9.18 times the 149,931 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.537.

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 one

population 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.

Bar chart comparing the sum of per-pharmacy population counts with the number of different county residents within 5 and 10 minutes of a pharmacy.
Per-facility counts overlap; their total measures overlap, not people.

Step-by-step solution

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 16,002 people within 1 km ranked fifth by 10-minute drive, and one with 815 people within 1 km reached 79,830 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 9.0%. 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,513 blocks × 22 pharmacies took 15 ms — and compare it with the threshold. The per-facility count is a matrix product. At 10 minutes the median pharmacy reached 72,589 people; at 5 minutes, 18,063.

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 438,052, but only 115,125 different residents — 68.4% of the county — were within 5 minutes of any. At 10 minutes the sum was 1,376,616 against 149,931 residents (89.1%), and 75.6% 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 9.4% (17,409 against 15,915). 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.537
10-min drive  vs  5 km buffer    0.904
 5-min drive  vs  1 km buffer    0.752
10-min drive  vs  5-min drive    0.819

If a decision rests on which facility serves the most people, show that it survives a second definition.

Bar chart of the median and largest per-pharmacy difference between block-point and area-weighted population counts for 1 km and 5 km buffers.
Point counts are fine for large buffers and unreliable for small ones.

Code examples

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,063; summed over facilities 438,052; different residents 115,125 (68.4%)
10 min: per facility median 72,589; summed over facilities 1,376,616; different residents 149,931 (89.1%)

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 9.0%
rank agreement, 10-min drive vs 1 km buffer: 0.537

Explanation

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.

Bar chart of rank correlations between pharmacy population counts under four pairs of reach definitions.
Small buffers and drive times disagree about which pharmacies serve the most people.

Edge cases or notes

  • 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.

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,376,616, while 149,931 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 9.0%.

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.537.

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.