A satellite image goes in. A forest map comes out.

The Cisokan watershed in West Java is a mosaic of paddy fields, plantations, forest, and buildings — and from space, it's hard to tell them apart. This pipeline reads satellite images and automatically labels every patch of ground with one of nine land-cover classes.

The result is a raster map where every pixel answers: what is this piece of land?

Trace what happens when you run it

Imagine you open a terminal and type python planetscope_10epoch_obia_v3.py. Here is the exact journey your data takes — from raw satellite pixels to a coloured forest map on screen.

Nine answers the map can give

The pipeline uses a two-level hierarchy. Level 1 sorts land into 7 broad categories. Level 2 zooms into "Dense Vegetation" and splits it into three forest subtypes.

How well does it work?

Before diving into the machinery, here's the headline: the pipeline achieves 89.5% overall accuracy for the 7-class L1 map — verified on 228 held-out segments the model never saw during training.

Think of each script like a specialist on a film crew

No single person on a film set does everything — the director doesn't also operate the camera. Same idea here. Each Python script has one focused responsibility, and they hand work to each other in a clear order.

The scripts break into four groups: data downloaders, feature builders, classifiers, and analysis tools.

Three eyes in the sky

The classifier doesn't rely on a single satellite. It fuses three completely different data types — each one sees something the others can't.

The scripts on game day

On a fresh setup, the scripts run in this order, handing work to each other. Here's how they'd describe it to each other:

Put your knowledge to the test

Before moving on, check your understanding of which scripts do what — and when to run them.

Why does the model need 591 numbers per patch?

A single pixel on a single date tells you very little. A paddy field and a wet dirt road look nearly identical in March. But watch them both across 10 dates, measure their radar response, check how tall the vegetation is — now they diverge completely.

The art of this pipeline is turning each 93-pixel segment into a rich "fingerprint" — 591 numbers that describe everything observable about that patch from space.

The zonal stats trick — going from pixels to patches

Think of each segment as a territory on a map. You want to describe that territory — not pixel by pixel, but as a whole. The pipeline uses a clever counting technique called zonal statistics.

For each of the 316,000 segments and each of ~266 pixel features, it computes two numbers: the mean (typical value inside this territory) and the standard deviation (how varied the values are). That gives 266 × 2 = 532 zonal features.

count = np.bincount(lbl, minlength=n_seg+1)
s     = np.bincount(lbl, weights=val, minlength=n_seg+1)
sq    = np.bincount(lbl, weights=val*val, minlength=n_seg+1)
mean  = s / count
std   = np.sqrt(sq/count - mean**2)
            

Why compute indices instead of just using the raw bands?

Raw reflectance values are sensitive to lighting, atmosphere, and sensor quirks. Spectral indices are ratios that cancel out much of that noise and amplify the signal you care about.

Seasonality is the secret weapon

A single date can be ambiguous — a harvested paddy field and a barren slope look identical in the dry season. But across 10 dates, paddy has a very specific signature: it goes green fast (planting), peaks, then goes brown fast (harvest). Forest stays green year-round.

The pipeline computes temporal summary features that capture this rhythm:

Match the feature to what it detects

Drag each feature group to the land-cover property it is best at detecting.

Why not just classify pixel by pixel?

The naïve approach would be: for each of the 300 million pixels in the Cisokan image, decide its class. The problem? A single pixel is 3 m × 3 m — about the size of a car. A single pixel of "forest" could easily be a sunlit gap, a shadow, or a branch tip. Noise overwhelms signal.

OBIA solves this by first grouping nearby similar pixels into segments — patches of 50–500 pixels that share similar colour, texture, and brightness — then classifying each patch as a unit.

LSMS: the pixel sorter

LSMS works like a crowd-sorting algorithm at a festival. Each pixel "walks" toward the nearest cluster of similar-coloured neighbours. After all pixels settle, the groups become segments.

Three parameters control how fine-grained the result is:

The Random Forest — a forest that classifies forests

A Random Forest is not one model but an ensemble of 100 decision trees, each independently trained on a random subset of the data. The final prediction is the majority vote across all 100 trees.

Here is the exact configuration used in this pipeline:

RandomForestClassifier(
  n_estimators=100,
  min_samples_leaf=1,
  max_samples=0.5,
  bootstrap=True,
  max_features="sqrt",
  random_state=42,
  n_jobs=-1,
)
            

Two classifiers, two questions

The pipeline runs two Random Forests back to back. Think of it as a two-round interview: the first screener narrows the field, then a specialist makes the final call.

Debugging the classifier

You've deployed the model. A colleague reports that Natural Forest segments are frequently mislabelled as Production Forest. Apply what you've learned to diagnose the issue.

Overall Accuracy and Kappa — what the headline numbers mean

After training, the pipeline withholds 30% of labeled segments for a final exam — segments the model never saw during training. It predicts the class for each, then compares to the ground truth labels.

F1 Score — the per-class honest grade

Overall accuracy hides class-level failures. F1 score catches them: a class with F1 = 0.50 is getting half its predictions wrong, even if OA looks fine.

What the model actually learned to look at

Feature importance is like a spotlight on the model's attention. After training, you can ask: "Out of 591 features, which ones did the trees lean on most?" The answers are surprising — and actionable.

The L2 problem — and the path to fixing it

The L2 forest-subtype classifier achieves 57.5% accuracy on 40 validation segments. The dominant failure mode: 7 of 14 Production Forest segments get labelled as Natural Forest.

Here's the full confusion matrix:

Production Forest and Natural Forest both have tall, dense canopy. The features that should separate them — tree_height_std, S1_VV_stdDev — are present, but the model has only 48 Production training examples to learn from. The fix is straightforward: collect 30–50 more Production Forest training points from areas where planting rows are visible in the satellite basemap.

Final challenge — reading a real scenario