SKATER vs Max-P Regionalization in Python

TL;DR: Use spopt.region.Skater(gdf, w, attrs_name, n_clusters=k) when the number of regions is fixed by the decision, and spopt.region.MaxPHeuristic(gdf, w, attrs_name, threshold_name="population", threshold=40_000) when a minimum region size is fixed by the data or by disclosure rules. SKATER gives you k and no size guarantee; Max-P gives you a size guarantee and picks k itself. Standardise the attributes first, or income wins.

Why This Matters

Both algorithms take the same three inputs — areal units, attributes, a contiguity graph — and both return contiguous regions. That surface similarity hides the fact that they answer different questions, and picking the wrong one produces a partition nobody can use. SKATER asks: given that I must end up with exactly nineteen regions, which nineteen are most internally homogeneous? Max-P asks: given that every region must hold at least forty thousand people, how many regions can I get away with? One of those is a constraint you inherit from a decision — nineteen delivery depots, nineteen sampling strata — and the other is a constraint you inherit from the data, from statistical disclosure control, or from a minimum sample size below which an estimate is not publishable.

Regionalization is usually a response to the modifiable areal unit problem: rather than inherit arbitrary administrative zones, you derive units from the data so the zoning is at least reproducible and stated. That only helps if the derived units are fit for the analysis that follows, which is exactly what the size distribution decides. A region holding 7,412 people cannot support a rate estimate that a region holding 135,819 people supports comfortably, and a partition that mixes the two produces a map whose reliability varies by a factor of eighteen across the study area. The wider treatment of these methods sits in Spatial Clustering & Regionalization, itself part of the Core Concepts of Spatial Statistics & Geostatistics.

Environment and Version Pinning

Both algorithms live in spopt, which sits on top of libpysal for the graph and scikit-learn for the dissimilarity metric. Pin spopt explicitly: the region module changed argument names between minor releases and older tutorials will not run against 0.6.

bash

pip install "spopt>=0.6.0" "libpysal>=4.9.0" "geopandas>=1.0" \
            "scikit-learn>=1.3.0" "numpy>=1.24" "pandas>=2.0" "shapely>=2.0"

python

import numpy
import pandas as pd
import geopandas as gpd
import libpysal
from shapely.geometry import box
from sklearn.metrics import pairwise as skm
from sklearn.preprocessing import StandardScaler
from sklearn.metrics import adjusted_rand_score
from spopt.region import Skater, MaxPHeuristic

Step-by-Step Implementation

1. Build the areal dataset

A 20 by 20 lattice of 500 m cells gives 400 units with a known attribute structure: income rises to the east, degree attainment rises to the north, and owner occupation rises along the diagonal. The three trends deliberately disagree, so the partition is not a trivial banding.

python

rng = numpy.random.default_rng(2026)
SIDE, CELL = 20, 500.0

cells, cx, cy = [], [], []
for r in range(SIDE):
    for c in range(SIDE):
        cells.append(box(c * CELL, r * CELL, (c + 1) * CELL, (r + 1) * CELL))
        cx.append((c + 0.5) * CELL)
        cy.append((r + 0.5) * CELL)

cx, cy = numpy.asarray(cx), numpy.asarray(cy)
u = cx / (SIDE * CELL)          # 0..1 west to east
v = cy / (SIDE * CELL)          # 0..1 south to north
n = SIDE ** 2

gdf = gpd.GeoDataFrame(
    {
        "median_income": 31_000 + 9_000 * (1.4 * u - 0.6 * v) + rng.normal(0, 2_600, n),
        "pct_degree":    28.0 + 12.2 * (1.1 * v - 0.4 * u) + rng.normal(0, 3.4, n),
        "pct_owner_occ": 52.0 + 15.4 * (0.9 * u + 0.7 * v - 0.8) + rng.normal(0, 4.1, n),
        "population":    rng.integers(600, 4_400, n),
    },
    geometry=cells,
    crs="EPSG:27700",
).reset_index(drop=True)

print(f"{len(gdf)} areas, total population {gdf['population'].sum():,}")

text

400 areas, total population 996,110

That total matters for the Max-P run: a floor of 40,000 people admits at most 996110 // 40000 = 24 regions, so any p the heuristic returns has to sit at or below 24.

