--- title: "Geodemographic classification in Python: k-means on Census 2021 Output Areas, and SKATER and max-p regionalisation" description: "Build a Census 2021 Output Area classification in Python: k-means with silhouette and stability checks, index scores, and SKATER and max-p in spopt." url: https://www.technicalgis.com/tutorials/geodemographic-classification-python/ author: TechnicalGIS Team publisher: TechnicalGIS (https://www.technicalgis.com/about/) published: 2026-10-27 updated: 2026-10-27 language: en-GB --- # Geodemographic classification in Python: k-means on Census 2021 Output Areas, and SKATER and max-p regionalisation > Build a Census 2021 Output Area classification in Python: k-means with silhouette and stability checks, index scores, and SKATER and max-p in spopt. - Level: Intermediate - Tested with: Python 3.12 - Source code: [Download scripts](https://www.technicalgis.com/code/geodemographic-classification-python/) ## Key takeaways 1. A geodemographic classification groups small areas by the make-up of their populations, so that areas of the same type can be far apart. Here, 25 Census 2021 variables for Cardiff's 1,115 Output Areas, transformed with the inverse hyperbolic sine and standardised, gave seven k-means clusters, from social-rented family estates to a student quarter. 2. Variable selection matters more than the algorithm. Four pairs of variables correlated at |r| ≥ 0.8, and the tenure shares sum to 100%, so two variables were dropped to stop one theme from counting twice. 3. There is rarely a 'correct' k. The elbow curve had no clear elbow, the mean silhouette peaked at k = 3 (0.241), and single-start k-means was most stable at k = 5. We chose k = 7 because it separates students and city-centre renters, which are distinct in Cardiff, and used 100 random starts because single starts at k = 7 agreed with the best solution only moderately (median ARI 0.708). 4. Profile clusters with index scores (cluster rate ÷ city rate × 100) and name them cautiously. Names are summaries for discussion; every cluster contains areas that do not fit its name, and 5.3% of OAs had a negative silhouette. 5. Adding a contiguity constraint costs homogeneity. The seven k-means clusters fall into 157 separate pieces and explain 52% of the variance; seven SKATER regions explain 34%, and 16 max-p regions of at least 20,000 residents explain 46%, against 61% for k-means with k = 16. Use contiguity when you need areas, not types. ## Before you start - Comfortable with pandas and GeoPandas (see the GeoPandas getting-started tutorial) - Python 3.10+ with scikit-learn, spopt, libpysal, networkx, geopandas, numpy, pandas, scipy, matplotlib and requests - An internet connection: the first run downloads about 57 MB of Census 2021 bulk tables from Nomis, 10 MB for the official OAC and the Cardiff Output Area boundaries from the ONS **Short answer:** to build a geodemographic classification in Python, download Census 2021 tables for Output Areas from Nomis, turn counts into percentages, transform skewed variables (the inverse hyperbolic sine, `np.arcsinh`, works well), standardise them to z-scores, drop near-duplicate variables, and run scikit-learn’s `KMeans` with many random starts. Choose k from the elbow, the silhouette and stability across seeds, together with how useful the clusters are, then profile them with index scores before you name them. If you need **contiguous regions** rather than area types, use `spopt`’s `Skater` or `MaxPHeuristic` on a queen contiguity graph. For Cardiff’s 1,115 Output Areas, **seven k-means clusters explained 52% of the variance in 25 variables but fell into 157 separate pieces; seven SKATER regions explained 34%**. Geodemographic classifications sort neighbourhoods into types such as “student quarter” or “detached owner-occupier suburbs”. The UK has an official one, the [2021/2 Output Area Classification (OAC)](https://data.geods.ac.uk/dataset/73a16168-f2b0-4732-b36a-3ec0253e34fe), and commercial ones such as Acorn and Mosaic. Building your own is useful when you want one city, different variables or a different purpose, and doing it once is the best way to learn how much the official product depends on judgement. This tutorial builds a small classification of Cardiff with real Census 2021 data and is explicit about every judgement it makes. You will learn how to: 1. download **Census 2021 tables for Output Areas** from Nomis and build 27 candidate variables; 2. **transform and standardise** them, and **check correlations** to remove duplicates; 3. **choose k** for k-means with the **elbow**, the **silhouette** and **stability across random seeds**; 4. **profile the clusters** with index scores, a heatmap and radar charts, and **name them cautiously**; 5. **compare** the result with the **official 2021/2 OAC** and see the **MAUP** at work at LSOA level; 6. build **contiguous regions** with **SKATER** and **max-p** in spopt, and decide **when you want contiguity**. ## Setup The full script, [`geodemographic_classification.py`](https://www.technicalgis.com/code/geodemographic-classification-python/), runs from top to bottom and writes seven figures, a GeoPackage and a JSON file of metrics to `output/`. It caches downloads in `data/`. Install the tested versions with: ```bash pip install scikit-learn==1.9.1 spopt==0.8.0 libpysal==4.15.1 networkx==3.7 geopandas==1.2.0 shapely==2.2.0 pyogrio==0.13.0 pyproj==3.8.0 numpy==2.5.4 pandas==3.1.0 scipy==1.18.1 matplotlib==3.11.2 requests==2.34.2 ``` The whole script takes about eight minutes on a laptop, almost all of it in max-p (Step 8). Everything else takes about a minute. ## Step 1: Output Areas for Cardiff Output Areas (OAs) are the smallest areas with published Census 2021 results. Cardiff has 1,115 of them, grouped into 218 LSOAs. The ONS Open Geography Portal serves the boundaries through an ArcGIS REST service; LSOA names start with the local authority name, so a `LIKE` filter selects Cardiff. `geodemographic_classification.py` ```python OA_URL = f"{ONS}/Output_Areas_2021_EW_BGC_V2/FeatureServer/0/query" feats, offset = [], 0 while True: page = requests.get(OA_URL, params={"where": "LSOA21NM LIKE 'Cardiff %'", "outFields": "OA21CD,LSOA21CD,LSOA21NM", "outSR": 27700, "f": "geojson", "resultOffset": offset, "resultRecordCount": 1000}, timeout=300).json() feats += page["features"]; offset += len(page["features"]) if len(page["features"]) < 1000: break ``` `Output — Output Areas` ```text 1,115 Output Areas in 218 LSOAs, 142.3 km2 ``` ## Step 2: Census 2021 variables from Nomis Nomis publishes every Census 2021 topic summary table as a bulk zip with one CSV per geography level, so a single download gives you all 188,880 OAs in England and Wales. The script downloads 15 tables, keeps Cardiff’s OAs, and turns counts into percentages of the right base: people for age and ethnicity, households for tenure and cars, residents aged 16 and over for occupation and qualifications. `geodemographic_classification.py` ```python NOMIS = "https://www.nomisweb.co.uk/output/census/2021/census2021-{}.zip" def table(code): z = zipfile.ZipFile(download(NOMIS.format(code), DATA / f"census2021-{code}.zip")) df = pd.read_csv(z.open(f"census2021-{code}-oa.csv")).set_index("geography code") df.columns = [c.split(": ", 1)[1].replace("; measures: Value", "") if ": " in c else c for c in df.columns] return df.loc[oa.OA21CD] def pct(df, num, den): num = [num] if isinstance(num, str) else num return 100 * df[num].sum(axis=1) / df[den] age, hh, ten = table("ts007a"), table("ts003"), table("ts054") # ... 12 more tables V = pd.DataFrame(index=oa.OA21CD) V["Aged 65+"] = pct(age, ["Aged 65 to 69 years", "Aged 70 to 74 years", "Aged 75 to 79 years", "Aged 80 to 84 years", "Aged 85 years and over"], "Total") V["Social rented"] = pct(ten, "Social rented", "Total: All households") V["Full-time students"] = pct(sec, "L15 Full-time students", sec_tot) V["Born outside UK"] = 100 - pct(cob, "Europe: United Kingdom", "Total") V["Welsh speakers"] = pct(welsh, [c for c in welsh.columns if c.startswith("Can speak")], "Total: All usual residents aged 3 years and over") # ... 22 more variables ``` The 27 candidates cover the themes that most classifications use: age, household type, tenure, housing type, overcrowding, cars, occupation (NS-SeC), qualifications, unemployment, health, ethnicity, country of birth, plus population density and, because this is Wales, the share of people who can speak Welsh. `Output — Census 2021 variables (excerpt)` ```text 27 candidate variables for 1,115 OAs; 0 missing values Cardiff: 362,289 usual residents; OA population median 313, range 101-1730 mean 50% min max Aged 20-24 9.3 5.6 0.0 79.4 Aged 65+ 15.4 13.6 0.0 91.1 Population density 7135.5 5802.5 53.9 47887.3 Social rented 17.1 7.2 0.0 91.1 Detached 14.4 4.3 0.0 98.4 Flats 25.5 15.1 0.0 100.0 Full-time students 12.4 7.4 0.0 95.2 Unemployed 2.9 2.6 0.0 37.2 Black 3.7 2.2 0.0 51.0 Welsh speakers 12.3 11.1 0.0 47.0 ``` ## Step 3: transform, standardise and check correlations The variables are on very different scales (density runs to 47,887 people per km², the percentages to 100) and many are skewed: most OAs have few students or Black residents, and a few have a lot. Without treatment, k-means would be dominated by density and by a handful of extreme OAs. The script applies the inverse hyperbolic sine, which behaves like a logarithm but is defined at zero, then converts every variable to a z-score (mean 0, standard deviation 1). `geodemographic_classification.py` ```python T = np.arcsinh(V) # inverse hyperbolic sine: log-like, but defined at zero Z = (T - T.mean()) / T.std(ddof=0) # z-scores C = Z.corr() pairs = [(a, b, C.loc[a, b]) for i, a in enumerate(C) for b in C.columns[i + 1:] if abs(C.loc[a, b]) >= 0.8] DROP = ["Degree level (L4+)", "Owns with mortgage"] Zk = Z.drop(columns=DROP) X = Zk.to_numpy() ``` `Output — transformation and correlation` ```text median |skewness|: raw 1.38 -> IHS 0.59; most skewed raw: Unemployed 4.4, Black 3.5, Aged 20-24 3.4 |r| >= 0.8: Higher professional ~ Degree level (L4+): r = 0.88 |r| >= 0.8: Routine occupations ~ No qualifications: r = 0.87 |r| >= 0.8: No qualifications ~ Bad health: r = 0.81 |r| >= 0.8: Aged 20-24 ~ Full-time students: r = 0.80 dropping ['Degree level (L4+)', 'Owns with mortgage']: degree level duplicates higher professional (r = 0.88); the four tenure shares sum to ~100%, so one is redundant 25 variables go into the clustering ``` ![A 27 by 27 correlation matrix of Census 2021 variables across Cardiff's Output Areas, coloured from dark blue for r = -1 to dark red for r = +1. Four pairs are outlined in black because |r| is at least 0.8: higher professional and degree level, routine occupations and no qualifications, no qualifications and bad health, and aged 20-24 and full-time students. Strong negative correlations appear between detached housing and no car, and between social renting and higher professional occupations.](https://www.technicalgis.com/_astro/01-correlation-matrix.BTeR9OYE_lqU15.webp) _Correlations between the 27 transformed, standardised variables. The boxed pairs have |r| ≥ 0.8. Strongly correlated variables measure the same thing twice and give that theme extra weight in k-means._ Correlation of 0.8 or more is a common rule of thumb for “measuring the same thing”. We dropped degree-level qualifications, which is almost the same signal as higher professional occupations, and owning with a mortgage, because the four tenure shares add up to 100% and any one is implied by the others. We kept the other three pairs: routine occupations, no qualifications and bad health overlap but are not the same concept, and students and 20–24-year-olds diverge in places with young workers. These are judgements, and a different analyst could defend different ones. ## Step 4: choosing k k-means needs the number of clusters in advance. Three diagnostics help, and none of them decides on its own: - the **elbow**: the within-cluster sum of squares falls as k grows; look for a bend after which extra clusters help little; - the **mean silhouette**: how much closer areas are to their own cluster than to the next nearest one; - **stability**: whether k-means finds the same solution from different random starts, measured with the adjusted Rand index (ARI). `geodemographic_classification.py` ```python for k in range(2, 13): km = KMeans(k, n_init=50, random_state=SEED).fit(X) print(k, km.inertia_, silhouette_score(X, km.labels_)) for k in (4, 5, 6, 7, 8, 9): runs = [KMeans(k, n_init=1, random_state=s).fit(X).labels_ for s in range(30)] best = KMeans(k, n_init=100, random_state=SEED).fit(X).labels_ vs_best = np.median([adjusted_rand_score(best, r) for r in runs]) pairwise = np.mean([adjusted_rand_score(runs[i], runs[j]) for i in range(30) for j in range(i + 1, 30)]) ``` `Output — choosing k` ```text k= 2: within-cluster sum of squares 21192 | silhouette 0.224 k= 3: within-cluster sum of squares 17401 | silhouette 0.241 k= 4: within-cluster sum of squares 15966 | silhouette 0.217 k= 5: within-cluster sum of squares 14763 | silhouette 0.169 k= 6: within-cluster sum of squares 13916 | silhouette 0.154 k= 7: within-cluster sum of squares 13286 | silhouette 0.161 k= 8: within-cluster sum of squares 12873 | silhouette 0.144 k= 9: within-cluster sum of squares 12474 | silhouette 0.128 k=10: within-cluster sum of squares 12152 | silhouette 0.130 k=11: within-cluster sum of squares 11856 | silhouette 0.122 k=12: within-cluster sum of squares 11550 | silhouette 0.116 k=4: 30 single-start runs vs the 100-start solution, median ARI 0.416; mean pairwise ARI 0.787 k=5: 30 single-start runs vs the 100-start solution, median ARI 0.959; mean pairwise ARI 0.833 k=6: 30 single-start runs vs the 100-start solution, median ARI 0.728; mean pairwise ARI 0.698 k=7: 30 single-start runs vs the 100-start solution, median ARI 0.708; mean pairwise ARI 0.710 k=8: 30 single-start runs vs the 100-start solution, median ARI 0.649; mean pairwise ARI 0.628 k=9: 30 single-start runs vs the 100-start solution, median ARI 0.641; mean pairwise ARI 0.613 chosen k = 7 ``` ![Three line charts against the number of clusters k. Left: the within-cluster sum of squares falls smoothly from about 21,200 at k = 2 to 11,550 at k = 12 with no sharp elbow. Middle: the mean silhouette peaks at 0.241 at k = 3, drops to 0.154 at k = 6, rises slightly to 0.161 at k = 7 and then falls to 0.116 at k = 12. Right: stability across 30 random seeds, measured by ARI, is highest at k = 5 (0.96 against the best solution) and about 0.7 at k = 7. A dotted vertical line marks the chosen k = 7 in each panel.](https://www.technicalgis.com/_astro/02-choosing-k.DeMmAY3I_1Furre.webp) _Three ways to choose k. The elbow is smooth, the silhouette prefers very few clusters, and k-means is most stable at k = 5. We chose k = 7 (dotted line) for interpretability, and used 100 random starts to get its best solution reliably._ This is a typical result for census data, and it is worth being honest about. There is no clear elbow. The silhouette is highest at k = 3, but three clusters (roughly “affluent suburbs”, “deprived estates” and “inner city”) are too coarse to be useful. Census variables vary continuously across a city, so clusters are convenient divisions of a continuum rather than natural groups, and silhouette values of 0.1–0.25 are normal. We chose **k = 7**. It has a small local peak in the silhouette and separates two groups that matter in Cardiff and merge at lower k: the student quarter and city-centre renters in flats. k = 5 was the most stable choice and would be a defensible alternative. Single-start runs at k = 7 agreed with the 100-start solution only moderately (median ARI 0.708), which is why the final model uses `n_init=100`. Always use many starts for the final run. ## Step 5: profile the clusters and name them cautiously The final model takes the best of 100 starts. To describe each cluster in units people understand, the script computes **index scores** from the untransformed percentages: the cluster’s rate (weighted by residents) divided by the Cardiff rate, times 100. `geodemographic_classification.py` ```python km = KMeans(7, n_init=100, random_state=SEED).fit(X) s = silhouette_samples(X, oa.kmeans) city = V.mul(residents, axis=0).sum() / residents.sum() idx = pd.DataFrame({c: V[oa.kmeans.to_numpy() == c].mul(residents[oa.kmeans == c], axis=0).sum() / residents[oa.kmeans == c].sum() for c in range(K)}).T index = (100 * idx / city).round(0) ``` `Output — final k-means classification` ```text k = 7: mean silhouette 0.161; OAs with negative silhouette 5.3% cluster 0 (250 OAs, 82,074 people): high Social rented 252, Lone parents 220, No qualifications 178, Routine occupations 171 | low Detached 33, Higher professional 37, Degree level (L4+) 54 cluster 1 (241 OAs, 77,080 people): high Detached 301, Owns outright 174, Higher professional 168, Aged 65+ 159 | low Social rented 11, Black 19, Flats 21 cluster 2 (188 OAs, 53,202 people): high Aged 65+ 160, Owns outright 137, Bad health 124, One-person households 116 | low Black 33, Full-time students 36, Aged 20-24 42 cluster 3 (160 OAs, 54,194 people): high Population density 294, Black 244, Asian 211, Born outside UK 190 | low Detached 24, Aged 65+ 56, Owns outright 59 cluster 4 (160 OAs, 50,693 people): high Terraced 162, Higher professional 150, Welsh speakers 142, Degree level (L4+) 136 | low Social rented 33, Overcrowded 51, Black 56 cluster 5 (58 OAs, 15,124 people): high Flats 348, Private rented 253, Born outside UK 195, Aged 25-44 187 | low Detached 4, Lone parents 24, Social rented 30 cluster 6 (58 OAs, 29,922 people): high Full-time students 506, Aged 20-24 459, Population density 287, Private rented 253 | low Detached 17, Aged 0-14 19, Aged 65+ 21 ``` The console numbers clusters from 0, ordered by size; the figures and the names add one. The heatmap shows all 27 index scores, including the two variables that were not used in the clustering. ![A heatmap with the seven named clusters as rows and 27 census variables as columns, each cell showing an index score against the Cardiff average of 100, coloured purple below 100 and orange above. Notable values: social-rented family estates 252 for social rented and 220 for lone parents; detached owner-occupier suburbs 301 for detached and 11 for social rented; dense multicultural inner city 294 for population density and 244 for Black residents; city-centre renters in flats 348 for flats; student quarter 506 for full-time students and 459 for aged 20-24.](https://www.technicalgis.com/_astro/04-index-scores.BfAN49tP_Z1Q4W88.webp) _Index scores for every variable and cluster (Cardiff = 100). Orange cells are above the city average, purple cells below; the scale is logarithmic, so 200 and 50 are equally far from 100._ ![Seven radar charts, one per cluster, each with 13 spokes for variables such as aged 0-14, aged 65+, owns outright, social rented, private rented, flats, no car, higher professional, routine occupations, full-time students, born outside UK and Welsh speakers. A dashed ring marks the Cardiff average. The student quarter has a long spike for full-time students and aged 20-24; the detached owner-occupier suburbs extend towards owns outright, aged 65+ and higher professional; the social-rented family estates extend towards social rented, routine occupations and no car.](https://www.technicalgis.com/_astro/05-radar-profiles.BLOXoJXS_ZvFNr0.webp) _The same profiles as radar charts of mean z-scores for 13 key variables (dashed ring = Cardiff average). Radar charts are good for spotting a cluster's shape at a glance, but the heatmap is easier to compare across clusters._ We named the clusters only after reading these profiles: 1. **Social-rented family estates** (250 OAs): 2.5 times the city’s share of social renting and 2.2 times its share of lone-parent households. 2. **Detached owner-occupier suburbs** (241 OAs): three times the share of detached houses; very little social renting. 3. **Older settled homeowners** (188 OAs): more people aged 65+, outright owners and people in bad health. 4. **Dense multicultural inner city** (160 OAs): about three times the city’s density and twice its shares of Black and Asian residents. 5. **Professional terraces** (160 OAs): terraced housing, professional occupations and the highest share of Welsh speakers. 6. **City-centre renters in flats** (58 OAs): 3.5 times the share of flats, 2.5 times the share of private renting, many people aged 25–44. 7. **Student quarter** (58 OAs, almost 30,000 people): five times the city’s share of full-time students. ![A map of Cardiff's Output Areas coloured by seven k-means clusters. The detached owner-occupier suburbs cover most of the north and west. Social-rented family estates form large blocks in the east (around Llanrumney and Trowbridge) and the south-west (Ely and Caerau). The student quarter is a compact area north of the city centre around Cathays, the city-centre renters in flats run south to Cardiff Bay, and the dense multicultural inner city lies around them in Riverside, Grangetown, Butetown and Adamsdown. A 2 km scale bar is in the lower left.](https://www.technicalgis.com/_astro/03-kmeans-map.CtTc_jtO_Z1TYW28.webp) _The seven-cluster classification of Cardiff's 1,115 Output Areas. Types repeat across the city: social-rented estates appear in both the east and the south-west._ The map shows the main property of a classification: areas of the same type need not be next to each other. Social-rented family estates appear in the east and the south-west, and professional terraces appear in several inner suburbs. The silhouette also tells you which areas fit badly: 5.3% of OAs are closer, on average, to another cluster than to their own. They are worth listing and checking on the ground before anyone uses the classification for decisions. ## Step 6: compare with the official 2021/2 OAC The UK has an official classification. The **2021/2 Output Area Classification**, built by UCL and the ONS and published by the Geographic Data Service, assigns every OA in the UK to one of 8 supergroups, 21 groups and 52 subgroups using 60 census variables, the inverse hyperbolic sine and k-means: broadly the same recipe as this tutorial, at national scale. It is the reference product to use when you need a classification that is comparable across the UK. The script downloads it and cross-tabulates our clusters against its supergroups. `geodemographic_classification.py` ```python oac = pd.read_csv(download(OAC_URL, DATA / "uk_oac_final.csv"), dtype=str).set_index("GeographyCode") oa["oac_sg"] = oac.Supergroup.reindex(oa.OA21CD).to_numpy() print(adjusted_rand_score(oa.oac_sg, oa.kmeans)) print(pd.crosstab(oa.kmeans, oa.oac_sg)) ``` `Output — comparison with the official 2021/2 OAC` ```text matched 1115 of 1115 OAs to the official classification official supergroups in Cardiff: 1 Retired Professionals 111; 2 Suburbanites and Peri-Urbanites 113; 3 Multicultural and Educated Urbanites 139; 4 Low-Skilled Migrant and Student Communities 141; 5 Ethnically Diverse Suburb-an Professionals 182; 6 Baseline UK 374; 7 Semi- and Un-Skilled Workforce 44; 8 Legacy Communities 11 ARI between our k-means and the official supergroups (8 classes, UK-wide): 0.282 ARI against the official groups (21 present): 0.323 ``` ![A cross-tabulation heatmap with our seven Cardiff clusters as rows and the eight official OAC supergroups as columns, with the number of Output Areas in each cell and darker blue for a larger share of the row. City-centre renters in flats are almost all Multicultural and Educated Urbanites (53 of 58). Social-rented family estates are mostly Baseline UK (183 of 250). Detached owner-occupier suburbs split between Retired Professionals (106), Ethnically Diverse Suburb-an Professionals (101) and Suburbanites and Peri-Urbanites (34). The dense multicultural inner city is mostly Low-Skilled Migrant and Student Communities (93 of 160).](https://www.technicalgis.com/_astro/07-official-oac-crosstab.BpgTCjfG_1006gH.webp) _Our Cardiff clusters against the official 2021/2 OAC supergroups. Numbers are OAs; shading is the share of each row. Supergroup names are the official ones._ The two classifications agree in broad terms but are far from identical (ARI 0.282 against supergroups, 0.323 against groups). That is expected. The OAC was built for the whole UK, so its clusters describe differences between places such as rural Wales, inner London and Cardiff; within one city, a UK-wide type such as “Baseline UK” (374 of Cardiff’s OAs) can hide differences that matter locally. Our classification was built only on Cardiff, so its types are relative to Cardiff. Some differences are also instructive: our “Student quarter” is mostly in the OAC’s “Multicultural and Educated Urbanites” supergroup (38 of 58 OAs) rather than its “Low-Skilled Migrant and Student Communities”, which in Cardiff matches our dense inner-city cluster better. Neither classification is wrong; they answer different questions. ## Step 7: the MAUP: what changes at LSOA level? Geodemographic results depend on the zones you use, which is one form of the **modifiable areal unit problem (MAUP)**. To see how much, the script aggregates the same variables to Cardiff’s 218 LSOAs (population-weighted means of the OA percentages, which approximate the LSOA rates) and repeats the whole recipe. `geodemographic_classification.py` ```python VL = V.drop(columns="Population density").mul(residents, axis=0).groupby(lsoa_key).sum().div( pd.Series(residents).groupby(lsoa_key).sum().to_numpy(), axis=0) TL = np.arcsinh(VL[V.columns].drop(columns=DROP)) ZL = ((TL - TL.mean()) / TL.std(ddof=0)).to_numpy() lab_l = KMeans(7, n_init=100, random_state=SEED).fit(ZL).labels_ oa_from_lsoa = pd.Series(lab_l, index=VL.index).reindex(lsoa_key).to_numpy() ``` `Output — MAUP at LSOA level` ```text 218 LSOAs, k = 7: silhouette 0.212; ARI with the OA classification (each OA takes its LSOA's cluster) 0.416 LSOAs whose OAs fall in more than one OA-level cluster: 70%; max-min range of 'Full-time students' within LSOAs, median 5.9 points ``` The LSOA classification has a higher silhouette (0.212 against 0.161), but that does not make it better: averaging over about five OAs smooths away the variation that makes OAs hard to separate. The two classifications agree only moderately (ARI 0.416), and 70% of LSOAs contain OAs from more than one OA-level cluster. Use the smallest zones your data support, report the zone you used, and do not apply a classification built for one geography to another without checking. ## Step 8: contiguous regions with SKATER and max-p A classification answers “what kind of place is this?” A **regionalisation** answers a different question: “how should I divide this city into a few internally similar, connected regions?” Each region must be a single contiguous piece. You want contiguity when the result will be used as a set of areas, for example: - **reporting or service areas**, such as neighbourhood management areas, catchments or sales territories, which must be contiguous to be managed; - **zones for small-area statistics**, where max-p’s population floor guarantees that every region has enough people for stable rates and for disclosure control; - **sampling frames and survey strata** that fieldworkers can cover area by area; - **spatial analysis that needs areas of similar size**, to reduce the effect of very small or very large zones. You do not want contiguity when you are targeting a type of area (a leaflet campaign to every student area, for instance) or comparing like with like across a city: a classification lets similar areas be far apart, and a contiguity constraint would force unlike areas together. spopt, the regionalisation module of PySAL, provides SKATER, max-p, AZP and others. They all need a **spatial weights** object that says which OAs touch. The script uses queen contiguity (sharing an edge or a corner) and checks that the graph is connected, because an island with no neighbours can never join a region. `geodemographic_classification.py` ```python w = libpysal.weights.Queen.from_dataframe(oa, use_index=False, silence_warnings=True) attrs = list(Zk.columns) reg = oa[["geometry"]].copy() reg[attrs] = X # the same 25 standardised variables as k-means reg["residents"] = residents sk = Skater(reg, w, attrs, n_clusters=7, floor=20, trace=False, islands="increase") sk.solve() np.random.seed(SEED); random.seed(SEED) # max-p's heuristic is randomised mp = MaxPHeuristic(reg, w, attrs, "residents", 20_000, top_n=2, max_iterations_construction=50) mp.solve() ``` **SKATER** builds a minimum spanning tree of the contiguity graph, with edge costs equal to the attribute dissimilarity between neighbours, and then cuts the tree into k pieces, choosing the cuts that most reduce within-region variation. Here we ask for seven regions with at least 20 OAs each. **max-p** works the other way round: you set a floor on a spatially extensive variable (here, at least 20,000 residents per region) and it finds the **largest number of regions** that meets the floor, then improves their homogeneity by moving OAs between regions. AZP (`spopt.region.AZP`) is a similar local-search alternative when you want a fixed number of regions. `Output — contiguity, SKATER and max-p` ```text queen contiguity: mean 6.0 neighbours, islands 0, connected components 1 SKATER with 7 regions, floor 20 OAs: 2.0 s; region sizes (OAs) [399, 242, 141, 129, 108, 53, 43] max-p, every region >= 20,000 residents: p = 16 regions in 411.0 s; region populations 20,025-26,734 ``` To compare the three results, the script computes the share of the variance in the 25 standardised variables that each grouping explains (1 − within-group sum of squares ÷ total sum of squares, an R²), and counts the contiguous pieces each grouping makes on the map. `Output — how much homogeneity does contiguity cost?` ```text k-means (no contiguity) : 7 groups | share of variance explained 0.523 | contiguous pieces on the map 157 SKATER : 7 groups | share of variance explained 0.344 | contiguous pieces on the map 7 max-p : 16 groups | share of variance explained 0.464 | contiguous pieces on the map 16 ARI k-means vs SKATER 0.266; k-means vs max-p 0.205 for reference, aspatial k-means with k = 16 explains 0.614 ``` ![Three maps of Cardiff. Left: the seven k-means clusters, which form a patchwork of 157 separate pieces. Middle: seven SKATER regions with black boundaries, including one very large region covering most of the north and west, a large region in the east, and smaller regions around the city centre and the south-west. Right: 16 max-p regions of at least 20,000 residents each, smaller and more even in size, with black boundaries. In the middle and right maps the colours only distinguish neighbouring regions.](https://www.technicalgis.com/_astro/06-skater-maxp.DI5xI5Uk_1X99LX.webp) _The same 25 variables grouped three ways. Left: k-means, with no spatial constraint (157 pieces, R² 0.52). Middle: SKATER, seven contiguous regions (R² 0.34). Right: max-p, 16 regions of at least 20,000 residents (R² 0.46)._ The contiguity constraint has a clear cost. With seven groups, SKATER explains 34% of the variance against 52% for k-means, because each region must include whatever lies between its similar areas. SKATER’s regions are also very unequal (399 OAs in the largest, 43 in the smallest), since it cuts the tree wherever that helps homogeneity most. max-p’s regions are more even by design (20,025–26,734 residents each), and its 16 regions explain 46% of the variance, against 61% for 16 unconstrained k-means clusters. In return you get areas you can draw, name and manage. ## Limitations and honest caveats - **Census 2021 was taken during the COVID-19 pandemic.** On 21 March 2021 some students were living at their family homes rather than at university, so student areas in Cardiff may be under-represented, and some people’s occupations and travel were affected by furlough and lockdown rules. - **The variables decide the result.** Different variables, or a different balance between themes, would give different clusters. The 27 candidates here are a reasonable small set, not a standard. - **k is a choice.** The diagnostics did not agree, and k = 5 would have been equally defensible. Report the alternatives you considered. - **Clusters are averages.** The ecological fallacy applies: a cluster’s profile does not describe every person, household or even every OA in it, and 5.3% of OAs fit their cluster poorly. - **Index scores need context.** An index of 500 for a variable with a city rate of 1% means 5%; always show rates alongside index scores when it matters. - **Census statistics have small-number noise.** The ONS applies targeted record swapping and cell key perturbation to protect confidentiality, so very small counts in single OAs are less reliable than the totals. - **The MAUP does not go away.** The results depend on the OA boundaries; a classification of LSOAs or of a grid would differ. - **Regionalisation results depend on the parameters.** SKATER’s floor and number of regions, and max-p’s threshold, are choices; so is queen rather than rook contiguity. ## Frequently asked questions ### How do I do k-means clustering on census data in Python? Turn counts into percentages of a suitable base, transform skewed variables (for example with `np.arcsinh`), standardise every variable to z-scores, drop near-duplicates, and fit `sklearn.cluster.KMeans(k, n_init=100, random_state=42)`. Choose k with the elbow, the silhouette and stability across random seeds, then profile the clusters with index scores against the study-area average. ### Why standardise variables before k-means? k-means uses Euclidean distance, so a variable measured on a larger scale, such as population density in thousands per km², would dominate percentages that range from 0 to 100. Converting each variable to a z-score gives every variable the same weight. Transforming skewed variables first (log, square root or inverse hyperbolic sine) stops a few extreme areas from dominating the standard deviation. ### What is the difference between classification and regionalisation? A classification groups areas by similarity only, so areas of the same type can be anywhere in the city. A regionalisation adds a contiguity constraint, so each group is one connected region. Use a classification to describe and target types of neighbourhood, and a regionalisation when you need a set of areas to manage, report on or sample from. ### Should I use SKATER or max-p? Use SKATER when you know how many regions you want and need a fast result; it can produce very unequal regions. Use max-p when every region must reach a minimum size, such as a population floor for stable rates or disclosure control, and you want the largest number of regions that meets it; it is randomised and much slower. AZP is a middle ground for a fixed number of regions. ### Should I use the official OAC instead of building my own? For most purposes, yes. The 2021/2 OAC covers the whole UK consistently, is documented and peer-reviewed, and is free under the Open Government Licence. Build your own when you need types relative to one city, variables the OAC does not use (such as the Welsh language here) or a classification for a specific purpose, and compare it with the OAC as in Step 6. ### Can I use LSOAs instead of Output Areas? Yes, and LSOAs are the only option for some data, such as the indices of deprivation. Expect smoother, more mixed clusters: in Cardiff, 70% of LSOAs contained OAs from more than one OA-level cluster, and the LSOA and OA classifications agreed only moderately (ARI 0.416). ## Next steps - Cluster points rather than areas, for example shops or GPS fixes, with our companion tutorial [DBSCAN and HDBSCAN spatial clustering in Python](https://www.technicalgis.com/tutorials/dbscan-spatial-clustering-python/). - Build spatial weights in more depth and test whether a single variable is spatially clustered with [hot spot analysis in Python (Getis-Ord Gi\*)](https://www.technicalgis.com/tutorials/hot-spot-analysis-getis-ord-python/). - Relate two census variables to each other across small areas with [bivariate correlation of census and crime data](https://www.technicalgis.com/2022/08/14/performing-bivariate-correlation-analysis-using-census-and-crime-data/). - Attach the clusters to points such as addresses or incidents with [spatial joins in GeoPandas](https://www.technicalgis.com/tutorials/spatial-join-geopandas/). The script writes `output/cardiff_oa_classification.gpkg` with the k-means, OAC, SKATER and max-p labels for every OA. - Make publication-quality maps of the clusters with [maps in Python with Matplotlib and Cartopy](https://www.technicalgis.com/tutorials/python-maps-matplotlib-cartopy/). ## Data and licences - **Census 2021**: Office for National Statistics, Census 2021 topic summary tables TS001, TS003, TS004, TS006, TS007A, TS021, TS032, TS037, TS044, TS045, TS052, TS054, TS062, TS066 and TS067 for Output Areas, bulk downloads from [Nomis](https://www.nomisweb.co.uk/sources/census_2021_bulk), downloaded on 27 October 2026. Source: Office for National Statistics licensed under the [Open Government Licence v3.0](https://www.nationalarchives.gov.uk/doc/open-government-licence/version/3/). - **Boundaries**: Office for National Statistics, _Output Areas (December 2021) Boundaries EW BGC (V2)_, from the Open Geography Portal. Source: Office for National Statistics licensed under the Open Government Licence v3.0. Contains OS data © Crown copyright and database right 2026. - **2021/2 Output Area Classification**: Geographic Data Service / UCL and the Office for National Statistics, UK OAC v3.2 (January 2025), from [data.geods.ac.uk](https://data.geods.ac.uk/dataset/73a16168-f2b0-4732-b36a-3ec0253e34fe), Open Government Licence v3.0. Reference: Wyszomierski, J., Longley, P. A., Singleton, A. D., Gale, C. and O’Brien, O. (2024), “A neighbourhood Output Area Classification from the 2021 and 2022 UK censuses”, _The Geographical Journal_ 190(2): e12550. - Further reading: Assunção, R. M., Neves, M. C., Câmara, G. and Freitas, C. da C. (2006), “Efficient regionalization techniques for socio-economic geographical units using minimum spanning trees”, _International Journal of Geographical Information Science_ 20(7): 797–811 (SKATER); Duque, J. C., Anselin, L. and Rey, S. J. (2012), “The max-p-regions problem”, _Journal of Regional Science_ 52(3): 397–419; Openshaw, S. (1984), _The Modifiable Areal Unit Problem_, CATMOG 38. - Code: [`geodemographic_classification.py`](https://www.technicalgis.com/code/geodemographic-classification-python/), tested with Python 3.12.0, scikit-learn 1.9.1, spopt 0.8.0, libpysal 4.15.1, networkx 3.7, geopandas 1.2.0, shapely 2.2.0, pyogrio 0.13.0, pyproj 3.8.0, numpy 2.5.4, pandas 3.1.0, scipy 1.18.1, matplotlib 3.11.2 and requests 2.34.2. ## Key terms - **Geodemographic classification**: A grouping of small areas, such as Output Areas or postcodes, into types based on many census and other variables, so that areas of the same type are similar wherever they are. - **Output Area (OA)**: The smallest geography for which Census 2021 results are published in England and Wales, usually 40–250 households and about 100–625 people. About five OAs make up a Lower layer Super Output Area (LSOA). - **Inverse hyperbolic sine (IHS)**: The transformation asinh(x) = ln(x + √(x² + 1)). For large x it behaves like a logarithm, which reduces skew, but unlike a logarithm it is defined at zero. The 2021 OAC uses it. - **Index score**: A cluster's value of a variable divided by the value for the whole study area, multiplied by 100. 100 means the same as the average; 200 means twice the average. - **Silhouette width**: For each area, (b − a) / max(a, b), where a is the mean distance to areas in its own cluster and b is the mean distance to areas in the nearest other cluster. It runs from −1 to 1; negative values mean an area is closer to another cluster. - **Adjusted Rand index (ARI)**: A measure of agreement between two partitions of the same items, corrected for chance: 1 for identical partitions and about 0 for unrelated ones. It does not depend on the cluster labels. - **Regionalisation**: Clustering with a spatial contiguity constraint, so that each group is a single connected region. SKATER, AZP and max-p are regionalisation methods. - **MAUP**: The modifiable areal unit problem: results of spatial analysis change when the same data are aggregated to different zones (the zoning effect) or to larger or smaller zones (the scale effect). ## Data for this tutorial - [Nomis](https://www.technicalgis.com/data-sources/united-kingdom/#src-nomis) · Census 2021 topic summary tables as bulk CSV downloads for every Output Area in England and Wales - [Open Geography Portal](https://www.technicalgis.com/data-sources/united-kingdom/#src-open-geography-portal) · ONS Output Areas (December 2021) generalised boundaries for Cardiff Licences and formats for each source are listed in the [GIS data sources catalogue](https://www.technicalgis.com/data-sources/). --- Source: https://www.technicalgis.com/tutorials/geodemographic-classification-python/