Geostatistical Modeling Spatial Statistics with Python

Spatially Constrained Clustering with Max-P Regions

TL;DR: Standardise your attributes, build a contiguity weights object with libpysal.weights.Queen.from_dataframe(gdf), then fit spopt.region.MaxPHeuristic(gdf, w, attrs_name=attrs, threshold_name="pop", threshold=50000) and call .solve(). Read model.labels_ for the region of each unit and model.p for the number of regions the heuristic discovered under your size floor.

Why This Matters

Max-P-Regions solves a problem that ordinary clustering cannot: build the maximum number of contiguous regions where each region meets a minimum-size requirement. That is exactly the constraint behind statistical reporting zones, minimum-population districts, and service territories with a viability floor. This walkthrough is the applied companion to the spatial clustering and regionalization guide, which develops the objective functions and contrasts Max-P with SKATER and Ward. Both sit within the broader Core Concepts of Spatial Statistics & Geostatistics framework. Here the focus is running spopt end to end.

Environment

bash
pip install "geopandas>=1.0" "libpysal>=4.9" "spopt>=0.6" "scikit-learn>=1.2" "numpy>=1.23" "matplotlib>=3.7"
python
import geopandas as gpd
import numpy as np
import libpysal
from spopt.region import MaxPHeuristic
from sklearn.preprocessing import StandardScaler
from shapely.geometry import box
import matplotlib.pyplot as plt

Step-by-Step Implementation

Step 1 — Build or load areal units in a projected CRS

python
# Synthetic 10x10 grid of areal units — replace with your own file:
# gdf = gpd.read_file("counties.gpkg").to_crs("EPSG:5070")
cells = [box(x, y, x + 1, y + 1) for y in range(10) for x in range(10)]
gdf = gpd.GeoDataFrame(geometry=cells, crs="EPSG:5070")
gdf = gdf.reset_index(drop=True)

Use a projected CRS and a reset integer index so weight rows align with attribute rows.

Step 2 — Create homogeneity attributes and a threshold variable

python
rng = np.random.default_rng(7)
cx = gdf.geometry.centroid.x.values
gdf["income"] = 40000 + cx * 2500 + rng.normal(0, 3000, len(gdf))  # spatial gradient
gdf["pct_college"] = 0.2 + cx * 0.04 + rng.normal(0, 0.03, len(gdf))
gdf["pop"] = rng.integers(8000, 22000, len(gdf)).astype(float)       # threshold var

The homogeneity attributes carry a west-to-east gradient so coherent regions exist. pop is the spatially extensive variable the size floor applies to.

Step 3 — Standardise the clustering attributes

python
attrs = ["income", "pct_college"]
gdf[attrs] = StandardScaler().fit_transform(gdf[attrs].values)

Standardisation is mandatory: without it, income (tens of thousands) would swamp pct_college (fractions) and regions would form along income alone.

Step 4 — Build a contiguity weights object

python
w = libpysal.weights.Queen.from_dataframe(gdf)
assert w.n_components == 1, "Max-P needs a single connected component"

Max-P merges only adjacent units, so the weight matrix must be fully connected. Resolve any islands before proceeding — see spatial weight matrices.

Step 5 — Fit MaxPHeuristic

python
model = MaxPHeuristic(
    gdf,
    w,
    attrs_name=attrs,       # variables minimised for within-region homogeneity
    threshold_name="pop",   # floor applies to this column
    threshold=60000,        # each region must total >= 60,000 population
    top_n=2,                # candidate solutions kept during local search
    max_iterations_construction=99,
)
model.solve()

MaxPHeuristic discovers the region count rather than taking it as input. The floor of 60,000 forces each region to aggregate roughly three to five cells here. The top_n argument controls how many candidate assignments the local-search phase keeps at each step: higher values explore more of the solution space at the cost of runtime, while max_iterations_construction bounds the greedy region-growing phase. Both are quality/speed levers, not the region count — that is fixed only by the data and the floor.

Step 6 — Extract labels, evaluate, and map

python
gdf["region"] = model.labels_
print(f"Max-P discovered {model.p} regions")

