Gravity and Huff Models Explained: Predicting Where People Go

On this page

Problem statement

A catchment says who is near a store. It does not say where people actually shop, because people do not always use the nearest store: a larger supermarket a few minutes further away draws custom past a smaller one next door. Gravity models describe that trade-off, and the Huff model is the version most used for retail: each household spreads its probability of using each store according to the store's attractiveness and the time it takes to get there.

The model has two parameters that people set by habit — a distance-decay exponent and an attractiveness exponent — and they behave very differently. Measured with 39 supermarkets in and around Chittenden County, Vermont, attractiveness taken as building footprint, and drive times from 2,239 populated census blocks:

The distance exponent barely changed market shares. The largest store's share of the county's shoppers stayed between 6.7% and 9.3% for every exponent from 0 to 5.
It completely changed the map of which store each household most likely uses. The largest store was the likeliest store for 100% of residents with no distance decay and for 12.7% with an exponent of 2.
The attractiveness exponent moved shares directly. Squaring floor area raised the largest store's share from 8.5% to 13.2%.

Quick answer

The Huff probability that households in area i use store j:

P(i, j) = A_j^α × T_ij^−β  /  Σ_k A_k^α × T_ik^−β

A is attractiveness (floor area, turnover, number of checkouts), T is travel time or distance, α weights attractiveness and β sets how fast the pull of a store fades with distance. A store's expected share of customers is the population-weighted sum of its probabilities.

import numpy as np

U = A[None, :] ** alpha * np.maximum(T, 0.5) ** -beta     # utility, homes x stores (minutes)
P = U / U.sum(axis=1, keepdims=True)                        # each row sums to 1
share = (P * population[:, None]).sum(axis=0) / population.sum()

Calibrate β and α against observed customer data where you have it. Where you do not, report the parameters, and show how the results you care about change across a plausible range.

Every home spreads one unit of probability across all stores; attractiveness pulls it one way and distance the other.

Attractiveness, travel time with a floor, then beta and alpha

1. Choose an attractiveness measure

Attractiveness stands in for everything that makes a store worth the trip. Floor area is the usual open proxy; here it came from OpenStreetMap building footprints, available for 24 of 39 supermarkets once two stores mapped twice were merged, with a median of 1,870 m² and a maximum of 6,499 m². Stores mapped only as points need a value, and how you fill the gap is itself a modelling choice (step 6).

The units do not matter. Footprints in square feet gave probabilities identical to those in square metres, to within floating-point rounding (under 10⁻¹⁵), because a constant factor cancels in the ratio.

2. Build a travel-time matrix

Rows are origins — homes or census blocks — and columns are stores. Network drive times are better than straight lines, for the same reasons as in any accessibility analysis. Include stores beyond the study area, since residents near the boundary shop across it.

3. Floor the travel times

Seven origin–store pairs had a travel time of exactly zero, because a block and a store snapped to the same network node. Raised to a negative power, zero becomes infinity, and those seven rows came out as NaN. A floor — half a minute here — keeps the model defined without changing anything else measurably.

4. Understand what β controls

With β = 0, distance plays no part: every home splits its custom in proportion to floor area, so every store's share equals its share of floor area — 6.7% for the largest — and the largest store is every household's single likeliest store.

As β grows, nearby stores take over each home's probabilities:

β      largest store's share   top 3 stores' share   homes whose likeliest store is the largest
0.0           6.7%                   20.2%                        100.0%
0.5           7.8%                   22.3%                         27.5%
1.0           8.4%                   22.7%                         14.5%
2.0           8.5%                   21.7%                         12.7%
5.0           9.3%                   23.2%                         11.3%

The shares hardly move because stores are spread around the county: a large exponent concentrates each home on its local store, and homes split among local stores roughly as they did among large ones. The household map is what β sets.

5. Understand what α controls

α changes how much size matters relative to distance, and the effect on shares is direct:

α      largest store's share   top 3 stores' share
0.5           6.2%                   18.7%
1.0           8.5%                   21.7%
2.0          13.2%                   32.5%

