Classify land cover with a random forest in Python

A land cover map answers a basic question: what is on the ground? It tells you where there is forest, farmland, water and buildings. Governments use these maps to plan, farmers to monitor, and researchers to measure change. In this tutorial you will make one yourself, by teaching a random forest model to classify every pixel of a satellite image.

A random forest is a good first model. It handles many input bands, needs little tuning, and gives you a sense of which inputs matter. The important lessons in this tutorial are less about the algorithm and more about preparing data and testing honestly, because that is where most classification projects go wrong.

What you need

• Python 3.9 or newer, with rasterio, numpy and scikit-learn installed.
• stack.tif: a multiband image where all bands share the same grid. For example Sentinel-2 blue, green, red and near-infrared at 10 m, stacked into one file.
• labels.tif: a raster of the same size and grid where 0 means unlabelled, and 1, 2, 3 and 4 mark your classes, for example 1 water, 2 forest, 3 cropland, 4 built-up.

pip install rasterio numpy scikit-learn

You can create labels in QGIS by drawing polygons over areas you are sure about, giving each a class number, and converting them to a raster that matches your image. Aim for a few hundred well-spread pixels per class at the very least. Quality matters more than quantity: pick pure, unambiguous areas, and spread them across the whole region, not just one corner.

Step 1: turn pixels into rows

Machine learning models expect a table: one row per example, one column per feature. A satellite image is a stack of grids, so we reshape it. Every pixel becomes a row, and its values in each band become the columns.

import numpy as np
import rasterio

with rasterio.open("stack.tif") as src:
    image = src.read().astype("float32")   # shape: (bands, rows, cols)
    profile = src.profile

with rasterio.open("labels.tif") as lab:
    labels = lab.read(1)                    # shape: (rows, cols)

bands, rows, cols = image.shape
X_all = image.reshape(bands, -1).T          # one row per pixel, one column per band
y_all = labels.ravel()
labelled = y_all > 0
X, y = X_all[labelled], y_all[labelled]

print(X.shape, "labelled pixels;", np.bincount(y)[1:], "per class")

Check the count per class. If one class has ten times more examples than another, the model will favour the common one. Collect more labels for the rare class, or use class_weight="balanced" in the next step.

Step 2: train a first model

from sklearn.ensemble import RandomForestClassifier
from sklearn.model_selection import train_test_split
from sklearn.metrics import classification_report

X_train, X_test, y_train, y_test = train_test_split(
    X, y, test_size=0.3, stratify=y, random_state=42
)

model = RandomForestClassifier(
    n_estimators=200, n_jobs=-1, random_state=42, class_weight="balanced"
)
model.fit(X_train, y_train)

print(classification_report(y_test, model.predict(X_test)))

You will probably see very high scores, perhaps 95 percent or more. Do not celebrate yet.

Step 3: test honestly with spatial blocks

The fix is to keep whole areas together. We divide the image into square blocks, give each pixel the number of its block, and make sure that a block is either used for training or for testing, never both. Scikit-learn’s GroupKFold does exactly this.

from sklearn.model_selection import GroupKFold

