Spatial Clustering & Regionalization

Regionalization is the art of drawing regions that are both internally homogeneous and geographically contiguous. Statistical agencies use it to build reporting zones with a minimum population; ecologists use it to delineate habitat regions; utilities use it to carve service territories from customer data. The defining constraint is contiguity: unlike ordinary clustering, which can scatter members of a group across the map, regionalization requires every region to be a single connected block of adjacent units. That constraint is enforced through the spatial weight matrix, which tells the algorithm which units are allowed to merge. This page develops the objective functions, the leading methods (Max-P-Regions, SKATER, and Ward with connectivity), and a runnable spopt implementation, all within the broader context of Core Concepts of Spatial Statistics & Geostatistics.

Prerequisites #

• Python 3.9+
• numpy>=1.22, geopandas>=0.13, libpysal>=4.8, spopt>=0.5, scikit-learn>=1.2
• A GeoDataFrame of areal units in a projected (metric) CRS with a reset integer index
• One or more numeric clustering variables, standardised to zero mean and unit variance
• A contiguity weight matrix with a single connected component — islands cannot join a region

Mathematical Core #

Regionalization partitions $n$ areal units into $p$ regions $\{R_1, \dots, R_p\}$ that minimise within-region heterogeneity subject to a contiguity constraint. Let $\mathbf{x}_i \in \mathbb{R}^d$ be the standardised attribute vector for unit $i$. A general objective is the total within-region sum of squared deviations from each region’s centroid:

$$\min \sum_{r=1}^{p} \sum_{i \in R_r} \lVert \mathbf{x}_i - \boldsymbol{\mu}_r \rVert^2, \qquad \boldsymbol{\mu}_r = \frac{1}{|R_r|} \sum_{i \in R_r} \mathbf{x}_i$$

subject to the constraint that every region $R_r$ induces a connected subgraph of the contiguity graph defined by the weight matrix $W$. It is this constraint that separates regionalization from ordinary k-means, which optimises the same objective with no spatial term.

Each symbol has a physical reading worth holding onto. $n$ is the number of areal units — counties, tracts, grid cells — and is fixed by the geography you were handed rather than chosen by the model. $p$ is the number of regions and $d$ the number of clustering variables, so $\mathbf{x}_i$ is one unit’s position in a $d$-dimensional attribute space after standardisation. The region centroid $\boldsymbol{\mu}_r$ is not a geographic centre: it is the mean attribute profile of the units assigned to region $r$, a point in attribute space that no unit need occupy. The squared distance $\lVert \mathbf{x}_i - \boldsymbol{\mu}_r \rVert^2$ therefore measures how atypical unit $i$ is of its own region, and the sum over all units is the total within-region sum of squares. For standardised variables that total starts at roughly $nd$ when every unit sits in a single region and falls towards zero as $p$ approaches $n$. Because it decreases monotonically with $p$, the objective can never select $p$ on its own; that number has to come from a floor constraint, an external design requirement, or the point at which extra regions stop buying a meaningful reduction.

SKATER: Minimum Spanning Tree Partitioning #

