How to Classify Land Cover from Satellite Imagery in Python
On this page
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 Snowdon — 22 September 2025, a 10 km window round the summit, classes taken from the scene's own classification layer — the median reflectances a classifier has to separate:
blue green red nir
vegetation 0.039 0.068 0.057 0.362
water 0.015 0.013 0.008 0.007
bare soil 0.106 0.123 0.126 0.285
cloud 0.336 0.330 0.327 0.454
cloud shadow 0.029 0.037 0.031 0.133
Water and cloud shadow are both dark in the visible bands and separate cleanly in the near-infrared — 0.007 against 0.133, a factor of nineteen. That is the kind of contrast a classifier exploits, and it only exists if you include the right bands. (The setup below prints this table.)
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. Example 2 runs this on the Snowdon scene, with and without the groups.
What you get — from 2. Build features beyond the raw bands
NDWI 3.02 sd
NDVI 2.76 sd
MNDWI 1.92 sd
NIR alone 1.43 sd
Mask first, richer features, polygon-based sampling
1. Mask before you sample
Cloud, cirrus and shadow are not land cover. On the single Sentinel-2 date over Snowdon used here, 57% 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 vs shadow
ndbi = (swir - nir) / (swir + nir) # built-up, but reads high over waterHow far apart water and cloud shadow sit on each, in pooled standard deviations:
def nd(a, b):
with np.errstate(divide="ignore", invalid="ignore"):
return (a - b) / (a + b)
water, shadow = scl == 6, scl == 3
for name, index in (("NDWI", nd(bands["green"], bands["nir"])), ("NDVI", nd(bands["nir"], bands["red"])),
("MNDWI", nd(bands["green"], bands["swir16"])), ("NIR alone", bands["nir"])):
a, b = index[water], index[shadow]
pooled = np.sqrt((np.nanvar(a) + np.nanvar(b)) / 2)
print(f"{name:9} {abs(np.nanmean(a) - np.nanmean(b)) / pooled:5.2f} sd")NDWI 3.02 sd
NDVI 2.76 sd
MNDWI 1.92 sd
NIR alone 1.43 sd
Adding the right index 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.
Features, training samples and a classifier
The data used on this page
One Sentinel-2 L2A scene, S2A_30UVD_20250922_0_L2A (22 September 2025), from Element 84's Earth Search, cut to the 10 × 10 km window round the summit of Snowdon (Yr Wyddfa) that the NDVI and radiometric levels guides use: six bands at 10 m (the shortwave-infrared bilinearly resampled from 20 m) and the scene classification (SCL). The training labels combine two sources, and say so: land classes from ESA WorldCover 2021 where the SCL calls the pixel clear, and cloud and cloud shadow from the SCL itself. Up to 40 connected patches of at least 30 pixels per class become the training polygons. WorldCover and the SCL are models' outputs, so the scores below measure how the method behaves, not the true accuracy of a Snowdonia land cover map.
import geopandas as gpd
import numpy as np
import pystac_client
import rasterio
from rasterio.enums import Resampling
from rasterio.features import shapes
from rasterio.transform import from_origin
from rasterio.warp import reproject
from rasterio.windows import from_bounds
from scipy import ndimage
from shapely.geometry import shape
BOUNDS = (422_880, 5_875_440, 432_880, 5_885_440)
TRANSFORM = from_origin(422_880, 5_885_440, 10, 10)
item = next(pystac_client.Client.open("https://earth-search.aws.element84.com/v1")
.search(ids=["S2A_30UVD_20250922_0_L2A"], collections=["sentinel-2-l2a"]).items())
bands = {}
for name in ("blue", "green", "red", "nir", "swir16", "swir22"):
with rasterio.open(item.assets[name].href) as src:
bands[name] = src.read(1, window=from_bounds(*BOUNDS, transform=src.transform), out_shape=(1000, 1000),
resampling=Resampling.bilinear if name.startswith("swir") else Resampling.nearest
).astype("float32") / 10_000
with rasterio.open(item.assets["scl"].href) as src:
scl = src.read(1, window=from_bounds(*BOUNDS, transform=src.transform), out_shape=(1000, 1000),
resampling=Resampling.nearest)
print(f"cloud, cirrus or shadow: {np.isin(scl, (3, 8, 9, 10)).mean():.1%} of the window")
print(" blue green red nir")
for label, codes in (("vegetation", [4]), ("water", [6]), ("bare soil", [5]), ("cloud", [8, 9]), ("cloud shadow", [3])):
m = np.isin(scl, codes)
print(f"{label:14} " + " ".join(f"{np.median(bands[b][m]):7.3f}" for b in ("blue", "green", "red", "nir")))
worldcover = np.zeros((1000, 1000), "uint8")
with rasterio.open("https://esa-worldcover.s3.eu-central-1.amazonaws.com/v200/2021/map/"
"ESA_WorldCover_10m_2021_v200_N51W006_Map.tif") as src:
reproject(rasterio.band(src, 1), worldcover, dst_transform=TRANSFORM, dst_crs="EPSG:32630",
resampling=Resampling.mode)
CLASSES = {1: "tree cover", 2: "grassland", 3: "built-up", 4: "bare", 5: "water", 6: "cloud", 7: "cloud shadow"}
clear = np.isin(scl, (4, 5, 6, 7, 11))
label = np.zeros((1000, 1000), "uint8")
for code, wc in ((1, 10), (2, 30), (3, 50), (4, 60), (5, 80)):
label[clear & (worldcover == wc)] = code
label[np.isin(scl, (8, 9, 10))] = 6
label[scl == 3] = 7
rng, rows = np.random.default_rng(0), []
for code, name in CLASSES.items():
patches, n = ndimage.label(label == code)
sizes = ndimage.sum_labels(np.ones_like(patches), patches, index=np.arange(1, n + 1))
keep = np.nonzero(sizes >= 30)[0] + 1
keep = rng.choice(keep, min(40, len(keep)), replace=False)
for geom, _ in shapes(patches.astype("int32"), mask=np.isin(patches, keep), transform=TRANSFORM):
rows.append({"class": name, "geometry": shape(geom)})
training_polygons = gpd.GeoDataFrame(rows, crs="EPSG:32630")
print(training_polygons["class"].value_counts().to_dict())cloud, cirrus or shadow: 57.3% of the window
blue green red nir
vegetation 0.039 0.068 0.057 0.362
water 0.015 0.013 0.008 0.007
bare soil 0.106 0.123 0.126 0.285
cloud 0.336 0.330 0.327 0.454
cloud shadow 0.029 0.037 0.031 0.133
{'tree cover': 40, 'grassland': 40, 'cloud shadow': 40, 'cloud': 31, 'built-up': 16, 'water': 11, 'bare': 5}
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, namesstack, feature_names = classification_features(bands) 11 features: ['blue', 'green', 'red', 'nir', 'swir16', 'swir22', 'ndvi', 'ndwi', 'mndwi', 'ndbi', 'nir_texture_5']
100.0% of cells complete
Texture 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, gX, y, groups = sample_training(training_polygons, stack, TRANSFORM) 8,705 pixels from 183 polygons
bare 5 polygons, 233 pixels
built-up 16 polygons, 695 pixels
cloud 31 polygons, 1,464 pixels
cloud shadow 40 polygons, 2,000 pixels
grassland 40 polygons, 1,894 pixels
tree cover 40 polygons, 1,869 pixels
water 11 polygons, 550 pixels
The quick answer's model, run on these samples — and again with the groups taken away:
from sklearn.ensemble import RandomForestClassifier
from sklearn.metrics import balanced_accuracy_score, classification_report
from sklearn.model_selection import StratifiedGroupKFold, StratifiedKFold
model = RandomForestClassifier(n_estimators=200, min_samples_leaf=2,
class_weight="balanced_subsample", n_jobs=4, random_state=0)
for label, splitter, split_groups in (("grouped by polygon", StratifiedGroupKFold(5, shuffle=True, random_state=0), groups),
("random pixels", StratifiedKFold(5, shuffle=True, random_state=0), None)):
predicted = np.empty(len(y), dtype=y.dtype)
for train, test in splitter.split(X, y, split_groups):
predicted[test] = model.fit(X[train], y[train]).predict(X[test])
print(f"--- {label}: balanced accuracy {balanced_accuracy_score(y, predicted):.3f}")
if split_groups is not None:
print(classification_report(y, predicted, digits=3))--- grouped by polygon: balanced accuracy 0.484
precision recall f1-score support
bare 0.305 0.288 0.296 233
built-up 0.355 0.430 0.389 695
cloud 0.694 0.647 0.669 1464
cloud shadow 0.506 0.549 0.527 2000
grassland 0.360 0.287 0.319 1894
tree cover 0.542 0.601 0.570 1869
water 0.621 0.589 0.604 550
accuracy 0.506 8705
macro avg 0.483 0.484 0.482 8705
weighted avg 0.503 0.506 0.503 8705
--- random pixels: balanced accuracy 0.719
The 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_pathprofile = {"driver": "GTiff", "width": 1000, "height": 1000, "crs": "EPSG:32630", "transform": TRANSFORM}
codes = {str(name): i + 1 for i, name in enumerate(np.unique(y))}
numeric = RandomForestClassifier(n_estimators=200, min_samples_leaf=2, class_weight="balanced_subsample",
n_jobs=4, random_state=0).fit(X, np.array([codes[v] for v in y]))
classify_scene(numeric, stack, profile, "landcover.tif", {v: k for k, v in codes.items()})
with rasterio.open("landcover.tif") as src:
classes = src.read(1)
print({name: f"{(classes == code).mean():.1%}" for name, code in codes.items()}) 100.0% of the scene classified
29.0% of classified pixels have a winning probability below 0.5
{'bare': '0.4%', 'built-up': '1.9%', 'cloud': '32.1%', 'cloud shadow': '28.2%', 'grassland': '15.6%', 'tree cover': '19.9%', 'water': '2.0%'}
The 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.
Why near-infrared does the work
Why the near-infrared does the work
The class table at the top shows why. Water and cloud shadow differ by a factor of two to four in the visible bands and nineteen in the near-infrared. Vegetation differs from bare soil by a factor of 2.2 in red, and its near-infrared to red ratio is 2.8 times bare soil's.
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: on the Snowdon samples, balanced accuracy was 0.719 with random pixel folds and 0.484 with each polygon kept whole.
Why to include or exclude clouds deliberately
Cloud and shadow are 57% of this Sentinel-2 date over Snowdon. 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. What it chooses is easy to check — train on the land classes only, then predict the pixels the SCL calls cloud and shadow:
land = ~np.isin(y, ["cloud", "cloud shadow"])
land_only = RandomForestClassifier(n_estimators=200, min_samples_leaf=2, class_weight="balanced_subsample",
n_jobs=4, random_state=0).fit(X[land], y[land])
flat = stack.reshape(stack.shape[0], -1).T
for name, codes_ in (("cloud", (8, 9)), ("cloud shadow", (3,))):
guess = land_only.predict(flat[np.isin(scl.ravel(), codes_)])
kinds, counts = np.unique(guess, return_counts=True)
print(f"{name:13} ->", ", ".join(f"{k} {c / counts.sum():.0%}" for k, c in sorted(zip(kinds, counts), key=lambda t: -t[1])[:3]))cloud -> water 47%, grassland 44%, tree cover 4%
cloud shadow -> grassland 37%, water 33%, tree cover 26%
Not the nearest-looking class you might guess. Cloud, far brighter than anything in the training data, came out as water or grassland — a random forest cannot extrapolate beyond the reflectances it was trained on, so an unfamiliar pixel lands wherever the trees' last splits happen to send it. Shadow was spread across grassland, water and tree cover. Either way a masked-out class does not disappear; it turns up in the map as something else.
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.
Masking, extra bands and indices
- 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.02 standard deviations, against 1.43 for the near-infrared alone.
- 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
- How to classify crop types from a satellite time series — why one date is not enough for crops
- How to monitor crops with Sentinel-1 in Python — radar features for land cover
- How to fit a species distribution model in Python — the presence-background version of the same classifier
- Land cover products explained — WorldCover, Impact Observatory and Hansen compared
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 nineteenfold 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.