Location Analysis Explained: Catchments, Accessibility and Site Selection

On this page

Problem statement

Location analysis is the part of GIS that answers practical questions about places people go to: which residents a pharmacy serves, who lives too far from one, where a new one could go, and where several should go so the whole county is served. Each question sounds like the others, and each has its own method with its own way of being wrong.

The differences are not academic. Measured on one county — Chittenden County, Vermont, with 168,323 residents in 2,241 populated census blocks, 27 pharmacies in and around it, and the OpenStreetMap drive network:

Choosing where to add one pharmacy gave two sites 0.8 km apart when the goal was "reach the most new people within 10 minutes" or "shorten the average trip", and a third site 16.0 km away when the goal was "shorten the longest trip".
Ranking tracts by how badly served they are with three standard accessibility measures put 8 tracts in the worst fifth each time — but only 6 tracts were in the worst fifth under all three.
Straight-line distance said 33.0% of residents were within 1 km of a pharmacy; along the road network it was 14.6%.

The method you pick is a statement about what the analysis values. This guide maps the questions to the methods, so the choice is deliberate.

Quick answer

question                                  method                          output
who does this facility serve?             catchment: buffer, isochrone,   population per facility
                                          Voronoi or network assignment
how well served is each place?            accessibility: nearest,         a score per home or area
                                          cumulative, gravity, 2SFCA
how much custom will a site draw?         gravity or Huff model           a share per store
where could a facility go at all?         site suitability                a score per cell
where should k facilities go together?    location-allocation: p-median,  a set of k sites
                                          coverage, p-centre

All of them start from the same inputs: where people are, where facilities are or could be, and a cost of travelling between them. Get those right first — a network travel time is usually worth the effort over a straight line — and every method downstream inherits the improvement.

Five questions share one set of inputs, which is why the travel-cost matrix deserves the most care.

Demand, facilities, travel cost, then the question asked

1. Represent demand as small populated units

People are not spread evenly across a county, and the unit you use for them decides every distance you compute. Census blocks — here 2,241 populated blocks, each as a point inside the block — are small enough that a home's distance to a facility is close to its block's. Tracts and block groups work, but each point then stands in for thousands of people spread over a large area.

2. Get facilities from a source you can check

OpenStreetMap is the usual open source for pharmacies, supermarkets, schools and fire stations. When the public Overpass API at overpass-api.de refused connections, the pharmacies here came through the maps.mail.ru Overpass mirror instead; a Geofabrik extract is the other fallback. Mapped facilities also include duplicates — a shop recorded as both a point and a building outline — which must be removed before any method that assumes one point per facility.

3. Build a travel cost, not a straight line

Straight-line distance is fast and consistently optimistic. For the county's residents, the population-weighted median distance to the nearest pharmacy was 1,381 m in a straight line and 2,342 m along the drive network, and for 13.4% of residents the nearest pharmacy by road was a different pharmacy. On the network, the median drive to the nearest pharmacy was 2.8 minutes and the longest 22.9 minutes.

4. Answer "who does it serve?" with a catchment

A catchment either counts everyone within reach of a facility — a buffer or an isochrone — or assigns every resident to one facility, as a Voronoi cell does. The median pharmacy served 2,066 people inside a 1 km circle, 601 within a 1 km walk and 4,904 in its Voronoi cell (Example 4). Say which one a number comes from.

5. Answer "how well served?" with an accessibility measure

Accessibility scores each place, not each facility. The simplest is the time to the nearest facility; a cumulative measure counts facilities within a time; a gravity measure adds them up with a weight that decays with time. They agree broadly and disagree at the margins that matter for targeting: 16,018 residents, 9.5%, had no pharmacy within 10 minutes, yet under a gravity measure their scores still ranged across a factor of 17.

6. Answer "where could it go?" with suitability

