Fixing a Gravity Model That Sends Everyone to the Biggest Store
On this page
Problem statement
You fit a Huff model to a set of stores and map each neighbourhood's most likely store. One store — the largest — wins almost everywhere, including in towns with their own supermarket a few minutes away. Or its predicted market share looks far too large. Nothing raised an error.
Measured with 39 supermarkets around Chittenden County, Vermont, mapped store outlines as attractiveness and network drive times from 2,239 populated census blocks, the largest store is genuinely the nearest store for 9.3% of residents. Five common mistakes each inflate it:
setting homes whose likeliest store is the largest
power decay, β = 2, minutes (reference) 12.7%
no distance decay, β = 0 100.0%
β = 0.5 27.5%
exponential λ = 0.2 per minute, times in hours 59.1%
power decay on (hours + 1) 39.0%
And two inflate its market share instead: raising attractiveness to the power 2 took its share from 8.5% to 13.2%, and filling missing footprints with 1 m² took it to 11.2%.
Quick answer
Check the parameters against the units and the data, in this order:
import numpy as np
A = footprints.fillna(footprints.median()).to_numpy() # a stated fill rule, not an accident
T = np.maximum(T_minutes, 0.5) # no zero travel times
P = A[None, :] ** 1.0 * T ** -2.0 # power decay: unit-free
P /= P.sum(axis=1, keepdims=True)
largest = A.argmax()
print(f"likeliest is largest: {pop[P.argmax(axis=1) == largest].sum() / pop.sum():.1%}")
print(f"nearest is largest: {pop[T.argmin(axis=1) == largest].sum() / pop.sum():.1%}")likeliest is largest: 12.7%
nearest is largest: 9.3%
If the first figure is far above the second, distance decay is too weak or in the wrong units. If the largest store's share is far above its share of floor area, α is too high or missing attractiveness is filled too low.
Measure against nearest-store, then check decay and units
1. Measure the symptom against the nearest-store baseline
"Everyone goes to the biggest store" needs a number. Compute the share of population whose likeliest store is the largest, and compare it with the share whose nearest store is the largest — 9.3% here. A well-specified model with reasonable decay lands somewhat above the nearest-store share, because a large store does pull custom past smaller ones: 12.7% at β = 2.
2. Check that distance decay is switched on
With β = 0, travel time drops out: every home's probabilities are proportional to floor area, so the largest store is everyone's likeliest store — 100% — and its share equals its share of floor area, 6.7%. This happens when β is read from a configuration as zero, or when a variable meant to hold travel time holds something constant.
3. Check that decay is strong enough
Small exponents leave distance too weak to overcome size. At β = 0.5 the largest store was the likeliest for 27.5% of residents. Unless you have calibrated β, test the range you consider plausible and show how the map changes across it:
β = 0.5 27.5%
β = 1.0 14.5%
β = 2.0 12.7%
β = 5.0 11.3%
4. Check the units of an exponential or offset decay
The power form,
exponential λ = 0.2 per minute likeliest is largest 18.5%
same λ applied to times in seconds 9.3% (decay 60× too strong)
same λ applied to times in hours 59.1% (decay 60× too weak)
An offset such as
5. Check the attractiveness exponent
α = 1 treats a store twice the size as twice as attractive. Larger α lets large stores capture custom far beyond their surroundings:
α = 0.5 largest store's share 6.2% top three 18.7%
α = 1.0 8.5% 21.7%
α = 2.0 13.2% 32.5%
Values above 1 need evidence. Without calibration, α = 1 is the neutral choice.
6. Check how missing attractiveness was filled
Fifteen of the 39 supermarkets were mapped only as points in OpenStreetMap, with no outline. The fill rule changes who wins:
missing footprints filled with the median largest store's share 8.5%
missing footprints filled with the maximum 5.9%
missing footprints filled with 1 m² 11.2%
Filling with a tiny value — or with zero, which removes the stores entirely — hands their customers to the stores that do have footprints, including the largest. Use a real measure for every store if at all possible.
7. Floor zero travel times instead of dropping rows
Seven origin–store pairs had a travel time of zero. Without a floor, those rows became NaN in the power form; if they are later dropped or filled with zeros, the result is biased in whichever direction the rows happened to be. Floor times at a small positive value in the same unit.
8. Recalibrate if you have observed data
If loyalty-card, survey or turnover data exist, fit β and α to them. Two uncalibrated parameters can reproduce almost any market story; calibrated ones are a model.
Diagnosing and sweeping the decay parameter
The data used on this page
Chittenden County, Vermont. Origins are the 2020 Census: the TIGER/Line 2020 blocks for Vermont, each populated block as a representative point with its POP20. Stores are OpenStreetMap shop=supermarket within 10 km of the county; attractiveness is the area of the store's own mapped outline, and stores mapped only as points are filled with the median outline (step 6 shows what the fill rule does). Drive times in minutes 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. One block (229 people) with no road route to any store 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 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 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-gravity-model-favours-biggest-site-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]
stores = ox.features_from_polygon(reach, tags={"shop": "supermarket"}).to_crs(CRS)
stores = stores[shapely.distance(stores.geometry.representative_point().values, county) <= 10_000]
stores = stores.sort_index().reset_index()
footprint = stores.geometry.area.where(stores.geom_type.isin(["Polygon", "MultiPolygon"]))
stores["geometry"] = stores.geometry.representative_point()
# one store mapped twice (two outlines, or a point and an outline) is one store: keep the outline
keep = []
for i in footprint.notna().sort_values(ascending=False, kind="stable").index:
if not any(stores.name[i] == stores.name[j] and stores.geometry[i].distance(stores.geometry[j]) < 50 for j in keep):
keep.append(i)
keep = sorted(keep)
print(f"{len(stores) - len(keep)} duplicate features dropped")
stores, footprint = stores.loc[keep].reset_index(drop=True), footprint.loc[keep].reset_index(drop=True)
A = footprint.fillna(footprint.median()).to_numpy()
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()
s_idx = [index[n] for n in ox.distance.nearest_nodes(G, stores.geometry.x, stores.geometry.y)]
h_idx = 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=s_idx)[:, h_idx].T / 60 # home x store, minutes
ok = np.isfinite(T_all).all(axis=1)
T, pop = T_all[ok], homes.loc[ok, "POP20"].to_numpy(float)
T_minutes, footprints = T, footprint
print(f"{len(stores)} supermarkets; {ok.sum():,} blocks, {pop.sum():,.0f} people "
f"(dropped {(~ok).sum()} blocks, {homes.POP20[~ok].sum():,} people)")2 duplicate features dropped
39 supermarkets; 2,239 blocks, 167,770 people (dropped 2 blocks, 553 people)
Example 1 — a diagnostic for the largest-store effect
import numpy as np
def largest_store_diagnostic(A, T, population, alpha=1.0, beta=2.0, floor=0.5):
"""Compare the model's pull of the largest store with the nearest-store baseline."""
A = np.asarray(A, float)
Tf = np.maximum(np.asarray(T, float), floor)
U = A[None, :] ** alpha * Tf ** -beta
P = U / U.sum(axis=1, keepdims=True)
largest = int(np.argmax(A))
w = population / population.sum()
report = {
"floor_area_share": A[largest] / A.sum(),
"market_share": float((P[:, largest] * w).sum()),
"likeliest_is_largest": float(w[P.argmax(axis=1) == largest].sum()),
"nearest_is_largest": float(w[Tf.argmin(axis=1) == largest].sum()),
}
print(" | ".join(f"{k} {v:.1%}" for k, v in report.items()))
return report
largest_store_diagnostic(A, T, pop)floor_area_share 6.7% | market_share 8.5% | likeliest_is_largest 12.7% | nearest_is_largest 9.3%
These four numbers are the sanity check. A likeliest-is-largest share several times the nearest-is-largest share, or a market share several times the floor-area share, points to one of the mistakes above.
Example 2 — sweeping the decay parameter
def decay_sweep(A, T, population, betas=(0.0, 0.5, 1.0, 2.0, 5.0)):
for beta in betas:
r = largest_store_diagnostic(A, T, population, beta=beta)
decay_sweep(A, T, pop)floor_area_share 6.7% | market_share 6.7% | likeliest_is_largest 100.0% | nearest_is_largest 9.3%
floor_area_share 6.7% | market_share 7.8% | likeliest_is_largest 27.5% | nearest_is_largest 9.3%
floor_area_share 6.7% | market_share 8.4% | likeliest_is_largest 14.5% | nearest_is_largest 9.3%
floor_area_share 6.7% | market_share 8.5% | likeliest_is_largest 12.7% | nearest_is_largest 9.3%
floor_area_share 6.7% | market_share 9.3% | likeliest_is_largest 11.3% | nearest_is_largest 9.3%
The rows are β = 0, 0.5, 1, 2 and 5. A curve that flattens between 1 and 2, as here, says that the exact value in that range matters little for the household map; a curve that is still falling steeply says it matters a lot.
Example 3 — refusing unit-dependent decay without a unit
def distance_decay(T, form="power", beta=None, lam=None, unit=None):
"""Distance decay that refuses to guess units for the forms that depend on them."""
T = np.asarray(T, float)
if form == "power":
return T ** -beta
if form == "exponential":
if unit not in ("s", "min", "h"):
raise ValueError("exponential decay needs the unit lambda was fitted in: 's', 'min' or 'h'")
return np.exp(-lam * T)
raise ValueError(f"unknown form {form!r}")Making the unit a required argument does not convert anything; it makes whoever calls the function state what λ means, which is where the 59.1% error came from.
Example 4 — the unit, exponent and fill checks
Tf = np.maximum(T, 0.5)
largest = int(np.argmax(A))
def likeliest_is_largest(U):
P = U / U.sum(axis=1, keepdims=True)
return pop[P.argmax(axis=1) == largest].sum() / pop.sum()
for unit, factor in (("minutes", 1), ("seconds", 60), ("hours", 1 / 60)):
print(f"exponential lambda 0.2 per minute, times in {unit:7}: {likeliest_is_largest(A * np.exp(-0.2 * Tf * factor)):.1%}")
for unit, factor in (("minutes", 1), ("hours", 1 / 60)):
print(f"power decay on (T + 1), T in {unit:7}: {likeliest_is_largest(A * (Tf * factor + 1) ** -2.0):.1%}")
for alpha in (0.5, 1.0, 2.0):
largest_store_diagnostic(A, T, pop, alpha=alpha)
for label, fill in (("median", footprints.median()), ("maximum", footprints.max()), ("1 m²", 1.0)):
a = footprints.fillna(fill).to_numpy()
U = a[None, :] * Tf ** -2.0
share = (U / U.sum(axis=1, keepdims=True) * pop[:, None]).sum(axis=0) / pop.sum()
print(f"missing footprints filled with the {label:7}: largest store's share {share[largest]:.1%}")
print(f"stores without a footprint: {footprints.isna().sum()} of {len(footprints)}; zero travel times: {(T_minutes == 0).sum()}")exponential lambda 0.2 per minute, times in minutes: 18.5%
exponential lambda 0.2 per minute, times in seconds: 9.3%
exponential lambda 0.2 per minute, times in hours : 59.1%
power decay on (T + 1), T in minutes: 12.9%
power decay on (T + 1), T in hours : 39.0%
floor_area_share 6.7% | market_share 6.2% | likeliest_is_largest 11.6% | nearest_is_largest 9.3%
floor_area_share 6.7% | market_share 8.5% | likeliest_is_largest 12.7% | nearest_is_largest 9.3%
floor_area_share 6.7% | market_share 13.2% | likeliest_is_largest 14.5% | nearest_is_largest 9.3%
missing footprints filled with the median : largest store's share 8.5%
missing footprints filled with the maximum: largest store's share 5.9%
missing footprints filled with the 1 m² : largest store's share 11.2%
stores without a footprint: 15 of 39; zero travel times: 7
Weak decay favours size
Why weak decay favours size
In the Huff model a store's pull is its attractiveness divided by a power of distance. When the power is small, distance barely discounts a store, so the store with the largest numerator wins every household's comparison. As β grows, nearby stores' denominators become much smaller than distant stores', and local stores win their own neighbourhoods.
Why exponential decay is fragile
exp(−λT) compares
Why missing values move shares
Attractiveness enters every household's normalisation. Giving 15 stores a near-zero attractiveness removes them from the competition for everyone, and their customers spread to the remaining stores in proportion to those stores' pull. The largest store takes the largest part of that windfall.
Why the share and the map can disagree
A store's share sums probabilities; the likeliest-store map takes each household's maximum. A store can be every household's second choice and still earn a large share, or be the likeliest for many households with modest probabilities each. Checking only one of the two lets the other go wrong unseen.
Zero attractiveness and degree distances
- Zero attractiveness removes a store. A store with
A = 0has zero probability for every household. - Distances in degrees behave like a unit error. Project coordinates or use network times.
- A single dominant store can be real. A regional hypermarket may win widely; compare with observed data before "fixing" it.
- Stores outside the study area must be included or in-area stores absorb their customers.
- Constrained models behave differently. If store capacity is enforced, the largest store cannot absorb unlimited demand.
- Aggregated origins hide the effect. Tract-level origins average away the household map; use small units.
- Report α, β, the decay form, the time unit and the fill rule with any published result.
Internal links
- Gravity and Huff models explained: predicting where people go — what the parameters mean
- How to estimate market share with a Huff model in Python — the complete workflow
- Accessibility measures explained: nearest, cumulative and gravity — the same decay in access scores
- How to measure distance to the nearest facility for every home — the nearest-store baseline
- Fixing accessibility scores that are wrong near the study area edge — missing stores beyond the boundary
- Catchment areas explained: buffers, isochrones and Voronoi compared — deterministic catchments for comparison
- How to handle missing and null values in spatial datasets — filling attractiveness honestly
- Location analysis explained: catchments, accessibility and site selection — where market share fits
FAQ
Why does my Huff model assign everyone to the largest store?
Distance decay is too weak or not working. With β = 0 the largest store was every household's likeliest store; at β = 0.5 it was for 27.5%; at β = 2 for 12.7%, close to the 9.3% for whom it is actually nearest.
Does the unit of travel time matter?
For power decay, no: minutes and seconds gave identical probabilities. For exponential decay, yes: a per-minute λ applied to hours made the largest store the likeliest for 59.1% of residents.
What attractiveness exponent should I use?
α = 1 unless data support something else. At α = 2 the largest store's share rose from 8.5% to 13.2%, and the top three stores took a third of the market.
How should I handle stores with no size data?
Fill them with a stated, defensible value, or find a real measure. Filling 15 missing footprints with 1 m² raised the largest store's share from 8.5% to 11.2%.
Why do some rows come out as NaN?
A travel time of zero raised to a negative power is infinite. Seven origin–store pairs had zero time; floor times at a small positive value.
How do I know the model is reasonable?
Compare the share of people whose likeliest store is the largest with the share for whom it is the nearest, and the largest store's market share with its share of floor area. Better still, calibrate against observed shopping data.