Remote-Sensing Scene Classification — Reading the Land from Above

Chapter 1: The Task & Its Datasets

Let us make the task precise, because precision is what lets us build instead of wave. Remote-sensing scene classification takes a fixed-size image tile — a small square crop of satellite or aerial imagery — and outputs a land-use or land-cover label (or, as we'll see, a set of labels) from a fixed taxonomy. The unit of analysis is the tile, not the pixel and not a whole scene; you slice a giant satellite image into a grid of tiles and classify each one.

A word on the two phrases you'll hear. Land cover is the physical material on the ground — water, trees, bare soil, concrete. Land use is the human function — residential, industrial, agricultural. They overlap but differ: a grass field (cover) might be a park, a lawn, or a golf course (use). Datasets vary in which they target, and this matters because land use is harder — it asks the model to infer human intent from physical appearance.

Three datasets define the field. Anchor on these three names; they recur in every paper.

EuroSAT is where almost everyone starts. It cuts 27,000 small tiles from Sentinel-2 imagery and labels each with one of ten covers: AnnualCrop, Forest, HerbaceousVegetation, Highway, Industrial, Pasture, PermanentCrop, Residential, River, SeaLake. Because each tile is forced to one label, EuroSAT is single-label and easy to score — perfect for building intuition, which is why we'll use its ten classes throughout this lesson.

NWPU-RESISC45 raises the stakes to 45 land-use classes from higher-resolution aerial RGB, with the fine distinctions that make scene use hard: it separates "freeway" from "overpass," "church" from "palace." It is the bird's-eye analogue of an indoor-scene benchmark — many classes, subtle differences.

BigEarthNet is the realistic one and the reason Chapter 6 exists. With ~590,000 Sentinel-2 tiles, it abandons the single-label fiction: each tile is annotated with all the land covers present, drawn from a 19- or 43-class nomenclature. A tile can be simultaneously "Coniferous forest," "Pastures," and "Water bodies." Train on BigEarthNet and you are forced to confront multi-label classification head-on.

The widget below lays out the ten EuroSAT classes as a small gallery of stylized overhead tiles. Tap any class to see a one-line definition and the kind of texture a satellite records for it. Notice how visually close some are — AnnualCrop, Pasture, and PermanentCrop are all "green-ish farmland," separated more by texture and season than by color. That closeness is the difficulty we'll quantify in Chapter 4.

Chapter 2: The Multispectral Input — Beyond RGB

Here is the heart of what makes remote sensing different from ordinary vision: the input is not a 3-channel RGB image but a multispectral stack of many channels, each measuring a different slice of the light spectrum. "Multispectral" means the sensor splits incoming light into several named bands — intervals of wavelength — and records a separate grayscale image for each. Sentinel-2 records 13 bands, spanning visible light, near-infrared, and shortwave infrared.

Why bother? Because different surfaces reflect different wavelengths in characteristic ways, and the most useful differences live outside the visible range. The signature example is the near-infrared (NIR) band, light just beyond red that your eyes can't detect. Healthy, chlorophyll-rich vegetation reflects NIR extremely strongly — far more than it reflects green — while reflecting red light weakly because chlorophyll absorbs red for photosynthesis. So in NIR, a thriving crop field is almost glowing, while bare soil, water, and concrete stay dim. The vegetation/non-vegetation question, ambiguous in RGB, becomes nearly trivial in NIR.

Think of the input as a channel stack: a tensor of shape (bands, height, width). For a 64×64 Sentinel-2 tile with all bands, that's roughly (13, 64, 64) — thirteen stacked grayscale images, perfectly aligned, one per wavelength. (In practice the bands come at different native resolutions — 10 m, 20 m, 60 m ground sampling — and are resampled to a common grid.) Where an RGB photo gives a network three numbers per pixel, a multispectral tile gives it thirteen, and the extra ten are where much of the discriminative power lives.

Rather than feed all thirteen raw bands to a network and hope, remote-sensing practice leans on spectral indices: simple arithmetic combinations of two bands that isolate a physical property. The most famous is the NDVI, the Normalized Difference Vegetation Index, a single number per pixel that measures "how vegetated is this." Its formula combines the NIR and red bands:

Read the structure. The numerator (NIR − RED) is large and positive exactly when NIR is high and red is low — the chlorophyll signature. The denominator (NIR + RED) normalizes by overall brightness, so a sunlit field and a shaded field of the same vegetation give similar NDVI — the index measures type, not brightness. NDVI ranges from −1 to +1: dense healthy vegetation lands around 0.6–0.9, bare soil near 0.1–0.2, and water goes negative because water reflects almost no NIR.

The widget makes NDVI tactile. Drag the NIR and red sliders and watch the resulting NDVI value and its color band. Push NIR up with red low and you climb toward dense vegetation; pull NIR below red and the index goes negative into the water regime. This single derived channel often separates EuroSAT's classes better than any raw band, and it is the substrate for one of our Code Labs.

Chapter 3: The CNN Baseline

Now we build a classifier. The remote-sensing baseline is the same workhorse as ordinary scene recognition — a convolutional neural network — with one surgical change at the input. Let us trace the data flow with real shapes, because the shapes show exactly where multispectral enters.

In ordinary vision the input is (3, 224, 224): three RGB channels. Here the input is the multispectral stack, say (13, 64, 64): thirteen bands over a 64×64 tile. The only structural difference is the first convolutional layer. A standard CNN's first conv expects 3 input channels; we widen it to accept 13 (or however many bands we use). Concretely, the first layer's weight tensor changes from shape (out, 3, k, k) to (out, 13, k, k) — each output filter now mixes thirteen input bands instead of three. Everything downstream is unchanged.

From there the network behaves like any CNN. Blocks of convolutions shrink the spatial size and grow the channel count, until you hold a feature map of shape, say, (256, 4, 4): 256 learned feature channels, each a tiny 4×4 spatial grid. Each channel fires where some learned pattern appears — "bright NIR texture here," "regular grid of building roofs there," "smooth dark water over here."

Then the move that matters for tiles: global average pooling (GAP). We collapse each spatial grid to a single number by averaging it. The (256, 4, 4) map becomes a (256,) vector. This vector says, for each of the 256 patterns, "how much of this pattern is present anywhere in the tile," discarding where. For overhead land cover this is exactly right: a forest tile is a forest tile whether the densest trees are top-left or bottom-right. Location within the tile is mostly irrelevant; presence is everything.

Finally a fully-connected layer maps the (256,) vector to K logits (one per land-cover class), and softmax turns logits into probabilities. That is the whole single-label classifier. Let us watch one toy example end to end.

The simulation below runs this pipeline live. Click cells in the feature grid to toggle which "patterns" are firing (think: NIR-bright cells, built-up cells), and watch GAP collapse the grid to a value, a tiny linear layer produce logits over four land covers, and softmax produce probabilities. Paint a "forest signature" or a "residential signature" and see what it classifies as.

Here is the same pipeline as code — first the multispectral first layer from scratch, then the library version, so you see they do the same thing:

python
import numpy as np

# tile: (C_in, H, W) multispectral stack, e.g. (13, 64, 64)
# the ONLY change vs an RGB model is C_in = 13 instead of 3
def global_average_pool(feature_map):
    return feature_map.mean(axis=(1, 2))   # (C,) one number per channel — location discarded

def softmax(z):
    z = z - z.max()                      # subtract max for numerical stability
    e = np.exp(z)
    return e / e.sum()

# W_fc: (K, C) classifier weights; b_fc: (K,) biases; K = number of land covers
def classify_tile(feature_map, W_fc, b_fc):
    v      = global_average_pool(feature_map)   # (C,) the tile descriptor
    logits = W_fc @ v + b_fc                     # (K,) raw class scores
    probs  = softmax(logits)                     # (K,) sum to 1 (single-label)
    return probs.argmax(), probs                 # top-1 index, full distribution

python
# Library version: take a standard ResNet and widen its first conv to 13 bands
import torch.nn as nn, torchvision.models as models
net = models.resnet18(num_classes=10)        # 10 EuroSAT classes
# swap the 3-channel stem for a 13-channel one — the whole multispectral trick:
net.conv1 = nn.Conv2d(13, 64, kernel_size=7, stride=2, padding=3, bias=False)
# resnet internally does conv blocks -> AdaptiveAvgPool2d(1) [== GAP] -> Linear
logits = net(tile_batch)                  # (B, 10)
probs  = logits.softmax(dim=1)            # same softmax we wrote by hand

See the correspondence: nn.Conv2d(13, …) is the multispectral input change, AdaptiveAvgPool2d(1) is our global_average_pool, nn.Linear is our W_fc @ v + b_fc. A plain ResNet-18 with a 13-band stem, trained from scratch, reaches well over 95% accuracy on EuroSAT — the task is genuinely learnable once you let the model see the spectral bands.

Chapter 4: Why It's Genuinely Hard

EuroSAT looks easy — 95%+ accuracy — but that number is a trap, because it comes from small, clean, single-label tiles. Real deployment is humbling, and three distinct forces make remote-sensing scene classification much harder than the toy benchmark suggests. Naming each one tells you what to engineer around.

Force 1: geographic intra-class variance. A "Residential" tile in suburban Germany — red-tile roofs, gardens, narrow streets — looks nothing like a "Residential" tile in Phoenix — flat tan roofs, pools, wide grids — or in Lagos, or rural India. The same class spans wildly different appearances across continents, climates, and building styles. A model trained on European tiles can collapse on African ones, not because it's weak but because it never saw that variety. The class concept is abstract; the pixels are parochial.

Force 2: resolution and ground sample distance. The key quantity is GSD, the Ground Sample Distance — how many metres of real ground each pixel covers. Sentinel-2's best bands are 10 m GSD: one pixel is a 10×10 m patch. At that scale, individual houses are a pixel or two and small features vanish into texture. Aerial imagery at sub-metre GSD shows individual cars. A model and a class definition that work at 10 m can be meaningless at 0.5 m and vice versa — "industrial" is a texture at coarse GSD and a collection of recognizable buildings at fine GSD. Resolution silently changes what each class even is.

Force 3: the multi-label reality. The single biggest gap between EuroSAT and the real world: a one-kilometre tile almost never contains one land cover. It contains forest and pasture and a road and some water. Forcing a single label is not just lossy — it makes the labels themselves inconsistent, because annotators must pick "the dominant one," and dominance is subjective. This is why BigEarthNet went multi-label, and it's the subject of Chapter 6.

Watch the geographic-shift force directly. The widget shows a classifier's confidence as you "fly" it from its training region (Europe) to increasingly different regions. As the geographic distance grows, the same "Residential" tiles look less and less like training data, and accuracy slides — the classic distribution-shift collapse. This is the number that benchmark accuracy hides.

Let us put a number on the multi-label problem so it isn't just a feeling. Suppose a real tile genuinely contains, by area, 50% Forest, 30% Pasture, 15% Water, 5% Highway. A single-label dataset forces the annotator to write "Forest." Your softmax model, trained on such labels, learns to output "Forest" too — and is marked correct. But it has been actively trained to suppress Pasture, Water, and Highway, even though all three are truly present. The single-label format doesn't just lose information; it teaches the model to deny reality. Multi-label fixes this by scoring each cover independently, which we build in Chapter 6.

Chapter 5: Geospatial Foundation Models

Here is the problem the field bumped into: labelling satellite tiles is expensive (you need experts, or ground truth), so labelled data is scarce — but unlabelled satellite imagery is essentially infinite. Petabytes of Sentinel-2 pixels stream down every week, free, with no labels. The obvious move: learn from the unlabelled flood first, then fine-tune on the small labelled set. That move is self-supervised pretraining, and it produced the geospatial foundation models that now define the state of the art.

The dominant recipe is the masked autoencoder (MAE). The idea is gorgeous in its simplicity: take a tile, hide a large fraction of it (mask out, say, 70% of the patches), feed the network only the visible part, and train it to reconstruct the hidden pixels. No labels needed — the supervision is the image itself. To fill in a masked patch correctly, the network is forced to learn what land covers look like, how they're spatially arranged, and how the spectral bands relate. Those learned representations then transfer beautifully to classification.

Three geospatial foundation models are worth knowing by name, each adapting MAE to a different axis of satellite data:

What makes these geospatial rather than generic MAEs is that they mask along axes a normal image model doesn't have. A photo MAE masks spatial patches. SatMAE additionally masks entire spectral bands (predict the NIR band from the others) and entire time steps (predict July's tile from May's). SpectralGPT goes further, treating the spectrum as a dimension to be masked and reconstructed in 3D cubes. Prithvi, from NASA and IBM, is a temporal Vision Transformer MAE pretrained on massive Sentinel/Harmonized-Landsat-Sentinel data and released openly. The unifying theme: exploit the structure satellites give you for free — many bands, repeated passes — as the pretraining signal.

Once pretrained, a geospatial foundation model is used through the same three modes as any foundation model — and the choice is a real engineering decision driven by how many labels you have:

The widget below sketches the payoff that makes pretraining worth the trouble: how accuracy scales with the number of labelled tiles, for a model trained from scratch versus a pretrained foundation model that's only fine-tuned. The pretrained curve starts far higher and stays ahead — especially in the low-label regime, which is exactly where real deployments live (you rarely have a million labels for your specific region). Drag the label budget and watch the gap.

Chapter 6: Multi-Label — A Tile Is Many Covers

This is the chapter the whole lesson has been building toward, because it's where the machinery genuinely changes. We've said a real one-kilometre tile contains several land covers at once. Now we make the model honor that, and it requires abandoning the most reflexive habit in classification: softmax.

Recall what softmax does: it forces the class probabilities to sum to 1. That is a statement — "exactly one class is the answer; raising one probability must lower the others." For a tile that is genuinely Forest and Pasture and Water, that constraint is a lie. If Forest's probability goes up, softmax mechanically pushes Water's down, even though both are truly present. Softmax models competition between classes. Multi-label land cover is not a competition.

The fix is to score each land cover independently with a sigmoid. A sigmoid squashes a single logit into [0,1] on its own, with no coupling to the other classes. So instead of one softmax over K classes, you run K separate sigmoids — each answering "is this cover present? yes/no" — and the outputs need not sum to 1. A tile can be 0.9 Forest, 0.8 Pasture, 0.7 Water simultaneously. The training loss changes to match: binary cross-entropy applied per class (each cover is its own little present/absent classifier), replacing the single cross-entropy over a softmax.

Read the contrast carefully. Softmax: P(class k) = e^z_k / Σ_j e^z_j — coupled, sums to 1, "pick one." Sigmoid: P(class k present) = 1 / (1 + e^−z_k) — independent, each in [0,1], "each present or not." This single swap is the entire multi-label shift.

With sigmoids, prediction has one more knob: a per-class threshold. Each cover is "predicted present" if its sigmoid output exceeds a threshold — commonly 0.5, but tunable per class. This matters because covers differ in prevalence: a rare cover (say "Beaches") may need a lower threshold to be caught at all, while a common one ("Pastures") can afford a higher one. Tuning thresholds per class is how you trade off catching every instance (recall) against avoiding false alarms (precision) for each land cover separately — impossible under a single softmax argmax.

The widget below shows it live. You set three logits and toggle between softmax and sigmoid. Under softmax, the three bars are chained — pushing one up drags the others down, and exactly one "wins." Under sigmoid, each bar moves on its own, and you draw a threshold line that decides the predicted set. This is the clearest way to feel why multi-label needs sigmoid.

Chapter 7: Showcase — A Live Multispectral Tile Classifier

Everything converges here. This is a working multi-label land-cover classifier you can drive and, more importantly, break. You set the tile's spectral character with sliders — near-infrared brightness, redness, water-darkness, the built-up texture, and an overall haze/noise knob — and the model computes a logit for each of five covers, then sigmoids them into independent presence probabilities. A threshold line decides the predicted set. The bars are the live per-cover probabilities; the panel shows which covers cross the threshold.

The five covers — Forest, Pasture, Water, Residential, Highway — each respond to the spectral sliders the way real surfaces do. Forest loves high NIR and a vegetation index; Water loves low NIR and darkness; Residential loves the built-up texture; Highway loves a thin built-up streak with low vegetation. No single slider decides a cover — the combination of spectral cues does, exactly as Chapters 2 and 3 argued.

Three experiments, each illustrating a chapter:

• Make a clean forest (high NIR, low red, low water-dark, low built-up). Watch Forest's bar dominate and cross the threshold alone — the easy, separable case from Chapter 3.
• Now raise built-up texture and keep NIR high. Watch Forest and Residential both cross the threshold — you've manufactured the multi-label reality of Chapter 6, live: a tile that is genuinely two things at once, which softmax could never express.
• Crank the haze/noise slider. This blurs the spectral signal the way clouds, shadow, or a hard tile would. Confident covers drift toward the threshold and the predicted set becomes unstable — you're watching the geographic/atmospheric difficulty of Chapter 4 erode the decision in real time.

If you removed this simulation, would you lose understanding? Yes — reading "multi-label means several covers can be present" is abstract until you watch Forest and Residential both light up by sliding one texture knob, with no winner forced. The napkin drawing of this whole lesson is exactly this: spectral cues in, an independent probability per cover out, a threshold drawing the predicted set, haze pulling it toward chaos.

Chapter 8: Evaluation — Scoring Multi-Label Honestly

You understand the models; here is how to judge them without being fooled — and multi-label evaluation has its own culture, because plain accuracy quietly lies. Three ideas cover most of it: why accuracy fails, mean average precision (mAP), and geographic generalization.

Why "accuracy" breaks for multi-label. With single-label tiles you compute accuracy = fraction predicted exactly right. But with 19 independent covers, you could score "accuracy" as fraction of the 19 yes/no decisions you got right — and that number is misleadingly high, because most covers are absent on any given tile. If a tile has 3 of 19 covers and your model predicts "all absent," it's wrong on only 3 of 19 = 84% "accurate" while having found nothing. The class imbalance (mostly-absent labels) makes accuracy a useless flatterer. We need a metric that focuses on getting the present covers right.

Mean Average Precision (mAP). The standard multi-label metric. For each class separately, you rank all tiles by that class's sigmoid score, then compute Average Precision (AP) — the area under that class's precision-recall curve, which rewards putting the truly-present tiles at the top of the ranking. Then you average AP across all classes to get mAP. Because it's computed per class and averaged, mAP treats a rare cover as seriously as a common one, and because it's threshold-free (it scores the whole ranking), it isn't fooled by a single bad threshold choice. BigEarthNet leaderboards report mAP precisely for these reasons.

The same metric in code, computed the way a benchmark script does — first the per-class AP from scratch, then the one-liner:

python
import numpy as np

# scores: (N,) model sigmoid output for one class; truth: (N,) 1 if present else 0
def average_precision(scores, truth):
    order = np.argsort(-scores)            # rank tiles best-first by this class's score
    t = truth[order]
    hits = np.cumsum(t)                    # correct-so-far at each rank
    ranks = np.arange(1, len(t) + 1)     # 1,2,3,...
    precision_at_rank = hits / ranks       # precision at every position
    return (precision_at_rank * t).sum() / max(t.sum(), 1)  # average precision AT the hits

# mAP = mean of per-class AP over all K classes
def mean_ap(scores_KxN, truth_KxN):
    return np.mean([average_precision(scores_KxN[k], truth_KxN[k]) for k in range(len(scores_KxN))])

# sklearn one-liner equivalent (micro/macro averaged):
#   from sklearn.metrics import average_precision_score
#   mAP = average_precision_score(truth_NxK, scores_NxK, average='macro')

Notice the ranking, not a threshold, drives AP — that's why mAP is robust to a poorly chosen cutoff. Geographic generalization is the silent killer, again. Just as in Chapter 4, the honest question is not "what's the mAP on the held-out split?" but "what's the mAP on tiles from a region the model never trained on?" A model can post 90% mAP on a random European split and far less on an African or Asian split. Serious remote-sensing evaluation reports cross-region or geographically-disjoint splits, where training and test tiles come from different continents, because that's what deployment actually demands. An in-distribution number is an optimistic ceiling.

The interactive panel below is a small multi-label scoreboard. Each row is a tile; each column a cover; green = truly present, and you toggle the model's predictions to see precision, recall, and per-class AP respond. Make the model predict "all absent" and watch accuracy stay high while recall collapses — the exact trap that makes accuracy useless here.

Chapter 9: Connections & Where It Lives

Remote-sensing scene classification is a hub, not an island. The exact pattern you built — read a whole tile, summarize it into a descriptor, score a fixed taxonomy, and accept that the world is multi-label — recurs across vision, audio, and 3D sensing. Naming those connections turns one lesson into a mental map.

The scene-classification family. This lesson is one sibling of a larger family that all ask "what is this whole signal?" Ordinary Visual Scene Classification names "kitchen" or "beach" in a ground-level photo — same holistic-descriptor machinery, different viewpoint. Acoustic Scene Classification asks "is this audio a park, a metro, an airport?" — whole-signal classification over a spectrogram. RGB-D scene classification adds a depth channel that sharply helps indoor layout. Each is "this lesson, new sensor": the descriptor-then-classify skeleton is shared; what changes is the input modality and what "a class" means.

The multispectral input is the distinctive twist here, and it connects to the broader idea of feeding non-RGB channels through a CNN — the same first-layer-widening trick you'd use for depth, thermal, or hyperspectral. The masked-autoencoder pretraining of SatMAE/SpectralGPT/Prithvi is the satellite cousin of the self-supervised pretraining you'd recognize from Vision Transformers and masked-image modelling generally. And the zero-shot branch of foundation models runs through Contrastive CLIP: remote-sensing CLIP variants (like RemoteCLIP) let you classify a tile by writing "a satellite photo of a forest" — no labels, editable taxonomy.

Cheat sheet

"We shall not cease from exploration, and the end of all our exploring will be to arrive where we started and know the place for the first time." — T. S. Eliot. Remote sensing is the discipline of knowing a place — from a viewpoint no human ever had.