Accessibility Measures Explained: Nearest, Cumulative and Gravity

On this page

Problem statement

"Which neighbourhoods have poor access to pharmacies?" sounds like a question with one answer. It has at least three, because "access" can be measured three standard ways:

Nearest: the travel time to the closest facility.
Cumulative opportunities: how many facilities lie within a set travel time.
Gravity: every facility counts, weighted down the further away it is.

Measured on the 2,240 reachable populated census blocks of Chittenden County, Vermont, with drive times to 27 pharmacies, the three measures agree broadly and part company exactly where a decision gets made. Across the county's 40 census tracts:

The rank correlation between nearest time and a gravity score was 0.906.
Yet 22 of the 40 tracts landed in a different fifth of the ranking under the two measures.
Each measure put 8 tracts in the worst-served fifth; only 6 tracts were in it under all three.

Picking a measure is picking which residents a programme reaches.

Quick answer

All three come from the same origins × facilities matrix of travel times:

import numpy as np

nearest = T.min(axis=1)                                   # minutes to the closest facility
within_10 = (T <= 10).sum(axis=1)                         # facilities within 10 minutes
gravity_5 = np.exp(-np.log(2) / 5 * T).sum(axis=1)        # half-life of 5 minutes

Use nearest when one facility is enough and people go to the closest — emergency services, a single clinic. Use cumulative when choice matters and a clear standard exists — "at least two supermarkets within 10 minutes". Use gravity when choice matters and no hard cutoff is defensible. Where facilities have limited capacity relative to demand, use a two-step floating catchment (2SFCA), which divides supply by the population competing for it.

Each measure throws away something the others keep, which is why they rank places differently.

One matrix, then nearest, cumulative and gravity

1. Build the travel-time matrix once

Every measure is a summary of one matrix: rows are homes or small areas, columns are facilities, values are travel times. For 2,241 blocks and 27 pharmacies on the county's drive network the matrix took 19 ms, running one shortest-path search per pharmacy on the reversed graph. One block could not reach any pharmacy and was dropped, which should always be counted and reported.

2. Compute the nearest-facility time

nearest = T.min(axis=1)
print(f"nearest drive time: pop-weighted median {wquantile(nearest, w, .5):.1f} min, "
      f"p90 {wquantile(nearest, w, .9):.1f}, max {nearest.max():.1f}")

nearest drive time: pop-weighted median 2.8 min, p90 9.8, max 22.9

It is the easiest measure to explain and the easiest to map. It ignores everything except the closest facility: a home with one pharmacy three minutes away scores the same as a home with ten.

3. Compute cumulative opportunities at a cutoff you can defend

from scipy.stats import spearmanr

within = {c: (T <= c).sum(axis=1) for c in (5, 10, 15)}
for c, n in within.items():
    print(f"within {c:>2} min: pop-weighted median {wquantile(n, w, .5):>2}, max {n.max()}, "
          f"people with 0: {w[n == 0].sum() / w.sum():5.1%}, distinct values {len(np.unique(n))}")
zero = within[10] == 0
print(f"no pharmacy within 10 min: {w[zero].sum():,.0f} people, nearest {nearest[zero].min():.1f} "
      f"to {nearest[zero].max():.1f} min; {(within[10] == within[10].max()).sum()} blocks tied at the top")
print(f"Spearman, within 5 ~ within 15: {spearmanr(within[5], within[15])[0]:.3f}")

within  5 min: pop-weighted median  2, max 11, people with 0: 31.5%, distinct values 11
within 10 min: pop-weighted median  7, max 18, people with 0:  9.5%, distinct values 19
within 15 min: pop-weighted median 16, max 21, people with 0:  3.0%, distinct values 22
no pharmacy within 10 min: 16,018 people, nearest 10.0 to 22.9 min; 43 blocks tied at the top
Spearman, within 5 ~ within 15: 0.810

A cumulative measure captures choice, but it is coarse. At 10 minutes it took only 19 distinct values across 2,240 blocks, and 43 blocks tied at the top value of 18. And it is blind below the cutoff: the 16,018 residents with no pharmacy within 10 minutes all score 0, although their nearest pharmacy ranged from 10.0 to 22.9 minutes away.