2. Build one contiguity graph and reuse it

The weights object is the constraint that makes this regionalization rather than plain clustering, and it must be identical across both models or the comparison is meaningless. Check the component count before going any further — both algorithms behave badly on a disconnected graph, and they behave badly in different ways.

python

rook = libpysal.weights.Rook.from_dataframe(gdf, use_index=True)
w = libpysal.weights.Queen.from_dataframe(gdf, use_index=True)

for name, obj in (("Rook", rook), ("Queen", w)):
    print(f"{name:<6} {obj.n} areas, {int(obj.s0 // 2):>5} joins, "
          f"{obj.n_components} component(s), {len(obj.islands)} island(s)")

text

Rook   400 areas,   760 joins, 1 component(s), 0 island(s)
Queen  400 areas,  1482 joins, 1 component(s), 0 island(s)

Note that s0 counts each join twice while the weights are still binary, hence the halving. Do not row-standardise here: unlike Moran’s I, the region models read only the neighbour lists, and a w.transform = "R" is at best wasted work. The choice between the two topologies is a real modelling decision, discussed further under spatial weight matrices; Queen nearly doubles the edge count and so gives the spanning tree far more cheap edges to work with.

3. Standardise the attributes

This is the step people skip, and skipping it silently reduces a three-variable regionalization to a one-variable one.

python

attrs = ["median_income", "pct_degree", "pct_owner_occ"]
raw_var = gdf[attrs].var(ddof=0)
print("variance share before scaling")
for a in attrs:
    print(f"  {a:<15} {raw_var[a] / raw_var.sum():.7f}")

attrs_name = ["z_income", "z_degree", "z_owner"]
gdf[attrs_name] = StandardScaler().fit_transform(gdf[attrs])

text

variance share before scaling
  median_income   0.9999968
  pct_degree      0.0000013
  pct_owner_occ   0.0000019

Income carries more than 99.999 per cent of the raw sum of squares purely because it is measured in pounds. After scaling, each of the three columns has unit variance, so the total sum of squares over 400 areas is exactly 3 x 400 = 1200 — a fixed denominator that makes the two partitions directly comparable later.

4. Fit SKATER with a fixed k

SKATER builds a minimum spanning tree over the contiguity graph, weighting each edge by the attribute dissimilarity of the two areas it joins, then removes the k - 1 edges whose removal most reduces within-region heterogeneity. Every cut splits one connected subtree into two, so contiguity is preserved by construction.

python

spanning_forest_kwds = dict(
    dissimilarity=skm.manhattan_distances,
    affinity=None,
    reduction=numpy.sum,
    center=numpy.mean,
)

skater = Skater(
    gdf,
    w,
    attrs_name,
    n_clusters=19,
    floor=3,                 # minimum AREAS per region, not minimum population
    trace=False,
    islands="increase",
    spanning_forest_kwds=spanning_forest_kwds,
)
skater.solve()
gdf["skater"] = skater.labels_
print(f"SKATER regions: {pd.Series(skater.labels_).nunique()}")

text

SKATER regions: 19

The floor argument is worth pausing on, because its name invites the wrong assumption. In Skater it is a quorum on the number of areas in a region, not a sum over an attribute column. Setting floor=3 guarantees no region is a singleton or a pair; it says nothing at all about population, and there is no way to make it say something about population.

5. Fit Max-P with a fixed floor

Max-P inverts the problem. Rather than fixing the region count and optimising homogeneity, it fixes a threshold on a summed attribute and maximises the number of regions that satisfy it, using homogeneity only as a secondary objective in the local search.

python

numpy.random.seed(123456)          # spopt's region heuristics read the global RNG

maxp = MaxPHeuristic(
    gdf,
    w,
    attrs_name,
    threshold_name="population",   # any numeric column in gdf
    threshold=40_000,              # every region must sum to at least this
    top_n=2,
)
maxp.solve()
gdf["maxp"] = maxp.labels_
print(f"Max-P found p = {maxp.p} regions at a 40,000 floor "
      f"(upper bound {gdf['population'].sum() // 40_000})")

text

Max-P found p = 19 regions at a 40,000 floor (upper bound 24)

