How to Classify Land Cover from Satellite Imagery in Python
Problem statement
Supervised land cover classification is four steps: sample training pixels, build features, fit a classifier, predict the scene. Each has a decision that matters more than the algorithm choice.
Measured on a real Sentinel-2 scene over Snowdonia, the class medians that a classifier has to separate:
blue green red nir
vegetation 0.039 0.067 0.056 0.361
water 0.015 0.013 0.008 0.006
bare soil 0.090 0.106 0.110 0.235
cloud 0.527 0.507 0.504 0.599
cloud shadow 0.030 0.039 0.032 0.134
Water and cloud shadow are both dark in the visible bands and separate cleanly in the near-infrared — 0.006 against 0.134, a factor of twenty. That is the kind of contrast a classifier exploits, and it only exists if you include the right bands.
Quick answer
import numpy as np
from sklearn.ensemble import RandomForestClassifier
from sklearn.metrics import balanced_accuracy_score, classification_report
from sklearn.model_selection import StratifiedGroupKFold
model = RandomForestClassifier(n_estimators=200, min_samples_leaf=2,
class_weight="balanced_subsample",
n_jobs=-1, random_state=0)
predicted = np.empty(len(y), dtype=y.dtype)
for train, test in StratifiedGroupKFold(5, shuffle=True,
random_state=0).split(X, y, groups):
predicted[test] = model.fit(X[train], y[train]).predict(X[test])
print(classification_report(y, predicted, digits=3))
print(f"balanced accuracy {balanced_accuracy_score(y, predicted):.3f}")The groups argument is what makes the score meaningful. Training polygons are clustered, and without grouping the model memorises polygons rather than classes.
Step-by-step solution
1. Mask before you sample
Cloud, cirrus and shadow are not land cover. Measured on a single Sentinel-2 date over Snowdonia, 51% of pixels were cloud, cirrus or cloud shadow — a classifier trained on those learns to classify weather.
usable = ~np.isin(scl, (0, 1, 3, 8, 9, 10))Alternatively, include cloud and shadow as explicit classes so the model can label them. Either is defensible; leaving them unlabelled in the training data is not.
2. Build features beyond the raw bands
Raw reflectance separates the obvious classes. Indices separate the hard ones:
ndvi = (nir - red) / (nir + red) # vegetation vs bare
ndwi = (green - nir) / (green + nir) # water, 3.79 sd from shadow
ndbi = (swir - nir) / (swir + nir) # built-up, but reads high over waterMeasured separation between water and cloud shadow: NDWI 3.79 pooled standard deviations, NDVI 3.12, MNDWI 2.60. Adding NDWI is worth more than any hyperparameter.
Texture features — a local standard deviation over a window — separate classes with similar spectra and different structure, such as forest and grassland.
3. Sample training pixels from polygons, and keep the polygon id
samples = sample_in_polygons(training_polygons, per_polygon=50,
label_column="class")Fifty pixels from one polygon are fifty views of one place. They must stay in the same cross-validation fold, which needs the polygon id carried through as a group.
Spread the sample across as many polygons per class as possible. A class drawn from one polygon teaches the model that class looks like that polygon.
4. Validate with grouped, stratified folds
StratifiedGroupKFold keeps each polygon whole and each fold's class balance similar. Without grouping, the reported accuracy is memorisation; without stratification, a rare class can be absent from a fold entirely.
5. Report per-class recall, not accuracy
A class covering 2% of the landscape can be ignored entirely while accuracy stays at 98%. Balanced accuracy — the mean per-class recall — does not hide it.
Code examples
Example 1 — building the feature stack
import numpy as np
def classification_features(bands, include_indices=True, texture_window=5):
"""Reflectance, indices and texture, as an aligned stack."""
features, names = [], []
for name in ("blue", "green", "red", "nir", "swir16", "swir22"):
if name in bands:
features.append(bands[name])
names.append(name)
if include_indices:
def nd(a, b):
total = a + b
with np.errstate(divide="ignore", invalid="ignore"):
return np.where(np.abs(total) < 1e-6, np.nan, (a - b) / total)
pairs = [("ndvi", "nir", "red"), ("ndwi", "green", "nir")]
if "swir16" in bands:
pairs += [("mndwi", "green", "swir16"), ("ndbi", "swir16", "nir")]
for label, first, second in pairs:
if first in bands and second in bands:
features.append(nd(bands[first], bands[second]))
names.append(label)
if texture_window:
from scipy.ndimage import uniform_filter
base = np.nan_to_num(bands["nir"], nan=float(np.nanmean(bands["nir"])))
mean = uniform_filter(base, size=texture_window)
squared = uniform_filter(base ** 2, size=texture_window)
features.append(np.sqrt(np.maximum(squared - mean ** 2, 0)))
names.append(f"nir_texture_{texture_window}")
stack = np.stack(features)
print(f" {len(names)} features: {names}")
print(f" {np.isfinite(stack).all(axis=0).mean():.1%} of cells complete")
return stack, namesTexture is the feature most often missing from a land cover model. Forest and grassland can have very similar mean reflectance and very different local variance, and no amount of spectral tuning separates them.
Example 2 — sampling from training polygons with groups
import numpy as np
import geopandas as gpd
import rasterio
def sample_training(polygons, stack, transform, label_column="class",
per_polygon=50, min_polygons=3, seed=0):
"""Pixels inside training polygons, with the polygon id as the group."""
from rasterio.features import geometry_mask
rng = np.random.default_rng(seed)
height, width = stack.shape[1:]
rows, labels, groups = [], [], []
per_class_polygons = polygons.groupby(label_column).size()
for idx, polygon in polygons.iterrows():
mask = ~geometry_mask([polygon.geometry], out_shape=(height, width),
transform=transform, invert=False)
valid = mask & np.isfinite(stack).all(axis=0)
candidates = np.nonzero(valid.ravel())[0]
if candidates.size == 0:
continue
take = rng.choice(candidates, min(per_polygon, candidates.size),
replace=False)
r, c = np.unravel_index(take, (height, width))
rows.append(stack[:, r, c].T)
labels.append(np.full(len(take), polygon[label_column]))
groups.append(np.full(len(take), idx))
X = np.vstack(rows)
y = np.concatenate(labels)
g = np.concatenate(groups)
print(f" {len(X):,} pixels from {len(np.unique(g))} polygons")
for label, count in per_class_polygons.items():
n_pixels = int((y == label).sum())
flag = " <- too few polygons to validate" if count < min_polygons else ""
print(f" {str(label):18} {count:3d} polygons, "
f"{n_pixels:6,} pixels{flag}")
return X, y, gThe polygon count per class is the number that decides whether the class can be validated. Two polygons means any split either trains on both or tests on ground the model has seen nothing like.
Example 3 — predicting the scene with a confidence band
import numpy as np
import rasterio
def classify_scene(model, stack, profile, out_path, class_names=None):
"""Class raster plus the winning probability, as two bands."""
height, width = stack.shape[1:]
valid = np.isfinite(stack).all(axis=0)
X = stack[:, valid].T
classes = model.predict(X)
probabilities = model.predict_proba(X).max(axis=1)
class_raster = np.full(height * width, 0, dtype="uint8")
confidence = np.full(height * width, np.nan, dtype="float32")
class_raster[valid.ravel()] = classes.astype("uint8")
confidence[valid.ravel()] = probabilities
profile = dict(profile) | {"count": 2, "dtype": "float32",
"nodata": np.nan, "compress": "deflate"}
with rasterio.open(out_path, "w", **profile) as dst:
dst.write(class_raster.reshape(height, width).astype("float32"), 1)
dst.write(confidence.reshape(height, width), 2)
dst.set_band_description(1, "class")
dst.set_band_description(2, "winning class probability")
if class_names:
dst.update_tags(classes=",".join(f"{k}={v}"
for k, v in class_names.items()))
low = confidence[valid.ravel()] < 0.5
print(f" {valid.mean():.1%} of the scene classified")
print(f" {low.mean():.1%} of classified pixels have a winning "
"probability below 0.5")
return out_pathThe confidence band is the land cover equivalent of an applicability mask. Pixels where the winning probability is barely above the next one are the model guessing between two classes, and they cluster along class boundaries and in unfamiliar terrain.
Explanation
Why the near-infrared does the work
The class table at the top shows why. Water and cloud shadow differ by a factor of five in the visible bands and a factor of twenty in the near-infrared. Vegetation differs from bare soil by a factor of 1.5 in red and 6 in the ratio.
A classifier given only visible bands has to separate classes that genuinely overlap there. Adding the near-infrared and shortwave-infrared bands turns several hard separations into easy ones.
This is why band selection matters more than model selection: the model can only use contrasts that exist in its inputs.
Why grouping by polygon is not optional
Training polygons are drawn around homogeneous patches, so the pixels inside one are nearly identical.
Under a random split, a held-out pixel has fifty near-duplicates in the training set. The model predicts it perfectly, and the reported accuracy measures memorisation of that patch.
Grouping by polygon forces the model to generalise from one patch to another, which is what classifying a scene actually requires. The score always drops, and the drop is the size of the illusion.
Why to include or exclude clouds deliberately
Cloud and shadow are 51% of a single Sentinel-2 date over Snowdonia. They must be handled explicitly.
Masking them out means the classifier never sees them and will assign them a land cover class at prediction time — usually bare soil for cloud and water for shadow, since those are the nearest spectral matches.
Including them as classes means the model can label them, at the cost of two classes that are not land cover. For a single-date classification that is usually the better choice; for a composite, masking is right because the composite should not contain cloud at all.
Why per-pixel classification produces speckle
Each pixel is classified independently, so a boundary pixel that is a mixture of two classes is assigned to whichever it resembles more. The result is a scatter of isolated pixels along every boundary.
Options: include texture features so the model sees context, apply a majority filter afterwards, or segment the image and classify segments rather than pixels.
A majority filter is the cheapest and it erodes genuine small features. Segmentation is more principled and adds a parameter — the segment size — that becomes part of the analysis.
Edge cases or notes
- Mask cloud and shadow before sampling, or include them as explicit classes.
- Include the near-infrared and shortwave-infrared bands. Several separations exist only there.
- Add indices; NDWI separated water from shadow by 3.79 standard deviations.
- Group by training polygon at the split, always.
- Three or more polygons per class, or the class cannot be validated.
- Report per-class recall and balanced accuracy, not overall accuracy.
- Write a confidence band alongside the class raster.
- Expect speckle from per-pixel classification; texture features or segmentation reduce it.
Internal links
- Spectral bands explained: what satellite imagery actually measures — the class signatures
- Spectral indices explained — the features that separate hard classes
- Cloud masking explained — handling the 51%
- How to sample training points from rasters and polygons — the grouped sampling
- How to handle class imbalance in a spatial classifier — rare classes
- How to do spatial cross-validation in Python — grouped folds
- How to turn model predictions back into a raster — writing the class map
- How to load Sentinel-2 bands into Python as an analysis-ready array — the input stack
FAQ
How do I classify land cover in Python?
Mask cloud, build a feature stack of bands plus indices and texture, sample pixels from training polygons keeping the polygon id, fit a classifier with grouped stratified folds, and predict the scene.
Which bands do I need?
At least the visible and near-infrared. Shortwave-infrared separates several classes that overlap elsewhere — water and shadow differ twentyfold in the near-infrared alone.
Why must I group by training polygon?
Because pixels inside one polygon are near-duplicates. Without grouping, the score measures memorisation of the patch rather than knowledge of the class.
Should clouds be a class or be masked?
Either, deliberately. Masking means the model assigns cloud a land cover class at prediction time; including it as a class means it can be labelled.
Why is my classified map speckled?
Per-pixel classification treats each pixel independently, so mixed boundary pixels scatter. Add texture features, apply a majority filter, or classify segments instead.
What accuracy should I report?
Per-class precision and recall plus balanced accuracy. Overall accuracy is dominated by the largest class.
How many training polygons per class?
At least three, and preferably many more. A class from one polygon teaches the model that polygon rather than the class.