Site suitability combines constraints and graded factors on a grid — slope, road access, nearby population — into a score. In the county, constraints removed 14.9% of the land, and summing factors in their raw units instead of rescaled 0–1 values changed 8 of the 10 best sites.

7. Answer "where should several go?" with location-allocation

Location-allocation chooses a set of sites together, because the best second site depends on where the first one is. Its objectives are distinct — minimise total travel, maximise the population within a standard, minimise the longest trip — and step 8 shows they do not agree even for a single site.

8. Write the objective down before choosing a site

Adding one pharmacy to the county, with every populated block as a candidate site:

objective                         new people within 10 min   mean time   longest trip
maximise coverage                        +4,437                3.98 min     22.9 min
minimise mean time                       +4,330                3.98 min     22.9 min
minimise the longest trip                +2,175                4.13 min     21.5 min

Coverage and mean time pointed to neighbouring blocks 0.8 km apart. Minimising the longest trip moved the site 16 km to the county's edge, where it helped fewer people and shortened the worst journey by 1.4 minutes. Each answer is optimal; they optimise different things.

Two efficient objectives agreed; the equity objective chose a site 16 km away.

Accessibility and the best extra site in Chittenden County

The data used on this page

Chittenden County, Vermont. Demand is the 2020 Census: the TIGER/Line 2020 blocks for Vermont, each populated block as a representative point with its POP20. Pharmacies are OpenStreetMap amenity=pharmacy within 10 km of the county. Travel 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; the walking network is used only for the 1 km walk catchments. One block, 229 people, whose node has no road route to any pharmacy is dropped. 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 networkx as nx
import numpy as np
import osmnx as ox
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-location-analysis-explained-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
homes = vt_blocks[(vt_blocks.COUNTYFP20 == "007") & (vt_blocks.POP20 > 0)].copy()
homes["geometry"] = homes.geometry.representative_point()
homes = homes.reset_index(drop=True)

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) <= 10_000].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, L = {}, {}
for u, v, d in G.edges(data=True):                      # shortest of any parallel edges
    k = (index[u], index[v])
    secs[k], L[k] = min(secs.get(k, np.inf), d["travel_time"]), min(L.get(k, np.inf), d["length"])
rows, cols = zip(*secs)
shape = (len(index), len(index))
R = sp.csr_array((list(secs.values()), (rows, cols)), shape=shape).T      # reversed: times TO a node
L = sp.coo_array((list(L.values()), (rows, cols)), shape=shape)
p_idx = [index[n] for n in ox.distance.nearest_nodes(G, pharmacies.geometry.x, pharmacies.geometry.y)]
h_all = np.array([index[n] for n in ox.distance.nearest_nodes(G, homes.geometry.x, homes.geometry.y)])
T_all = dijkstra(R, directed=True, indices=p_idx)[:, h_all].T / 60       # home x pharmacy, minutes

ok = np.isfinite(T_all.min(axis=1))
print(f"dropping {(~ok).sum()} block ({homes.POP20[~ok].sum():,} people) with no road route to a pharmacy")
homes, T, h_idx = homes[ok].reset_index(drop=True), T_all[ok], list(h_all[ok])
w = homes["POP20"].to_numpy(float)
W = ox.project_graph(ox.graph_from_polygon(reach, network_type="walk"), to_crs=CRS)
print(f"{len(homes):,} blocks, {w.sum():,.0f} people; {len(pharmacies)} pharmacies; "
      f"drive graph {len(G):,} nodes, walk graph {len(W):,} nodes")

dropping 1 block (229 people) with no road route to a pharmacy
2,240 blocks, 168,094 people; 27 pharmacies; drive graph 10,055 nodes, walk graph 51,202 nodes

Example 1 — three accessibility measures from one travel-time matrix

import numpy as np
import pandas as pd


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