The gap between 19 and the arithmetic bound of 24 is the price of contiguity. Population is not distributed so that every neighbourhood of areas summing to 40,000 is also connected, and the heuristic has to overshoot the floor in places to keep regions joined up. A p that lands well below the bound is normal; a p that lands at the bound usually means your floor is far too small to bind.

6. Compare the two partitions

With k set equal to the p that Max-P returned, the two partitions have the same number of regions over the same areas, so any difference is attributable to the objective rather than to granularity. Compare them on within-region sum of squares — the quantity SKATER explicitly minimises — and on the population distribution.

python

Z = gdf[attrs_name].to_numpy()
TSS = float(((Z - Z.mean(axis=0)) ** 2).sum())     # exactly 1200.0 after scaling

def summarise(col, floor=40_000):
    lab = gdf[col].to_numpy()
    pop = gdf.groupby(col)["population"].sum()
    cnt = gdf.groupby(col).size()
    twss = sum(((Z[lab == r] - Z[lab == r].mean(axis=0)) ** 2).sum()
               for r in numpy.unique(lab))
    return pd.Series({
        "p": len(numpy.unique(lab)),
        "units_min": cnt.min(), "units_med": int(cnt.median()), "units_max": cnt.max(),
        "pop_min": pop.min(), "pop_med": int(pop.median()), "pop_max": pop.max(),
        "below_floor": int((pop < floor).sum()),
        "twss": round(twss, 2),
    })

comparison = pd.DataFrame({"SKATER": summarise("skater"), "Max-P": summarise("maxp")})
print(comparison.to_string())
print(f"\ntotal SS = {TSS:.1f}")

text

             SKATER     Max-P
p                19        19
units_min         3        13
units_med        17        20
units_max        54        34
pop_min        7412     40318
pop_med       43180     48630
pop_max      135819     79142
below_floor       9         0
twss         486.31    531.74

total SS = 1200.0

7. Check run-to-run stability

SKATER is deterministic for a fixed graph and dissimilarity, up to ties in the edge costs. Max-P is not: it seeds regions at random and refines them with a randomised local search, so the same call under a different seed can return a different partition and occasionally a different p.

python

runs = []
for seed in range(10):
    numpy.random.seed(seed)
    m = MaxPHeuristic(gdf, w, attrs_name, "population", 40_000, top_n=2)
    m.solve()
    lab = numpy.asarray(m.labels_)
    twss = sum(((Z[lab == r] - Z[lab == r].mean(axis=0)) ** 2).sum()
               for r in numpy.unique(lab))
    runs.append({"seed": seed, "p": m.p, "twss": round(twss, 2)})

print(pd.DataFrame(runs).to_string(index=False))

text

 seed  p    twss
    0 19  531.74
    1 19  536.02
    2 18  548.61
    3 19  529.95
    4 19  533.18
    5 18  544.70
    6 19  530.41
    7 19  538.86
    8 18  551.23
    9 19  527.60

Interpreting the Output

Read the comparison table as two claims, one about homogeneity and one about usability. SKATER’s within-region sum of squares of 486.31 against a total of 1,200 means it accounts for 59.5 per cent of the attribute variation between regions; Max-P’s 531.74 accounts for 55.7 per cent. SKATER wins, and it should — homogeneity is its objective function, whereas for Max-P it is a tie-breaker applied after the floor is satisfied. Four percentage points is the cost of the guarantee.

The second claim is the one that usually decides the matter. SKATER’s regions run from 7,412 to 135,819 people, and nine of the nineteen hold fewer than 40,000. If the next step is publishing a rate per region, nine of those rates are built on samples too small to be stable, and the map will show apparent extremes that are entirely sampling noise. Max-P’s regions run from 40,318 to 79,142, a spread of 1.96 to 1 rather than 18.3 to 1, and every one of them clears the floor by construction. Note where the Max-P minimum sits: 40,318 is barely above the threshold, which is exactly what a working constraint looks like. If your Max-P minimum were 68,000 against a 40,000 floor, the floor is not binding and you should raise it.

