How to Handle Class Imbalance in a Spatial Classifier
On this page
Problem statement
Land cover classes are never evenly distributed. In the stretch of the Norfolk Broads used on this page, ESA WorldCover puts 2.7% of the land in herbaceous wetland and 0.42% in bare ground. A class covering 2.7% of the landscape gets 2.7% of a random sample, which is rarely enough to learn — and a classifier that ignores it can still post a high overall accuracy.
Spatial data adds a complication that generic imbalance advice misses: rare classes are clustered. A rare wetland class is not 2% of every neighbourhood; it is 100% of a few places and 0% everywhere else.
That means a random split can put an entire class in one fold, and resampling techniques that assume independent rows create near-duplicates that leak straight across folds.
Quick answer
Stratify the sample, weight the classes, and never report accuracy alone:
from sklearn.ensemble import RandomForestClassifier
from sklearn.metrics import balanced_accuracy_score, classification_report
from sklearn.model_selection import StratifiedGroupKFold, cross_val_predict
model = RandomForestClassifier(
n_estimators=200, min_samples_leaf=2,
class_weight="balanced_subsample", n_jobs=4, random_state=0)
sample = np.random.default_rng(0).choice(len(y_all), 2800, replace=False) # a plain random sample
y_true = y_all[sample]
y_pred = cross_val_predict(model, X_all[sample], y_true, groups=g_all[sample],
cv=StratifiedGroupKFold(5, shuffle=True, random_state=0))
print(classification_report(y_true, y_pred, digits=3, target_names=NAMES.values()))
print(f"balanced accuracy {balanced_accuracy_score(y_true, y_pred):.3f}") precision recall f1-score support
tree cover 0.685 0.601 0.640 243
grassland 0.810 0.793 0.801 686
cropland 0.806 0.861 0.833 837
built-up 0.653 0.756 0.701 82
bare 0.538 0.636 0.583 11
water 0.994 0.985 0.989 857
herbaceous wetland 0.484 0.369 0.419 84
accuracy 0.841 2800
macro avg 0.710 0.715 0.710 2800
weighted avg 0.839 0.841 0.839 2800
balanced accuracy 0.715
Balanced accuracy is the mean per-class recall, so a class the model never predicts contributes zero rather than being hidden. Here the weighted forest is right on 84% of pixels overall and misses most of the wetland; the balanced figure says so, the overall one does not.
Stratify, keep groups, weight, report per class
1. Stratify the training sample
per_class = 300
for label in classes:
idx = np.nonzero(labels == label)[0]
picked.append(rng.choice(idx, min(per_class, len(idx)), replace=False))Equal counts per class give the model enough examples of each to learn a boundary. It also changes the class prior, so predicted probabilities no longer reflect landscape frequencies — a trade to make deliberately.
Set a minimum: a class with fewer than about 30 samples cannot be learned reliably, and including it produces a class the model predicts almost at random.
2. Keep spatial groups intact
A rare class often comes from a handful of polygons. Random splitting puts samples from one polygon in both folds, and the model memorises the polygon rather than the class.
Group by polygon, plot or survey unit, and use GroupKFold or StratifiedGroupKFold. Where a class exists in only two polygons, no split can validate it honestly — say so rather than reporting a number.
3. Prefer class weights to resampling
class_weight="balanced" reweights the loss so each class contributes equally, without duplicating rows.
Oversampling duplicates rows, which under any random split puts identical rows in different folds. Synthetic oversampling such as SMOTE interpolates between neighbours in feature space, creating rows that are spatially meaningless and leak in the same way. What duplication does to the score is easy to measure:
from sklearn.metrics import recall_score
from sklearn.model_selection import StratifiedKFold
rng = np.random.default_rng(0)
base = rng.choice(len(y_all), 2800, replace=False)
rare = [c for c in (60, 90)]
dup = np.concatenate([base] + [rng.choice(base[y_all[base] == c], 400) for c in rare]) # copies
for label, rows in (("as sampled", base), ("rare classes duplicated to 400", dup)):
for cv_name, cv, groups in (("random folds", StratifiedKFold(5, shuffle=True, random_state=0), None),
("patch-grouped folds", StratifiedGroupKFold(5, shuffle=True, random_state=0), g_all[rows])):
pred = cross_val_predict(RandomForestClassifier(200, min_samples_leaf=2, n_jobs=4, random_state=0),
X_all[rows], y_all[rows], cv=cv, groups=groups)
r = recall_score(y_all[rows], pred, labels=rare, average=None)
print(f"{label:31} {cv_name:20} bare {r[0]:.2f} wetland {r[1]:.2f}")as sampled random folds bare 0.82 wetland 0.31
as sampled patch-grouped folds bare 0.55 wetland 0.26
rare classes duplicated to 400 random folds bare 1.00 wetland 0.99
rare classes duplicated to 400 patch-grouped folds bare 0.63 wetland 0.49
With random folds, duplicating the rare rows makes their recall look excellent: most test rows have an identical twin in training. Grouped by patch, the copies stay together and the score falls back towards what the model can actually do. The duplication taught it nothing; it only taught the validation to agree.
For spatial data, weighting is almost always the better tool.
4. Report per-class metrics
Recall is what matters for a rare class. In the quick answer's report, the wetland row is the one to read: most wetland pixels were missed, and that is invisible in an overall accuracy in the 0.8s.
5. Set the threshold from the cost, not from the default
For a two-class problem, the 0.5 default is arbitrary. If missing the rare class is expensive and a false positive is cheap, lower it — and report the threshold with the results.
Grouped sampling and class weights on the Norfolk Broads
The data used on this page
A 17 × 20 km window of the Norfolk Broads — Hickling and Horsey broads, the dunes between Sea Palling and Winterton, and the coast down to Great Yarmouth — on a 20 m grid. The labels are ESA WorldCover 2021 (CC BY 4.0), resampled by majority; the features are a median of three clear Sentinel-2 L2A passes from June to August 2021 through Element 84's Earth Search, masked with each scene's classification band: six bands plus NDVI and NDWI. WorldCover is itself a model's output, not ground truth, so the numbers below describe how the methods behave, not the true accuracy of a Broads land cover map. The spatial groups are connected patches of each class.
import numpy as np
import pystac_client
import rasterio
from rasterio.warp import Resampling, reproject, transform_bounds
from rasterio.windows import bounds as window_bounds, from_bounds
from scipy import ndimage
BBOX = (1.53, 52.60, 1.78, 52.78)
NAMES = {10: "tree cover", 30: "grassland", 40: "cropland", 50: "built-up",
60: "bare", 80: "water", 90: "herbaceous wetland"}
IDS = ["S2A_31UDU_20210531_0_L2A", "S2A_31UDU_20210717_0_L2A", "S2B_31UDU_20210811_0_L2A"]
BANDS = ["blue", "green", "red", "nir", "swir16", "swir22"]
items = {i.id: i for i in pystac_client.Client.open("https://earth-search.aws.element84.com/v1").search(
collections=["sentinel-2-l2a"], ids=IDS).items()}
passes = []
for item_id in IDS:
item = items[item_id]
with rasterio.open(item.assets["swir16"].href) as src: # the 20 m grid
win = from_bounds(*transform_bounds(4326, src.crs, *BBOX), src.transform).round_offsets().round_lengths()
shape, grid, crs, box = (int(win.height), int(win.width)), src.window_transform(win), src.crs, window_bounds(win, src.transform)
layers = []
for band in BANDS:
with rasterio.open(item.assets[band].href) as src:
layers.append(src.read(1, window=from_bounds(*box, src.transform), out_shape=shape,
resampling=Resampling.average).astype("float32"))
with rasterio.open(item.assets["scl"].href) as src:
clear = np.isin(src.read(1, window=win), [4, 5, 6, 7, 11])
stack = np.stack(layers); stack[:, ~clear] = np.nan
passes.append(stack)
with np.errstate(all="ignore"):
comp = np.nanmedian(np.stack(passes), axis=0)
ndvi = (comp[3] - comp[2]) / (comp[3] + comp[2])
ndwi = (comp[1] - comp[3]) / (comp[1] + comp[3])
labels = np.zeros(shape, "uint8")
with rasterio.open("https://esa-worldcover.s3.eu-central-1.amazonaws.com/v200/2021/map/"
"ESA_WorldCover_10m_2021_v200_N51E000_Map.tif") as src:
reproject(rasterio.band(src, 1), labels, dst_transform=grid, dst_crs=crs, resampling=Resampling.mode)
patches, offset = np.zeros(shape, "int64"), 0
for c in NAMES:
lab, n = ndimage.label(labels == c, structure=np.ones((3, 3)))
patches[lab > 0] = lab[lab > 0] + offset
offset += n
edge = np.zeros(shape, bool) # a pixel with a different class next to it
for c in NAMES:
inside = labels == c
edge |= inside & ~ndimage.binary_erosion(inside, structure=np.ones((3, 3)))
X_all = np.concatenate([comp, ndvi[None], ndwi[None]]).reshape(8, -1).T
y_all, g_all, e_all = labels.ravel(), patches.ravel(), edge.ravel()
ok = np.isfinite(X_all).all(axis=1) & np.isin(y_all, list(NAMES))
X_all, y_all, g_all, e_all = X_all[ok], y_all[ok], g_all[ok], e_all[ok]
for c, name in NAMES.items():
print(f"{name:19} {(y_all == c).mean():7.2%} of {len(y_all):,} pixels, in {len(np.unique(g_all[y_all == c])):5,} patches")tree cover 7.59% of 876,400 pixels, in 1,595 patches
grassland 24.85% of 876,400 pixels, in 2,349 patches
cropland 28.95% of 876,400 pixels, in 247 patches
built-up 3.01% of 876,400 pixels, in 977 patches
bare 0.42% of 876,400 pixels, in 33 patches
water 32.51% of 876,400 pixels, in 128 patches
herbaceous wetland 2.66% of 876,400 pixels, in 552 patches
Example 1 — stratified, grouped sampling
import numpy as np
def stratified_grouped_sample(labels, groups, per_class=300, min_class=30,
min_groups=3, seed=0):
"""Equal-ish rows per class, drawn from as many groups as possible."""
rng = np.random.default_rng(seed)
picked, notes = [], []
for label in np.unique(labels):
idx = np.nonzero(labels == label)[0]
available_groups = np.unique(groups[idx])
if len(idx) < min_class:
notes.append(f"class {label}: only {len(idx)} samples — dropped")
continue
if len(available_groups) < min_groups:
notes.append(f"class {label}: only {len(available_groups)} "
"spatial groups — cannot be validated honestly")
per_group = -(-per_class // len(available_groups)) # ceiling, then trim below
chosen = []
for group in available_groups:
in_group = idx[groups[idx] == group]
chosen.append(rng.choice(in_group,
min(per_group, len(in_group)),
replace=False))
chosen = np.concatenate(chosen)
if len(chosen) > per_class:
chosen = rng.choice(chosen, per_class, replace=False)
picked.append(chosen)
print(f" class {label}: {len(idx):7,} available in "
f"{len(available_groups):3d} groups -> {len(chosen):5,} taken")
for note in notes:
print(f" ! {note}")
return np.concatenate(picked)strat = stratified_grouped_sample(y_all, g_all, per_class=400) class 10: 66,562 available in 1595 groups -> 400 taken
class 30: 217,820 available in 2349 groups -> 400 taken
class 40: 253,699 available in 247 groups -> 400 taken
class 50: 26,409 available in 977 groups -> 400 taken
class 60: 3,664 available in 33 groups -> 124 taken
class 80: 284,934 available in 128 groups -> 366 taken
class 90: 23,312 available in 552 groups -> 400 taken
The ceiling in per_group is a fix. The first version used per_class // len(groups), which for cropland — 247 patches — came to one pixel per patch: 247 pixels instead of 400, and a cropland recall of 0.34 in Example 2. Rounding up and trimming afterwards gives every class its full quota whenever the pixels exist.
Spreading the sample across groups is what stops a class being represented by one polygon. A class drawn entirely from one wetland teaches the model that wetland looks like that wetland.
Example 2 — weighting rather than resampling
import numpy as np
from sklearn.ensemble import RandomForestClassifier
from sklearn.metrics import balanced_accuracy_score, classification_report
from sklearn.model_selection import StratifiedGroupKFold
def fit_balanced(X, y, groups, class_names=None, n_splits=5, seed=0):
"""Class-weighted forest, validated with groups kept together."""
classes, counts = np.unique(y, return_counts=True)
print(f" {len(classes)} classes, "
f"imbalance ratio {counts.max() / counts.min():.1f}:1")
for label, count in zip(classes, counts):
name = class_names.get(label, label) if class_names else label
print(f" {name:18} {count:7,} ({count / len(y):6.2%})")
predicted = np.empty(len(y), dtype=y.dtype)
splitter = StratifiedGroupKFold(n_splits=n_splits, shuffle=True,
random_state=seed)
for train, test in splitter.split(X, y, groups):
model = RandomForestClassifier(
n_estimators=200, min_samples_leaf=2,
class_weight="balanced_subsample", n_jobs=-1, random_state=seed)
predicted[test] = model.fit(X[train], y[train]).predict(X[test])
print(classification_report(y, predicted, digits=3,
target_names=[str(class_names.get(c, c))
for c in classes]
if class_names else None))
print(f" balanced accuracy {balanced_accuracy_score(y, predicted):.3f}")
print(f" overall accuracy {(predicted == y).mean():.3f}")
return predictedprint("A plain random sample of 2,800 pixels")
random_pred = fit_balanced(X_all[sample], y_all[sample], g_all[sample], NAMES)
print("\nThe stratified, grouped sample from Example 1")
strat_pred = fit_balanced(X_all[strat], y_all[strat], g_all[strat], NAMES)
for c in (30, 40):
r, st = sample[y_all[sample] == c], strat[y_all[strat] == c]
print(f"{NAMES[c]:10} pixels on a patch edge: random {e_all[r].mean():.0%}, stratified-grouped {e_all[st].mean():.0%}")A plain random sample of 2,800 pixels
7 classes, imbalance ratio 77.9:1
tree cover 243 ( 8.68%)
grassland 686 (24.50%)
cropland 837 (29.89%)
built-up 82 ( 2.93%)
bare 11 ( 0.39%)
water 857 (30.61%)
herbaceous wetland 84 ( 3.00%)
precision recall f1-score support
tree cover 0.685 0.601 0.640 243
grassland 0.810 0.793 0.801 686
cropland 0.806 0.861 0.833 837
built-up 0.653 0.756 0.701 82
bare 0.538 0.636 0.583 11
water 0.994 0.985 0.989 857
herbaceous wetland 0.484 0.369 0.419 84
accuracy 0.841 2800
macro avg 0.710 0.715 0.710 2800
weighted avg 0.839 0.841 0.839 2800
balanced accuracy 0.715
overall accuracy 0.841
The stratified, grouped sample from Example 1
7 classes, imbalance ratio 3.2:1
tree cover 400 (16.06%)
grassland 400 (16.06%)
cropland 400 (16.06%)
built-up 400 (16.06%)
bare 124 ( 4.98%)
water 366 (14.70%)
herbaceous wetland 400 (16.06%)
precision recall f1-score support
tree cover 0.563 0.545 0.554 400
grassland 0.510 0.532 0.521 400
cropland 0.599 0.438 0.506 400
built-up 0.635 0.743 0.684 400
bare 0.828 0.895 0.860 124
water 0.848 0.899 0.873 366
herbaceous wetland 0.635 0.640 0.638 400
accuracy 0.642 2490
macro avg 0.660 0.670 0.662 2490
weighted avg 0.639 0.642 0.637 2490
balanced accuracy 0.670
overall accuracy 0.642
grassland pixels on a patch edge: random 35%, stratified-grouped 95%
cropland pixels on a patch edge: random 21%, stratified-grouped 71%
The two overall accuracies are not comparable — the stratified sample's test rows are balanced by construction — but the per-class recalls are. The rare classes gain: wetland recall rises from 0.37 to 0.64 and bare ground from 0.64 to 0.90. The common classes lose, and the reason is the grouping, not the balance: one pixel from each of 2,349 grassland patches is mostly the edge of a small patch — 95% of them sit next to a different class, against 35% in the random sample — and an edge pixel is a mixed pixel, hard for anyone. Spreading a sample across patches makes it harder and more honest at the same time; it is the rare-class recall it is there to protect.
Printing both accuracies side by side makes the imbalance visible in one line. A large gap between them means the model is doing well on the common classes and poorly on the rare ones.
class_weight="balanced_subsample" recomputes the weights for each tree's bootstrap sample, which suits a forest better than a single global weighting.
Example 3 — a confusion matrix that shows where the errors go
import numpy as np
from sklearn.metrics import confusion_matrix
def confusion_report(y_true, y_pred, class_names=None, normalise="true"):
"""Row-normalised confusion, so recall is readable directly."""
labels = np.unique(np.concatenate([y_true, y_pred]))
matrix = confusion_matrix(y_true, y_pred, labels=labels,
normalize=normalise)
names = [str(class_names.get(l, l)) if class_names else str(l)
for l in labels]
width = max(len(n) for n in names) + 1
print(" " * (width + 2) + " ".join(f"{n[:6]:>7}" for n in names))
for i, name in enumerate(names):
row = " ".join(f"{matrix[i, j]:7.2f}" for j in range(len(labels)))
print(f" {name:{width}} {row}")
for i, name in enumerate(names):
confused_with = np.argsort(-matrix[i])
worst = [j for j in confused_with if j != i][0]
if matrix[i, worst] > 0.15:
print(f" ! {name} is confused with {names[worst]} "
f"{matrix[i, worst]:.0%} of the time")
return matrix_ = confusion_report(y_all[strat], strat_pred, NAMES) tree c grassl cropla built- bare water herbac
tree cover 0.55 0.18 0.04 0.10 0.00 0.01 0.13
grassland 0.16 0.53 0.13 0.12 0.00 0.00 0.06
cropland 0.08 0.19 0.44 0.16 0.01 0.01 0.12
built-up 0.05 0.08 0.06 0.74 0.04 0.03 0.01
bare 0.00 0.00 0.02 0.05 0.90 0.04 0.00
water 0.01 0.00 0.00 0.03 0.01 0.90 0.05
herbaceous wetland 0.12 0.07 0.06 0.01 0.00 0.09 0.64
! tree cover is confused with grassland 18% of the time
! grassland is confused with tree cover 16% of the time
! cropland is confused with grassland 19% of the time
Row normalisation makes each row's diagonal the recall for that class, and the off-diagonals show exactly which class the errors go to. That is far more actionable than a single score: "wetland is being called water 40% of the time" points at a feature that could separate them.
Rare classes come in patches
Why rare classes are clustered, and why that matters
A rare land cover class is not scattered uniformly at 2% density. It occupies a few patches where it is the dominant class.
That has two consequences generic imbalance advice does not address.
A random split puts samples from one patch in both training and test folds, so the model memorises the patch and the validation rewards it. The reported recall for the rare class is then optimistic by a large margin.
And a class present in only two or three patches cannot be validated at all: any split either trains on all of them or tests on a patch the model has never seen anything like.
Why weighting beats resampling here
Oversampling duplicates rows. Under GroupKFold the duplicates stay together, but under any random element they can separate, and a duplicated row in both folds is perfect leakage.
Synthetic methods such as SMOTE interpolate between neighbours in feature space to create new rows. Those rows have no location, so they cannot be assigned to a spatial group, and they are typically dropped from spatial splitting or assigned arbitrarily.
Class weighting changes the loss without changing the rows, so the spatial structure — and therefore the validation — stays intact.
Why accuracy is the wrong headline
With one class at 90% of the landscape, a model that always predicts it scores 0.90. That is a strong-sounding number for a model that has learned nothing.
Balanced accuracy — the mean recall across classes — gives that model 1/n, where n is the number of classes. For five classes that is 0.20, which is a fairer description.
Per-class recall is more informative still, because it says which classes are being missed, and that usually points at a specific confusion worth fixing.
Why stratification changes what the probabilities mean
A classifier learns the class frequencies of its training data. Training on 300 samples of each of five classes teaches it that all five are equally likely.
The decision boundaries usually improve, because the rare classes get enough examples. The predicted probabilities become conditional on the sampling design rather than on the landscape.
Where the probabilities matter — a suitability surface, a risk map — correct them back to the landscape priors, or weight rather than resample so the priors are preserved.
Group counts, balanced accuracy and normalised confusion
- Rare classes are clustered. Group by patch and use
StratifiedGroupKFold. - A class in fewer than about three groups cannot be validated honestly.
- Prefer
class_weightto resampling — duplicates and synthetic rows leak. - Report balanced accuracy and per-class recall, never accuracy alone.
- Stratifying changes the prior, so probabilities need correcting.
- Set a minimum class size, around 30 samples.
- Row-normalise the confusion matrix so the diagonal is recall.
- Report the decision threshold if you moved it from the default.
Internal links
- How to sample training points from rasters and polygons — stratified, grouped sampling
- How to classify land cover from satellite imagery in Python — the classification workflow
- How to do spatial cross-validation in Python — keeping groups intact
- How to evaluate a spatial model honestly — reporting beyond one number
- Spatial machine learning explained — why rows are not independent
- Spatial leakage explained: why random cross-validation lies — what duplicated rows do
- My model scores 0.95 in testing and fails in the field — the same problem in regression
- How to build spatial features for a machine learning model — features that separate confused classes
FAQ
How do I handle class imbalance in a land cover classifier?
Stratify the sample so each class has enough rows, use class_weight="balanced" rather than resampling, keep spatial groups intact when splitting, and report per-class recall.
Why not use SMOTE or oversampling?
Both create rows that leak across spatial folds — duplicates directly, synthetic rows because they have no location and cannot be grouped.
Why is my accuracy high but the map wrong?
Overall accuracy is dominated by the majority class. A model ignoring a 2% class still scores 98%. Report balanced accuracy and per-class recall.
How many samples does a rare class need?
At least about 30, and from at least three separate spatial groups. Fewer than that and it cannot be validated honestly.
Does stratified sampling bias my model?
It changes the class prior, so predicted probabilities no longer reflect landscape frequencies. Decision boundaries usually improve. Correct the priors if the probabilities matter.
What metric should I report?
Per-class precision and recall, plus balanced accuracy. A row-normalised confusion matrix shows which classes are being confused.
Should I move the decision threshold?
If the costs of the two error types differ, yes — and report the threshold with the results.