def access_measures(T, w, cutoffs=(5, 10, 15), half_lives=(2.5, 5, 10)):
    """Nearest, cumulative and gravity accessibility from an origins x facilities matrix of minutes."""
    nearest = T.min(axis=1)
    out = {"nearest": nearest}
    for c in cutoffs:
        out[f"within_{c}"] = (T <= c).sum(axis=1)
    for h in half_lives:
        out[f"gravity_{h}"] = np.exp(-np.log(2) / h * T).sum(axis=1)
    none = out["within_10"] == 0
    print(f"nearest: median {wquantile(nearest, w, .5):.1f} min, p90 {wquantile(nearest, w, .9):.1f}, "
          f"max {nearest.max():.1f}")
    print(f"nothing within 10 min: {w[none].sum():,.0f} people ({w[none].sum() / w.sum():.1%})")
    return pd.DataFrame(out)


measures = access_measures(T, w)

nearest: median 2.8 min, p90 9.8, max 22.9
nothing within 10 min: 16,018 people (9.5%)

T has one row per home and one column per pharmacy, in minutes of drive time; w is the population of each home. A gravity score with a 5-minute half-life counts a pharmacy 5 minutes away as half a pharmacy.

Example 2 — the best extra site under three objectives

def best_new_site(nearest, C, w, ids, threshold=10):
    """Best candidate for one extra facility under coverage, mean-time and worst-time objectives."""
    new = np.minimum(nearest[:, None], C)                       # homes x candidates
    gain = (w[:, None] * ((new <= threshold) & (nearest[:, None] > threshold))).sum(axis=0)
    mean_time = (w[:, None] * new).sum(axis=0) / w.sum()
    worst = new.max(axis=0)
    picks = {"coverage": int(np.argmax(gain)), "mean time": int(np.argmin(mean_time)),
             "worst time": int(np.argmin(worst))}
    for name, k in picks.items():
        print(f"{name:10} {ids[k]}: +{gain[k]:,.0f} within {threshold} min, "
              f"mean {mean_time[k]:.2f} min, worst {worst[k]:.1f} min")
    return picks


nearest = T.min(axis=1)
start = time.perf_counter()
C = dijkstra(R, directed=True, indices=h_idx)[:, h_idx].T / 60          # home x candidate block, minutes
print(f"candidate matrix {C.shape[0]:,} x {C.shape[1]:,} in {time.perf_counter() - start:.1f} s")
picks = best_new_site(nearest, C, w, homes["GEOID20"].to_numpy())
xy = np.column_stack([homes.geometry.x, homes.geometry.y])
print(f"coverage site to mean-time site: {np.linalg.norm(xy[picks['coverage']] - xy[picks['mean time']]) / 1000:.1f} km; "
      f"to worst-time site: {np.linalg.norm(xy[picks['coverage']] - xy[picks['worst time']]) / 1000:.1f} km")

candidate matrix 2,240 x 2,240 in 1.5 s
coverage   500070028004011: +4,437 within 10 min, mean 3.98 min, worst 22.9 min
mean time  500070028004019: +4,330 within 10 min, mean 3.98 min, worst 22.9 min
worst time 500070030001049: +2,175 within 10 min, mean 4.13 min, worst 21.5 min
coverage site to mean-time site: 0.8 km; to worst-time site: 16.0 km

C is the travel time from every home to every candidate site — here the 2,240 reachable populated blocks themselves, a 2,240 × 2,240 matrix that took about 1.5 s to compute. Evaluating all candidates is exhaustive for one site; for several sites at once the combinations explode, which is what location-allocation solvers are for.

Example 3 — do the measures agree on who is worst served?

def worst_fifth_agreement(areas):
    """Areas in the worst-served fifth under each measure, and how many agree."""
    q = pd.DataFrame({
        "nearest": pd.qcut(-areas["nearest"], 5, labels=False),
        "within_10": pd.qcut(areas["within_10"].rank(method="first"), 5, labels=False),
        "gravity_5": pd.qcut(areas["gravity_5"], 5, labels=False),
    })
    worst = {k: set(q.index[q[k] == 0]) for k in q}
    print(f"worst fifth per measure: {[len(v) for v in worst.values()]}; "
          f"under all three: {len(set.intersection(*worst.values()))}; "
          f"under at least one: {len(set.union(*worst.values()))} of {len(areas)}")
    return worst