The cutoff also changes the map. The rank correlation between counts within 5 and within 15 minutes was 0.810.

4. Compute a gravity score with a stated decay

A negative exponential with a half-life is the easiest decay to explain: a facility one half-life away counts as half, two half-lives as a quarter.

gravity = {h: np.exp(-np.log(2) / h * T).sum(axis=1) for h in (2.5, 5, 10)}
for h, g in gravity.items():
    print(f"half-life {h:>3} min: median {wquantile(g, w, .5):.2f}, max {g.max():5.2f}, "
          f"distinct values {len(np.unique(g.round(3))):,}")   # to the thousandth: smaller gaps come and go between runs
print(f"no pharmacy within 10 min: gravity(5) from {gravity[5][zero].min():.3f} to {gravity[5][zero].max():.3f}; "
      f"Spearman, half-lives 2.5 and 10: {spearmanr(gravity[2.5], gravity[10])[0]:.2f}")

half-life 2.5 min: median 1.72, max  5.48, distinct values 1,399
half-life   5 min: median 5.15, max  9.96, distinct values 1,719
half-life  10 min: median 10.55, max 14.78, distinct values 1,741
no pharmacy within 10 min: gravity(5) from 0.132 to 2.306; Spearman, half-lives 2.5 and 10: 0.98

Gravity scores separate almost every block, including the ones cumulative measures lump together: among residents with no pharmacy within 10 minutes, the 5-minute gravity score still ranged from 0.132 to 2.306. The half-life matters less than you might fear here — scores with 2.5- and 10-minute half-lives had a rank correlation of 0.98 — but it has to be chosen and reported.

5. Compare the rankings before choosing

print("Spearman rank correlation across blocks")
for label, a, b in (("nearest ~ within 10", -nearest, within[10]), ("nearest ~ gravity 5", -nearest, gravity[5]),
                    ("within 10 ~ gravity 5", within[10], gravity[5])):
    print(f"  {label:22} {spearmanr(a, b)[0]:.3f}")

Spearman rank correlation across blocks
  nearest ~ within 10    0.833
  nearest ~ gravity 5    0.858
  within 10 ~ gravity 5  0.965

Cumulative and gravity measures agree with each other much more than either agrees with nearest time, because both count several facilities.

6. Aggregate to the areas you will report

Programmes are usually targeted at tracts, not blocks. Take population-weighted means within each tract and compare quintiles:

import pandas as pd

tract = homes["GEOID20"].str[:11].to_numpy()
per_block = pd.DataFrame({"nearest": nearest, "within_10": within[10], "gravity_5": gravity[5]})
sums = per_block.mul(w, axis=0).assign(w=w, tract=tract).groupby("tract").sum()
tract_nearest, tract_within_10, tract_gravity_5 = (sums[c] / sums["w"] for c in ("nearest", "within_10", "gravity_5"))
fifth = {k: pd.qcut(v.rank(method="first"), 5, labels=False)
         for k, v in (("nearest", -tract_nearest), ("within 10", tract_within_10), ("gravity 5", tract_gravity_5))}
values = {"nearest": -tract_nearest, "within 10": tract_within_10, "gravity 5": tract_gravity_5}
for a, b in (("nearest", "within 10"), ("nearest", "gravity 5"), ("within 10", "gravity 5")):
    print(f"tract Spearman {a} ~ {b}: {spearmanr(values[a], values[b])[0]:.3f}; "
          f"tracts in a different fifth: {(fifth[a] != fifth[b]).sum()} of {len(sums)}")
worst = [set(f.index[f == 0]) for f in fifth.values()]
print(f"worst fifth: {[len(x) for x in worst]} tracts; under all three: {len(set.intersection(*worst))}")

tract Spearman nearest ~ within 10: 0.897; tracts in a different fifth: 19 of 40
tract Spearman nearest ~ gravity 5: 0.906; tracts in a different fifth: 22 of 40
tract Spearman within 10 ~ gravity 5: 0.971; tracts in a different fifth: 6 of 40
worst fifth: [8, 8, 8] tracts; under all three: 6