The stability table is the third thing to read. Seven of ten seeds returned p = 19 and three returned 18, with within-region sum of squares between 527.60 and 551.23. That spread of 23.6 is about half the gap between the two algorithms, which is a useful calibration: quoting a single Max-P solution as though it were the answer overstates the precision by roughly the amount that separates it from SKATER in the first place. Good behaviour is a modal p that holds across most seeds with a narrow sum-of-squares band. A warning sign is p bouncing between three or more values, which means the floor is close to a structural limit of the graph and small random choices tip regions over or under it.

Critical Best Practices

Never run either algorithm on unstandardised attributes

The dissimilarity is a distance in raw attribute space, so a column measured in pounds annihilates one measured in percentage points. The demonstration is quick: build a partition from income alone, then compare it against the multivariate partitions with and without scaling.

python

income_only = Skater(gdf, w, ["median_income"], n_clusters=19,
                     spanning_forest_kwds=spanning_forest_kwds)
income_only.solve()

raw = Skater(gdf, w, attrs, n_clusters=19,
             spanning_forest_kwds=spanning_forest_kwds)
raw.solve()

print(f"ARI, raw attributes vs income only     : "
      f"{adjusted_rand_score(raw.labels_, income_only.labels_):.2f}")
print(f"ARI, z-scored attributes vs income only: "
      f"{adjusted_rand_score(gdf['skater'], income_only.labels_):.2f}")

text

ARI, raw attributes vs income only     : 0.94
ARI, z-scored attributes vs income only: 0.31

An adjusted Rand index of 0.94 means the unscaled three-variable run reproduced the income-only partition almost exactly. The other two columns contributed nothing.

Set the seed before `solve()`, and report a distribution

MaxPHeuristic reads NumPy’s global random state, so numpy.random.seed(...) immediately before the call is what makes a run reproducible — a default_rng generator held elsewhere in your script will not do it. Beyond reproducibility, run at least ten seeds and publish the modal p with the sum-of-squares band, as above. A single heuristic solution presented without its variability is a claim you cannot defend when a reviewer runs the notebook.

Match the weights to the geography, then keep them fixed

Rook and Queen contiguity gave 760 and 1,482 joins on the same lattice. That is not a detail: SKATER’s tree is built from those edges, and Max-P’s feasible merges are exactly those edges. On real polygon data, sliver geometries and coastlines create spurious or missing joins that change both answers, so inspect w.islands and w.n_components before solving and repair the graph rather than letting the models absorb the damage silently.

Do not compare within-region sums of squares across different k

The within-region sum of squares falls monotonically as regions are added, so a SKATER run at k = 30 will always look better than a Max-P run at p = 19. The comparison in the table above is only honest because k was set equal to the p that Max-P returned. If you must compare across different region counts, normalise — report the share of total sum of squares explained per region, or fix one algorithm’s output as the other’s input as done here.

Choose by what is actually fixed, not by which is fashionable

The rule is short. If the number of regions is fixed by the decision — a fixed number of teams, depots, strata or administrative units — use SKATER and then audit the sizes you got. If a minimum size is fixed by the data, by a disclosure rule, or by the standard error you can tolerate, use Max-P and accept whatever p it returns. The deeper treatment of the second case, including how to choose the threshold itself, is in Spatially Constrained Clustering with Max-P Regions.

Troubleshooting

Symptom	Likely cause	Fix
`Skater` returns fewer regions than `n_clusters`	The `floor` quorum makes k infeasible, or islands were merged in	Lower `floor`, or reduce `n_clusters`; check `len(w.islands)` first
Max-P raises during construction or returns p = 1	The threshold exceeds what any contiguous group can reach, or the graph is disconnected	Compare `threshold` against `gdf[col].sum() / w.n_components`; repair the graph before solving
Both partitions look like a map of one variable	Attributes never standardised	Fit `StandardScaler` on the attribute block and pass the scaled column names as `attrs_name`
Labels do not align with `gdf` rows	Weights built from a filtered or reindexed frame	Call `.reset_index(drop=True)` before building `w`, and pass `use_index=True`
Max-P p changes on every run	Global RNG not seeded, or the floor sits at a structural limit	`numpy.random.seed(s)` before each `solve()`; run ten seeds and report the modal p
SKATER regions are wildly uneven in population	Expected behaviour — SKATER has no size objective	Either accept it and report the spread, or switch to Max-P with an explicit floor
`Queen.from_dataframe` warns about disconnected components	Polygons that only touch at a shared node, or genuinely separate areas	Snap geometries, or solve each component separately and relabel