tract = homes["GEOID20"].str[:11].to_numpy()
frame = measures[["nearest", "within_10", "gravity_5"]].mul(w, axis=0).assign(w=w, tract=tract)
areas = frame.groupby("tract").sum()
areas = areas[["nearest", "within_10", "gravity_5"]].div(areas["w"], axis=0)
worst = worst_fifth_agreement(areas)

worst fifth per measure: [8, 8, 8]; under all three: 6; under at least one: 11 of 40

areas holds population-weighted means of each measure per census tract. A programme targeted at the worst-served fifth reaches a different set of tracts depending on which measure was chosen, so publish the measure with the list.

Example 4 — the other figures quoted above

import shapely

# step 3: straight line against road, from the nearest-facility guide's method
d_line = shapely.distance(shapely.points(xy)[:, None], np.array(list(pharmacies.geometry))[None, :]).min(axis=1)
# road distance plus the walk from each home to the network and from the network to the pharmacy
_, h_snap = ox.distance.nearest_nodes(G, homes.geometry.x, homes.geometry.y, return_dist=True)
_, p_snap = ox.distance.nearest_nodes(G, pharmacies.geometry.x, pharmacies.geometry.y, return_dist=True)
snap_at = {}
for node, d in zip(p_idx, p_snap):
    snap_at[node] = min(d, snap_at.get(node, np.inf))
road, _, source = dijkstra(sp.csr_array(L).T, directed=True, indices=p_idx, min_only=True,
                           return_predecessors=True)
D = road[h_idx] + h_snap + np.array([snap_at[s] for s in source[h_idx]])
print(f"nearest pharmacy, population-weighted median: straight {wquantile(d_line, w, .5):,.0f} m, "
      f"road {wquantile(D, w, .5):,.0f} m; within 1 km: {w[d_line <= 1000].sum() / w.sum():.1%} "
      f"against {w[D <= 1000].sum() / w.sum():.1%}")

# step 4: three catchments per pharmacy
circle = [w[shapely.distance(shapely.points(xy), g) <= 1000].sum() for g in pharmacies.geometry]
walk_idx = {n: i for i, n in enumerate(W.nodes)}
home_walk = ox.distance.nearest_nodes(W, homes.geometry.x, homes.geometry.y)
walked = []
for node in ox.distance.nearest_nodes(W, pharmacies.geometry.x, pharmacies.geometry.y):
    reached = nx.single_source_dijkstra_path_length(W, node, cutoff=1000, weight="length")
    walked.append(w[np.isin(home_walk, list(reached))].sum())
cells = shapely.voronoi_polygons(shapely.MultiPoint(list(pharmacies.geometry)), extend_to=county.buffer(20_000),
                                 ordered=True)                    # cell i belongs to pharmacy i
voronoi = [w[shapely.contains(cell, shapely.points(xy))].sum() for cell in cells.geoms]
print(f"median pharmacy serves {np.median(circle):,.0f} in a 1 km circle, {np.median(walked):,.0f} within a "
      f"1 km walk, {np.median(voronoi):,.0f} in its Voronoi cell")

# step 5: the spread of gravity scores among people with nothing within 10 minutes
none = measures["within_10"].to_numpy() == 0
g = measures["gravity_5"].to_numpy()[none]
print(f"gravity scores with no pharmacy within 10 min: {g.min():.4f} to {g.max():.4f} ({g.max() / g.min():.0f}x)")

nearest pharmacy, population-weighted median: straight 1,381 m, road 2,342 m; within 1 km: 33.0% against 14.6%
median pharmacy serves 2,066 in a 1 km circle, 601 within a 1 km walk, 4,904 in its Voronoi cell
gravity scores with no pharmacy within 10 min: 0.1323 to 2.3059 (17x)