A correlation around 0.9 sounds like agreement. Half the tracts still changed class.

7. Choose the measure from the service, and say so

State the measure, the travel mode and speeds, the cutoff or half-life, and the facility list with its date. A published "access score" without those is unrepeatable.

High correlations still move half the tracts across a class boundary.

Cumulative counts, half-lives and how much the choice matters

The data used on this page

Chittenden County, Vermont. Demand is the 2020 Census: 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. 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 9 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 scipy.sparse as sp
import shapely
from scipy.sparse.csgraph import dijkstra

# OpenStreetMap as Overpass returned it on 9 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://www.spatialworkflow.io/data/osm-cache-accessibility-measures-explained-2026-10-09.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 = {}
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)]
h_all = np.array([index[n] for n in ox.distance.nearest_nodes(G, homes.geometry.x, homes.geometry.y)])
start = time.perf_counter()
T_all = dijkstra(R, directed=True, indices=p_idx)[:, h_all].T / 60       # home x pharmacy, minutes
print(f"matrix {T_all.shape[0]:,} x {T_all.shape[1]} in {(time.perf_counter() - start) * 1000:.0f} ms")

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 = homes[ok].reset_index(drop=True), T_all[ok]
w = homes["POP20"].to_numpy(float)


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

matrix 2,241 x 27 in 16 ms
dropping 1 block (229 people) with no road route to a pharmacy

Example 1 — cumulative counts, with the ties they create

import numpy as np


def cumulative_access(T, w, cutoffs=(5, 10, 15)):
    """Facilities within each cutoff, and how coarse each count is."""
    out = {}
    for c in cutoffs:
        counts = (T <= c).sum(axis=1)
        none = w[counts == 0].sum() / w.sum()
        top_ties = int((counts == counts.max()).sum())
        print(f"within {c:>2} min: {len(np.unique(counts))} distinct values, "
              f"{none:.1%} of people with none, {top_ties} areas tied at the top ({counts.max()})")
        out[c] = counts
    return out


counts = cumulative_access(T, w)

within  5 min: 11 distinct values, 31.5% of people with none, 3 areas tied at the top (11)
within 10 min: 19 distinct values, 9.5% of people with none, 43 areas tied at the top (18)
within 15 min: 22 distinct values, 3.0% of people with none, 33 areas tied at the top (21)

Few distinct values and large ties mean a quantile classification will split tied areas arbitrarily. That is a property of the measure, not of the classification.

Example 2 — gravity scores with a half-life

def gravity_access(T, half_life):
    """Sum of facilities weighted by exp(-ln 2 * t / half_life); a facility one half-life away counts 0.5."""
    return np.exp(-np.log(2) / half_life * T).sum(axis=1)

scores = {h: gravity_access(T, h) for h in (2.5, 5, 10)}
zero = (T <= 10).sum(axis=1) == 0
print(f"no pharmacy within 10 min: gravity(5) from {scores[5][zero].min():.3f} to {scores[5][zero].max():.3f}")

no pharmacy within 10 min: gravity(5) from 0.132 to 2.306

Choosing a half-life is choosing a behaviour: 2.5 minutes says people barely consider a second pharmacy; 10 minutes says a pharmacy across town is nearly as good as the one down the road.

Example 3 — how much the choice of measure moves a tract

import pandas as pd
from scipy.stats import spearmanr


def measure_agreement(areas, pairs):
    """Rank correlation and fifth-of-ranking changes between pairs of measures (higher = better access)."""
    fifths = {m: pd.qcut(areas[m].rank(method="first"), 5, labels=False) for m in areas}
    for a, b in pairs:
        rho = spearmanr(areas[a], areas[b])[0]
        moved = int((fifths[a] != fifths[b]).sum())
        print(f"{a} ~ {b}: rho {rho:.3f}; {moved} of {len(areas)} areas in a different fifth")

tracts = pd.DataFrame({"nearest": -tract_nearest, "within_10": tract_within_10, "gravity_5": tract_gravity_5})
measure_agreement(tracts, [("nearest", "within_10"), ("nearest", "gravity_5"), ("within_10", "gravity_5")])