With β = 2, doubling α concentrated a third of all shopping in three stores (Example 4). If the question is how much a new large store would take from existing ones, α matters more than β.

6. Decide how to fill missing attractiveness

Filling the 15 footprint-less stores with the median footprint gave the largest store 8.5%. Filling them with the maximum footprint gave it 5.9%; filling them with 1 m² gave it 11.2% and pushed the top three to 29.1%. Missing data in A is not a detail. Prefer a real measure for every store, or at least report the fill rule.

7. Use the right form of distance decay for your units

The power form, $T^{-\beta}$ , is scale-free: travel times in minutes or seconds gave identical probabilities. The exponential form, exp(−λT), is not. λ = 0.2 per minute applied to times in hours made the largest store the likeliest store for 59.1% of residents instead of 18.5%. Any additive term, such as $(T + 1)^{-\beta}$ , also depends on units.

Market shares barely change with β; which store each household most likely uses changes completely.

Huff probabilities and calibration

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). Two stores were in OpenStreetMap twice — the same Hannaford as two identical outlines, and a market as both a point and an outline — and are merged. 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-huff-models-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]
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)
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 — Huff probabilities with the safeguards

import numpy as np


def huff_probabilities(A, T, alpha=1.0, beta=2.0, floor=0.5):
    """Probability of each origin using each store; rows sum to 1. T in any unit, floored."""
    A = np.asarray(A, dtype=float)
    if np.isnan(A).any():
        raise ValueError(f"{int(np.isnan(A).sum())} stores have no attractiveness value")
    T = np.maximum(np.asarray(T, dtype=float), floor)
    utility = A[None, :] ** alpha * T ** -beta
    return utility / utility.sum(axis=1, keepdims=True)

Refusing missing attractiveness forces the fill rule into the open, and the floor removes the zero-time NaN rows. floor must be in the same unit as T.

A duplicate is invisible on a map and very visible in a share.

Example 2 — shares and the likeliest store together

def huff_summary(P, population, attractiveness, label):
    """Expected shares, plus how often the largest store is each area's likeliest."""
    share = (P * population[:, None]).sum(axis=0) / population.sum()
    largest = int(np.argmax(attractiveness))
    top3 = np.sort(share)[-3:].sum()
    likeliest_is_largest = population[P.argmax(axis=1) == largest].sum() / population.sum()
    print(f"{label:12} largest {share[largest]:6.1%} | top-3 {top3:6.1%} | "
          f"likeliest is largest {likeliest_is_largest:6.1%}")
    return share

for beta in (0.0, 0.5, 1.0, 2.0, 5.0):
    huff_summary(huff_probabilities(A, T, beta=beta), pop, A, f"beta {beta}")

beta 0.0     largest   6.7% | top-3  20.2% | likeliest is largest 100.0%
beta 0.5     largest   7.8% | top-3  22.3% | likeliest is largest  27.5%
beta 1.0     largest   8.4% | top-3  22.7% | likeliest is largest  14.5%
beta 2.0     largest   8.5% | top-3  21.7% | likeliest is largest  12.7%
beta 5.0     largest   9.3% | top-3  23.3% | likeliest is largest  11.3%

Report both columns. A model can have sensible shares and an absurd household map, or the reverse.

Example 3 — calibrating β against observed shares

from scipy.optimize import minimize_scalar


def calibrate_beta(A, T, population, observed_share, alpha=1.0, bounds=(0.0, 5.0)):
    """The beta whose predicted store shares are closest to observed shares (sum of squares)."""
    def loss(beta):
        P = huff_probabilities(A, T, alpha=alpha, beta=beta)
        predicted = (P * population[:, None]).sum(axis=0) / population.sum()
        return float(((predicted - observed_share) ** 2).sum())
    result = minimize_scalar(loss, bounds=bounds, method="bounded")
    return result.x, result.fun