Everything depends on the cost matrix

Why everything depends on the cost matrix

Every method in this family is arithmetic on distances between demand and facilities. An error in those distances — straight lines where roads bend, a snap to the wrong side of a river, missing facilities across a boundary — flows unchanged into catchments, scores and chosen sites. The 2,240-candidate matrix took under two seconds. The expensive part is building and checking the network, not using it.

Why the objectives disagree

Minimising mean travel time rewards serving many people a little better, so it favours dense areas. Maximising coverage rewards bringing people inside a threshold, which also favours populated edges of existing catchments. Minimising the longest trip cares only about the single worst-served home, however few people live there. The first two are efficiency objectives and usually agree; the third is an equity objective and usually does not.

Why accessibility measures disagree at the margins

Nearest-facility time ignores choice: a home with one pharmacy 3 minutes away scores the same as one with ten. Cumulative counts ignore distance within the cutoff and jump at its edge. Gravity scores are smooth in both but depend on a half-life that someone has to choose. Across blocks the cumulative and gravity measures correlated at 0.965, while nearest time and cumulative correlated at 0.834, and the worst-served lists diverged accordingly.

Why straight lines are optimistic in a specific way

Roads rarely run directly from a home to a facility, so network distance is longer than straight-line distance almost everywhere. The error is not uniform: it is largest where water, ridges or sparse rural roads force detours, which is where people are already worst served. A straight-line analysis therefore understates exactly the gaps it is usually commissioned to find.

The gap between the two is widest at short distances, where most service standards are set.

Demand and facilities beyond the boundary

Demand outside the study area still counts. Residents beyond a boundary use facilities inside it, and vice versa.
Facilities outside the boundary matter too. For 3.4% of the county's residents, the nearest pharmacy by road was outside the county.
Free-flow speeds understate travel times at busy hours; state the speed model.
One unreachable home can break a solver. Remove or fix homes whose snapped node is disconnected before building a matrix.
Population is not demand. Pharmacy use depends on age and health; weight by the relevant population where you can.
OSM completeness varies by facility type. Check facility lists against an official register before publishing.
Results are tied to one year of data. Population, roads and facilities all change; record the vintages.

Internal links

Catchment areas explained: buffers, isochrones and Voronoi compared — who a facility serves
Accessibility measures explained: nearest, cumulative and gravity — how well each place is served
Gravity and Huff models explained: predicting where people go — market share between competing sites
Site suitability explained: constraints, factors and weights — where a facility could go
Location-allocation explained: p-median, coverage and what each optimises — where several facilities should go
How to measure distance to the nearest facility for every home — building the cost that everything uses
How to find the areas nobody can reach within a travel time — coverage gaps on a network
Street networks as graphs explained — the network underneath the costs

FAQ

What is location analysis in GIS?

The analysis of where facilities are, who they serve and where new ones should go. It covers catchments, accessibility, market share, site suitability and location-allocation, all built on travel costs between people and facilities.

Which method answers which question?

Catchments say who a facility serves; accessibility says how well each place is served; gravity models estimate market share; suitability says where a site could go; location-allocation chooses several sites together.

Do I need a road network, or is straight-line distance enough?

Use the network if you can. Straight-line distance put 33.0% of residents within 1 km of a pharmacy; the road network put 14.6% there.

Why do different optimisation objectives choose different sites?

They value different things. For one extra pharmacy, maximising coverage and minimising mean time chose sites 0.8 km apart, while minimising the longest trip chose a site 16 km away.

How do I identify the worst-served areas?

Compute more than one accessibility measure and compare. Three standard measures each flagged 8 of 40 tracts as the worst fifth, but agreed on only 6.

What data do I need to start?

Population in small units such as census blocks, a facility list you have checked, and a road network with travel times. Everything else is computed from those three.