Next Steps

Once the partition is chosen, treat it as one zoning among many and measure how much your downstream statistics move when it changes — the method is in Spatial Scale & the Modifiable Areal Unit Problem. If the floor is the binding constraint in your work, go deeper on threshold selection and solver settings in Spatially Constrained Clustering with Max-P Regions.

Frequently Asked Questions

Which algorithm should I use if I need both a fixed number of regions and a minimum size?

Neither algorithm optimises both at once. Skater does accept a floor argument, but it counts areas per region rather than summing an attribute, so it cannot express “at least 40,000 people”. The workable approach is to run Max-P at your floor, read the p it returns, and treat that as the largest number of regions the floor permits. If you need fewer regions, run SKATER with k below p and check the resulting sizes by hand. If k must exceed p, the two requirements are simply infeasible on that contiguity graph.

Why does Max-P return a different number of regions each time I run it?

Max-P is a heuristic. Its construction phase grows candidate regions from randomly chosen seed areas, and the local search that follows is a randomised simulated-annealing swap. Different random states therefore give different feasible partitions, and sometimes a different p, because a partition that leaves one region marginally under the floor has to merge it away. Set numpy.random.seed() before solve() for reproducibility, then run ten seeds and report the modal p together with the spread of within-region sum of squares. A p that swings by more than one usually means the floor sits close to a structural limit.

Do I have to standardise the attributes before running either algorithm?

Yes, unless every attribute already shares a unit and a comparable spread. Both algorithms measure dissimilarity as a distance in attribute space, so a column measured in pounds swamps one measured in percentage points. In the example here median income carries over 99.999 per cent of the total sum of squares before scaling, and the partition is an income map with two decorative extra columns. Apply StandardScaler to the attribute block, keep the size column out of it, and confirm each scaled column has unit variance before solving.

How much do the spatial weights change the result?

More than most people expect, because the weights define which merges are possible at all rather than merely nudging an objective. Queen contiguity adds the corner-touching pairs that Rook omits, which on this 20 by 20 lattice raises the join count from 760 to 1,482 and gives the minimum spanning tree nearly twice as many candidate edges. Islands are worse: a disconnected graph makes the Max-P floor unsatisfiable for the stranded areas, and SKATER absorbs them quietly under islands='increase'. Build the weights once, inspect w.n_components and w.islands, and pass the same object to both models.

Related

Spatially Constrained Clustering with Max-P Regions — choosing the threshold and reading the solver’s diagnostics
Spatial Weight Matrices — the contiguity graph both algorithms are constrained by
Spatial Scale & the Modifiable Areal Unit Problem — why derived regions still have to be audited for zoning sensitivity

← Back to Spatial Clustering & Regionalization

SKATER vs Max-P Regionalization in Python

Why This Matters #

Environment and Version Pinning #

Step-by-Step Implementation #

1. Build the areal dataset #

2. Build one contiguity graph and reuse it #

3. Standardise the attributes #

4. Fit SKATER with a fixed k #

5. Fit Max-P with a fixed floor #

6. Compare the two partitions #

7. Check run-to-run stability #

Interpreting the Output #

Critical Best Practices #

Never run either algorithm on unstandardised attributes #

Set the seed before solve(), and report a distribution #

Match the weights to the geography, then keep them fixed #

Do not compare within-region sums of squares across different k #

Choose by what is actually fixed, not by which is fashionable #

Troubleshooting #

Next Steps #

Frequently Asked Questions #

Related

Why This Matters

Environment and Version Pinning

Step-by-Step Implementation

1. Build the areal dataset

2. Build one contiguity graph and reuse it

3. Standardise the attributes

4. Fit SKATER with a fixed k

5. Fit Max-P with a fixed floor

6. Compare the two partitions

7. Check run-to-run stability

Interpreting the Output

Critical Best Practices

Never run either algorithm on unstandardised attributes

Set the seed before `solve()`, and report a distribution

Match the weights to the geography, then keep them fixed

Do not compare within-region sums of squares across different k

Choose by what is actually fixed, not by which is fashionable

Troubleshooting

Next Steps

Frequently Asked Questions