P2 = huff_probabilities(A, T, beta=2.0)
observed = (P2 * pop[:, None]).sum(axis=0) / pop.sum()          # a test: shares the model itself made
fitted, _ = calibrate_beta(A, T, pop, observed)
for b in (1.5, 2.5):
    P = huff_probabilities(A, T, beta=b)
    loss = float((((P * pop[:, None]).sum(axis=0) / pop.sum() - observed) ** 2).sum())
    print(f"loss at beta {b}: {loss:.1e}")
print(f"fitted beta: {fitted:.3f}")

loss at beta 1.5: 4.7e-04
loss at beta 2.5: 3.6e-04
fitted beta: 2.000

observed_share comes from loyalty-card data, a survey or turnover figures, in the order of A. As a test, shares generated by the model with β = 2 were fitted back to exactly 2.000 — but the loss was only 4.7 × 10⁻⁴ at β = 1.5 and 3.6 × 10⁻⁴ at β = 2.5. That shallow minimum is what share insensitivity looks like: noise in real share data can move the fitted β a long way, and data on where individual households shop pins it down far better.

Example 4 — the other measurements on this page

print(f"footprints for {footprint.notna().sum()} of {len(stores)} stores: median {footprint.median():,.0f} m², "
      f"max {footprint.max():,.0f} m²")
with np.errstate(divide="ignore", invalid="ignore"):
    unfloored = A[None, :] * T ** -2.0                  # zero time -> infinite utility
    unfloored = unfloored / unfloored.sum(axis=1, keepdims=True)
print(f"origin-store pairs with zero travel time: {(T == 0).sum()}; "
      f"rows that come out NaN without a floor: {np.isnan(unfloored).any(axis=1).sum()}")

for alpha in (0.5, 1.0, 2.0):
    huff_summary(huff_probabilities(A, T, alpha=alpha), pop, A, f"alpha {alpha}")

filled = {"median": footprint.fillna(footprint.median()), "maximum": footprint.fillna(footprint.max()),
          "1 m²": footprint.fillna(1.0)}
for label, a in filled.items():
    huff_summary(huff_probabilities(a.to_numpy(), T), pop, A, f"fill {label}")

feet = huff_probabilities(A * 10.7639, T)
secs = huff_probabilities(A, T * 60, floor=30)
print(f"square feet vs metres: max difference {np.abs(feet - huff_probabilities(A, T)).max():.1e}; "
      f"seconds vs minutes: {np.abs(secs - huff_probabilities(A, T)).max():.1e}")
largest = int(np.argmax(A))
for unit, t in (("minutes", np.maximum(T, 0.5)), ("hours", np.maximum(T, 0.5) / 60)):
    U = A[None, :] * np.exp(-0.2 * t)
    P = U / U.sum(axis=1, keepdims=True)
    print(f"exponential decay, lambda 0.2 applied to {unit}: likeliest store is the largest for "
          f"{pop[P.argmax(axis=1) == largest].sum() / pop.sum():.1%} of residents")

footprints for 24 of 39 stores: median 1,870 m², max 6,499 m²
origin-store pairs with zero travel time: 7; rows that come out NaN without a floor: 7
alpha 0.5    largest   6.2% | top-3  18.7% | likeliest is largest  11.6%
alpha 1.0    largest   8.5% | top-3  21.7% | likeliest is largest  12.7%
alpha 2.0    largest  13.2% | top-3  32.5% | likeliest is largest  14.5%
fill median  largest   8.5% | top-3  21.7% | likeliest is largest  12.7%
fill maximum largest   5.9% | top-3  21.6% | likeliest is largest  10.5%
fill 1 m²    largest  11.2% | top-3  29.1% | likeliest is largest  12.9%
square feet vs metres: max difference 4.4e-16; seconds vs minutes: 5.6e-16
exponential decay, lambda 0.2 applied to minutes: likeliest store is the largest for 18.5% of residents
exponential decay, lambda 0.2 applied to hours: likeliest store is the largest for 59.1% of residents