nearest ~ within_10: rho 0.897; 19 of 40 areas in a different fifth
nearest ~ gravity_5: rho 0.906; 22 of 40 areas in a different fifth
within_10 ~ gravity_5: rho 0.971; 6 of 40 areas in a different fifth

Negate nearest time so that every column reads "higher is better". Run it before any list of priority areas is published.

What each measure throws away

Nearest time keeps distance and discards choice. Cumulative counts keep choice and discard distance, both inside the cutoff and beyond it. Gravity keeps both, at the price of a decay parameter that has no natural value. None is neutral: each assumes something about how people choose where to go.

Why cumulative measures clump

A count of facilities within a cutoff can only take whole-number values up to the number of facilities in reach. In a county with 27 pharmacies, most blocks share a handful of values, and whole towns tie. Gravity scores are continuous, so almost every block gets its own value — 1,719 distinct scores to the thousandth with a 5-minute half-life, against 19 distinct counts.

Why block correlation overstates tract agreement

A rank correlation weighs every pair of areas equally, and most pairs are far apart in any sensible ranking, so the correlation stays high. Classification cares about the areas near each class boundary. With 40 tracts in five classes, a small reordering near each boundary moves many tracts, which is why a correlation of 0.906 coexisted with 22 class changes.

When none of these is enough

All three measures treat a facility as equally available to everyone who can reach it. A pharmacy in a dense neighbourhood serves more people than one in a village, so the same travel time buys less service. Measures that divide supply by the population competing for it — the two-step floating catchment method — address that, and they disagree with the simpler measures most in dense areas with few facilities.

The service decides the measure: whether people need one facility, a choice, or a share of limited capacity.

Unreachable areas, cliffs and half-lives

Report unreachable areas separately. One block could not reach any pharmacy; leaving it as infinity breaks means and rankings.
Cutoffs create cliffs. A home at 10.1 minutes scores like a home at 30 minutes under a 10-minute cumulative measure.
Half-lives are not universal. Driving, walking and public transport need different decays.
Facilities outside the study area count. Omitting them understates access near every boundary.
Population weighting changes medians. Report both area-weighted and population-weighted summaries if areas differ in size.
Speeds are a model. Free-flow network speeds understate peak travel times, which lowers every measure at busy hours.
Gravity sums grow with the number of facilities. Compare scores only within one facility set.

Internal links

Location analysis explained: catchments, accessibility and site selection — where accessibility fits
How to measure distance to the nearest facility for every home — computing the nearest measure on a network
How to measure access to services with the two-step floating catchment method — accessibility with capacity
Catchment areas explained: buffers, isochrones and Voronoi compared — the facility-side view
Fixing accessibility scores that are wrong near the study area edge — boundary effects on every measure
Gravity and Huff models explained: predicting where people go — distance decay used for market share
How to measure network distance for many origin–destination pairs — building the travel-time matrix
Choropleth classification explained — why quintile boundaries move areas
Transit accessibility explained — why transit measures need a departure window

FAQ

What is an accessibility measure?

A score for each place describing how easily its residents can reach a type of facility. The standard forms are time to the nearest facility, the number within a cutoff, and a distance-weighted sum of all facilities.

Which accessibility measure should I use?

Match it to the service. Nearest time suits services where people use the closest one; cumulative counts suit standards such as two supermarkets within 10 minutes; gravity suits choice with no hard cutoff.

Do different measures give different results?

At the margins, yes. For 40 tracts, nearest time and a gravity score correlated at 0.906, yet 22 tracts fell in a different fifth of the ranking.

What is a gravity accessibility score?

A sum over all facilities in which each is weighted by a decay of travel time. With a half-life of 5 minutes, a pharmacy 5 minutes away counts as half and one 10 minutes away as a quarter.

Why do cumulative measures have so many ties?

They are whole-number counts. Across 2,240 blocks the count of pharmacies within 10 minutes took only 19 values, and 43 blocks tied at the maximum.

How do I choose a half-life or cutoff?

From travel behaviour or a service standard, and then test the sensitivity. Here gravity scores with 2.5- and 10-minute half-lives had a rank correlation of 0.98, while counts within 5 and 15 minutes correlated at 0.810.