row_idx, col_idx = np.indices((rows, cols))
block_id = (row_idx // 128) * 1000 + (col_idx // 128)   # 128 x 128 pixel blocks
groups = block_id.ravel()[labelled]

scores = []
for train_i, test_i in GroupKFold(n_splits=5).split(X, y, groups):
    m = RandomForestClassifier(n_estimators=200, n_jobs=-1, random_state=42,
                               class_weight="balanced")
    m.fit(X[train_i], y[train_i])
    scores.append(m.score(X[test_i], y[test_i]))

print(f"Spatial CV accuracy: {np.mean(scores):.1%} (+/- {np.std(scores):.1%})")

This number is usually lower than the first one, and it is the one to believe. If the gap is large, your model is relying on local memorisation, and you need more varied labels.

Step 4: look at where it makes mistakes

from sklearn.metrics import confusion_matrix

y_pred = model.predict(X_test)
print(confusion_matrix(y_test, y_pred))

Each row shows what the true class was, and each column what the model said. Large numbers off the diagonal show which classes are confused. Cropland and grass often mix, as do built-up land and bare soil. That tells you where to collect more labels or add extra information.

Step 5: map the whole image

Predicting every pixel of a large image in one go can exhaust your memory. We do it in strips of rows instead.

prediction = np.zeros((rows, cols), dtype="uint8")
step = 512

for r in range(0, rows, step):
    strip = image[:, r:r + step, :]                 # (bands, h, cols)
    h = strip.shape[1]
    flat = strip.reshape(bands, -1).T
    prediction[r:r + h, :] = model.predict(flat).reshape(h, cols)

profile.update(dtype="uint8", count=1, nodata=0)
with rasterio.open("landcover.tif", "w", **profile) as dst:
    dst.write(prediction, 1)

Open landcover.tif in QGIS, give each class a colour, and compare it with the original image. Your eyes are a powerful check, so look closely at places you know.

Which bands matter?

names = ["blue", "green", "red", "nir"]
for name, score in sorted(zip(names, model.feature_importances_), key=lambda x: -x[1]):
    print(f"{name:6s} {score:.2f}")

Feature importance shows which inputs the model leaned on. If one band dominates, adding related information, such as a vegetation index, often helps.

How to improve the map

• Add NDVI, a water index or texture measures as extra columns.
• Collect more labels where the confusion matrix shows mistakes, and in places that look different from your first samples.
• Use images from more than one date. Crops that look alike in one month often differ clearly across a season.
• Smooth the output with a majority filter to remove isolated stray pixels.
• For very detailed results, move to a deep learning model that understands shape and context, such as a U-Net.

Reporting accuracy responsibly

When you share a map, share how good it is. Report overall accuracy from spatial testing, and give per-class numbers, because a single average can hide a class that fails. Say which date and sensor the map comes from, and where it was checked on the ground. A map with a clear accuracy statement is far more trustworthy than one with a bigger number and no explanation.

Handling imbalance and noisy labels

Real label sets are rarely tidy. Some classes have thousands of examples, and others have a handful. Some labels are simply wrong. Both problems quietly cap your accuracy.

• Use class weights, as we did, or collect more examples of the rare class. Do not simply duplicate the few you have.
• Review a sample of labels by eye against the imagery. Even five percent of wrong labels can mislead a model.
• Merge classes that humans cannot reliably separate. Two classes that the labellers confuse will also confuse the model.
• Report results per class, so a poor result in a small class is not hidden by a good average.

Where deep learning fits

A random forest sees each pixel on its own. It knows nothing about what is next to it, so it cannot tell a rooftop from a parking lot by their shape. Convolutional neural networks look at patches of pixels, and learn shapes and textures, which is why they excel at detecting buildings, roads and fine boundaries. They need more labelled data, more computing power and more care in testing.

A sensible path is to start with a random forest, which gives you a baseline and teaches you about your data, then move to deep learning if the baseline is not good enough and you have the labels to support it.

Frequently asked questions

How many labelled pixels do I need?

Start with a few hundred per class, spread across the whole area. Add more where the confusion matrix shows the model struggling.

Should I scale my bands?

Random forests do not need it, because they split on thresholds. Models such as support vector machines and neural networks do.

Why does my map look noisy, with isolated pixels?

Per-pixel classification ignores neighbours. A majority filter, or an object-based approach that classifies groups of similar pixels, smooths the result.

Can I use the same model on another date or place?

Only after testing it there. Differences in season, sensor and land can reduce accuracy sharply. Retrain, or at least add labels from the new conditions.

You now know the full workflow: prepare labelled pixels, train, test honestly, map and improve. Real projects repeat this loop many times, and most of the gain comes from better labels, not a fancier model.