Why gravity, and what each exponent does

Why gravity

Newton's law of gravitation made the analogy: attraction grows with mass and falls with distance. Reilly applied it to retail trade between towns in the 1930s, and Huff turned it into probabilities in the 1960s, so that one household can use several stores with different likelihoods. The form survives because it is simple, gives every store a non-zero share, and has two parameters with a clear meaning.

Why shares and household maps respond differently to β

A store's share sums probabilities over every household. With weak decay, each household spreads its custom by store size, so large stores get a bit of everyone. With strong decay, each household concentrates on nearby stores, so every store gets most of its local households. Across a county with stores in every town, those two pictures add up to similar totals. For each household, though, the likeliest store switches from the largest to the nearest. Any analysis at the household or neighbourhood level depends heavily on β.

Why α is the lever for size

α raises attractiveness to a power before comparing stores. At α = 2, a store with twice the floor area has four times the pull, and large stores capture custom well beyond their neighbourhoods. That is why impact studies for a proposed large store are more sensitive to α than to β, and why α is the parameter most worth calibrating.

Why power decay is safer than exponential decay

A power function of time is scale-invariant: multiplying every time by 60 multiplies every utility by the same factor, which cancels. An exponential of time is not; λ has units of inverse time, and applying a per-minute λ to seconds or hours makes decay 60 times stronger or weaker. Exponential decay is perfectly good when λ is stated with its unit and fitted to data; it is dangerous when copied from a paper that used a different unit.

Where β reshapes the household map, α reshapes the market itself.

Zero travel times, missing values and units

A store with zero travel time breaks the power form. Floor the times.
Missing attractiveness must be filled deliberately. Median, maximum and 1 m² fills gave the largest store 8.5%, 5.9% and 11.2%.
Attractiveness units cancel; decay units only cancel for the power form.
Stores beyond the boundary take custom. Leaving them out inflates every in-area store's share.
Probabilities are not visits. Multiply by spending or trip rates to get a volume.
Uncalibrated parameters are assumptions. β = 2 is a convention, not a fact about any town.
The model has no capacity. A store predicted to take 30% of a county's shoppers does not get busier; add a constraint if congestion matters.

Internal links

How to estimate market share with a Huff model in Python — the full workflow
Fixing a gravity model that sends everyone to the biggest store — when the parameters go wrong
Accessibility measures explained: nearest, cumulative and gravity — distance decay used for access
Catchment areas explained: buffers, isochrones and Voronoi compared — the deterministic alternative
Location analysis explained: catchments, accessibility and site selection — where gravity models fit
How to measure network distance for many origin–destination pairs — the travel-time matrix
Location-allocation explained: p-median, coverage and what each optimises — choosing sites once demand is known
Spatial interpolation explained — another use of inverse-distance weighting

FAQ

What is the Huff model?

A probabilistic retail model in which each household's chance of using a store rises with the store's attractiveness and falls with travel time, normalised so that the probabilities across all stores sum to one.

What does the distance-decay exponent β do?

It sets how quickly a store's pull fades with distance. For Chittenden County supermarkets it changed each household's likeliest store dramatically — the largest store was likeliest for 100% of residents at β = 0 and 12.7% at β = 2 — while store market shares stayed within 6.7–9.3%.

What does the attractiveness exponent α do?

It sets how much size matters. Raising α from 1 to 2 increased the largest store's share from 8.5% to 13.2%.

What value of β should I use?

Calibrate it against observed shopping data. Without data, β around 1–2 for drive times is common in practice, but report the value and test how sensitive your conclusions are to it.

Does it matter whether travel time is in minutes or seconds?

Not for the power form, where the units cancel. For exponential decay it matters a great deal: a per-minute λ applied to hours made the largest store the likeliest for 59.1% of residents.

What should I use for store attractiveness?

Whatever best predicts patronage and exists for every store: floor area, turnover, checkouts or product range. Fill missing values deliberately, because the fill rule changed the largest store's share by nearly half.