# Within-region variance of the standardised attributes (lower = tighter)
wrv = sum(grp[attrs].var(ddof=0).sum() * len(grp)
          for _, grp in gdf.groupby("region")) / len(gdf)
print(f"Mean within-region variance: {wrv:.3f}")

fig, ax = plt.subplots(figsize=(7, 7))
gdf.plot(column="region", categorical=True, cmap="tab20",
         edgecolor="0.7", linewidth=0.3, ax=ax)
ax.set_title(f"Max-P regions (p = {model.p}, floor = 60,000 pop)")
ax.set_axis_off()
fig.savefig("maxp_regions.png", dpi=150, bbox_inches="tight")

Interpreting the Output

  • model.labels_ — the integer region index for each unit, in GeoDataFrame row order. Attach it as a column and map it.
  • model.p — the number of regions the heuristic found. It is an output, not an input; a higher floor yields fewer, larger regions.
  • Within-region variance — the mean spread of the standardised attributes inside regions. Because each standardised variable has unit variance, the global reference is the number of attributes (here 2.0); a well below-reference value means regions are internally homogeneous.
  • The map — every region should be one contiguous block. A fragmented region signals a weights or island problem, not a valid Max-P output.

Because the synthetic attributes carry a west-to-east income gradient, expect regions to elongate roughly north-south, stacking cells of similar income into vertical bands while satisfying the population floor. Compare the within-region variance against the global reference of 2.0 (two standardised attributes): a value around 0.3 to 0.6 indicates the partition captured most of the attribute structure, whereas a value approaching 2.0 means contiguity forced together units with little in common — usually a sign the attributes lack spatial structure, or that the floor is so high it overrides homogeneity. You can also inspect the population total per region to confirm the floor held:

python
region_pop = gdf.groupby("region")["pop"].sum()
print(region_pop)
assert (region_pop >= 60000).all(), "Floor violated — check the solver run"

Critical Best Practices

Set the floor from a real requirement

The threshold is a policy input, not a tuning knob. Derive it from the actual minimum-size rule (minimum population, minimum count of events) rather than trial and error. Too high a floor collapses everything into one or two giant regions; too low a floor makes the constraint inert and the result resembles unconstrained clustering.

Always standardise, and pick attributes deliberately

Regionalization minimises variance across attrs_name, so unstandardised or redundant variables distort the partition. Standardise every attribute and include only variables you genuinely want the regions to be homogeneous in. Two attributes that are strongly correlated effectively double-weight the same dimension, so drop or combine near-duplicates before fitting. If you have many candidate variables, reduce them first (for example with a principal component step) so the objective is not dominated by an accidental cluster of collinear features.

Confirm a single connected component

MaxPHeuristic cannot bridge disconnected graphs. Assert w.n_components == 1 before solving; if it fails, switch to KNN weights or bridge islands manually as described in spatial weight matrices.

Record the run for reproducibility

The heuristic uses randomised construction and local search, so results can vary between runs. Pin spopt and libpysal versions, and where the API exposes a seed, fix it so a reported partition can be regenerated.

Troubleshooting

Symptom Likely cause Fix
Connectivity / graph error on solve Weight matrix has multiple components assert w.n_components == 1; resolve islands first
model.p is 1 or 2 Floor too high for total threshold sum Lower threshold; verify units of threshold_name
Regions split along one variable Attributes not standardised Apply StandardScaler before fitting
Partition changes every run Randomised heuristic, no seed pinned Pin versions; fix the seed where available
Very long runtime Large nn with tight floor and high top_n Reduce top_n; pre-aggregate units
ValueError on labels_ length Index not reset before weights gdf.reset_index(drop=True) before building w

Validate against an unconstrained baseline

Fit a plain k-means on the same standardised attributes and compare its within-cluster variance to Max-P’s within-region variance. The gap between them is the price the contiguity constraint imposes: a small gap means geography and attributes agree, while a large gap warns that forcing contiguity has badly compromised homogeneity, which may mean a different floor or a different variable set is needed.

Next Steps

Return to the spatial clustering and regionalization guide to compare Max-P with SKATER and Ward, and revisit spatial weight matrices to make sure the contiguity graph underpinning the constraint is sound.

← Back to Spatial Clustering & Regionalization

Related