SKATER (Spatial 'K’luster Analysis by Tree Edge Removal) first builds a minimum spanning tree (MST) over the contiguity graph, where each edge $(i, j)$ carries a cost equal to the attribute dissimilarity between adjacent units:

$$c_{ij} = \lVert \mathbf{x}_i - \mathbf{x}_j \rVert$$

The MST reduces the full connectivity graph to $n - 1$ edges while preserving connectivity. SKATER then removes $p - 1$ edges to split the tree into $p$ connected subtrees, choosing at each step the edge whose removal yields the greatest reduction in within-region variance. Because every subtree is by construction a connected component of the original graph, each region is guaranteed contiguous.

Max-P-Regions: Maximising Regions Under a Floor #

Max-P-Regions inverts the question. Rather than fixing $p$, it fixes a minimum-size threshold and finds the maximum number of regions $p$ such that each region satisfies a floor constraint on a spatially extensive variable (population, area, count):

$$\sum_{i \in R_r} z_i \ge \text{floor} \quad \text{for every region } R_r$$

where $z_i$ is the threshold variable for unit $i$. Among all partitions satisfying the floor, Max-P selects the one that also minimises within-region attribute heterogeneity. This is the right formulation when the design requirement is “every region must contain at least 50,000 people” rather than “give me exactly eight regions”. Max-P is NP-hard and is solved with a heuristic (construction followed by local search).

Ward with Connectivity #

Agglomerative Ward clustering can be made spatial by supplying a connectivity matrix derived from $W$. At each merge step Ward joins the two adjacent clusters whose union produces the smallest increase in within-cluster variance, and the connectivity constraint restricts merges to spatially adjacent clusters. It produces a full hierarchy (a dendrogram) rather than a single partition, which is useful when you want to explore regionalizations at multiple scales.

The constraint changes Ward’s behaviour in a way worth anticipating. Unconstrained Ward merges the globally cheapest pair at every step, so merge heights rise monotonically and the dendrogram reads as a calibrated distance scale. Once connectivity is supplied, the cheapest available merge is often forbidden because the two clusters do not touch, and the algorithm is forced to take a costlier adjacent merge instead. Merge heights can then fall from one step to the next — a dendrogram inversion — and the familiar habit of cutting the tree at a large vertical gap becomes unreliable. Cut by requested cluster count instead, and read the dendrogram as a nesting structure rather than as a measure of how far apart the merged groups really are.

Annotated Implementation #

The spopt.region module implements SKATER and Max-P-Regions against a libpysal weights object. The pattern is: standardise attributes, build contiguity weights, fit the solver, extract labels.

1. Attribute Standardisation and Contiguity Weights #

import geopandas as gpd
import numpy as np
import libpysal
from sklearn.preprocessing import StandardScaler

# Load areal units in a projected CRS and reset the index
gdf = gpd.read_file("counties.gpkg").to_crs("EPSG:5070")
gdf = gdf.reset_index(drop=True)

# Clustering variables — standardise so no feature dominates the metric
attrs = ["median_income", "pct_college", "unemployment"]
gdf[attrs] = StandardScaler().fit_transform(gdf[attrs].values)

# Queen contiguity defines which units may merge into a shared region
w = libpysal.weights.Queen.from_dataframe(gdf)
assert w.n_components == 1, "Weight matrix must be a single connected component"

The n_components == 1 assertion is critical: regionalization cannot bridge disconnected components, so islands must be resolved (via KNN or manual bridging) before fitting.

2. SKATER for a Fixed Number of Regions #

from spopt.region import Skater

skater = Skater(
    gdf,
    w,
    attrs_name=attrs,
    n_clusters=6,          # target number of contiguous regions
    floor=5,               # minimum units per region
    trace=False,
)
skater.solve()

gdf["skater_region"] = skater.labels_
print(gdf["skater_region"].value_counts())

Skater takes the GeoDataFrame, the weights, the attribute names, and n_clusters. After solve(), labels_ holds the region index for each unit.

3. Max-P-Regions Under a Size Floor #

from spopt.region import MaxPHeuristic

# A spatially extensive threshold variable — e.g. population per unit
gdf["pop"] = gdf["population"].astype(float)

maxp = MaxPHeuristic(
    gdf,
    w,
    attrs_name=attrs,      # homogeneity variables
    threshold_name="pop",  # the floor is applied to this column
    threshold=50000,       # every region must contain >= 50,000 people
    top_n=2,
)
maxp.solve()

gdf["maxp_region"] = maxp.labels_
print(f"Max-P found {maxp.p} regions")

MaxPHeuristic discovers the number of regions rather than taking it as input. maxp.p reports how many regions satisfied the floor.

Output Interpretation #

Evaluate a regionalization on two axes: contiguity (are regions connected?) and homogeneity (are units within a region similar?).

def within_region_variance(gdf, label_col, attrs):
    """Mean within-region variance across clustering variables (lower = tighter)."""
    total = 0.0
    for _, grp in gdf.groupby(label_col):
        total += grp[attrs].var(ddof=0).sum() * len(grp)
    return total / len(gdf)

wrv = within_region_variance(gdf, "skater_region", attrs)
print(f"Mean within-region variance: {wrv:.3f}")

Good signs: each region label maps to a solid connected block, and within-region variance is well below the global variance of the standardised attributes (which is $d$, the number of variables, since each has unit variance). Warning signs: a region consisting of a single unit (the floor was too low, or an outlier could not merge), or within-region variance close to the global value, which means the contiguity constraint has overwhelmed attribute similarity — often a sign that the attributes carry little spatial structure to begin with. Whether that structure exists is exactly what the spatial autocorrelation metrics quantify: variables with near-zero Moran’s I have no coherent regions to find.

Production Considerations #

Complexity. Max-P-Regions is NP-hard; the heuristic’s runtime grows with the number of units and the tightness of the floor. For $n$ in the tens of thousands, expect minutes and set top_n modestly. SKATER is cheaper because the MST construction is near-linear in the number of edges and edge removal is $O(p)$.

Determinism. Both heuristics involve randomised initialisation. Pass a fixed random seed where the API exposes one, and record it, so a reported partition can be reproduced.

Islands. Resolve disconnected components before fitting. A single island forces the solver either to fail or to emit a singleton region; neither is usually acceptable. Repair geometry, switch to KNN, or bridge the island manually — the same remedies covered in the spatial weight matrices guide.

Scaling variables. Never skip standardisation. Because the objective is a variance across features, an unstandardised high-magnitude variable will silently dominate every region boundary.

Weak spatial structure. Every method here assumes the clustering variables vary smoothly across space. As that assumption weakens the failure is quiet rather than loud: the solver still returns perfectly contiguous regions, because contiguity is a hard constraint and not a fitted quantity, but the boundaries increasingly trace the shape of the adjacency graph instead of anything in the data. The measurable symptom is within-region variance sitting close to the global value $d$, and the pre-flight check is the univariate Moran’s I of each clustering variable. A variable near zero adds noise to the objective without adding signal, and is better dropped than standardised alongside the rest.

Troubleshooting #

Next Steps #

For a focused, copy-pasteable walkthrough of the Max-P heuristic — passing a libpysal weights object, setting the threshold variable and floor, extracting labels, and mapping the regions — follow Spatially Constrained Clustering with Max-P Regions. Because every method here depends on the connectivity graph, revisit spatial weight matrices to get the contiguity structure right, and use the spatial autocorrelation metrics to confirm your variables have the spatial structure worth regionalizing in the first place.

Frequently Asked Questions #

What is the difference between clustering and regionalization?

Attribute-space clustering such as k-means groups observations by feature similarity alone and ignores geography, so a resulting cluster can be scattered across the map. Regionalization adds a contiguity constraint through a spatial weight matrix, requiring every region to form a single connected block of adjacent units. Regionalization trades some attribute homogeneity for spatial contiguity, which is what makes the output usable as administrative or planning regions.

When should I use Max-P-Regions instead of SKATER?

Use SKATER when you know how many regions you want; it partitions a minimum spanning tree into a fixed number of contiguous groups. Use Max-P-Regions when you instead have a minimum-size requirement, such as a minimum population per region, and want the algorithm to find the maximum number of regions that all satisfy that floor while remaining internally homogeneous.

Why must clustering variables be standardised?

Regionalization objectives are built on a distance or variance metric across the clustering variables. If one variable is measured in the thousands and another in fractions, the large-magnitude variable dominates the metric and the regions form almost entirely along that single dimension. Standardising each variable to zero mean and unit variance gives every feature comparable influence on the partition.

Does regionalization guarantee contiguous regions?

Yes, provided the spatial weight matrix is a single connected component. SKATER and Max-P-Regions operate on the connectivity graph and can only merge units that are adjacent, so every output region is contiguous. If the weight matrix contains islands or disconnected components, those units cannot join a region and the solver will fail or isolate them.

How do I choose the number of regions for SKATER?

The within-region sum of squares falls monotonically as p rises, so the objective itself cannot pick p for you. Plot mean within-region variance against p across a range of candidate values and look for the point where the curve flattens, then check that candidate against the practical constraint: the number of reporting zones you need, the crews you can staff, or the minimum unit count implied by the floor. If no flattening appears, the variables probably lack the spatial structure that would make any particular p natural, and Max-P-Regions with an explicit size floor is the more defensible formulation.

Related

• Spatially Constrained Clustering with Max-P Regions — the step-by-step spopt MaxPHeuristic workflow
• Spatial Weight Matrices — the contiguity graph that enforces the region constraint
• Spatial Autocorrelation Metrics — confirming variables carry the spatial structure worth regionalizing
• Hot Spot Analysis with Getis-Ord Statistics — locating intense clusters as a complement to partitioning

← Back to Core Concepts of Spatial Statistics & Geostatistics