Building an Image Classification Pipeline for Satellite Imagery

Land-cover classification on satellite imagery is a small model and a very large loop. Doing the tile arithmetic before you write the loop is what stops a scene turning into a week of GPU time.

The tile arithmetic first

A Sentinel-2 Level-1C granule is, per the Copernicus Sentinel-2 product documentation, a 110 km × 110 km ortho-image in UTM, with four bands at 10 m, six at 20 m and three at 60 m. Take the 10 m bands and a 64-pixel tile, the patch size used by the EuroSAT land-use dataset:

scene edge    = 110 km / 10 m per pixel = 11,000 px  (nominal; the
                delivered raster is slightly smaller, tiles overlap)
tiles per row = floor(11,000 / 64) = 171
tiles         = 171 x 171 = 29,241 tiles per granule

with 50% overlap (stride 32):
tiles per row = floor((11,000 - 64) / 32) + 1 = 342
tiles         = 342 x 342 = 116,964 tiles per granule  (4.0x)

Now put a cost on it. If your classifier runs at an assumed 900 tiles/second on one GPU — an assumption you must replace with your own measurement, not a number to trust from a web page — a non-overlapping pass is 29,241 / 900 = 32 seconds of compute and the overlapped pass is 130 seconds. That is nothing. Multiply by a country: continental France is roughly 550,000 km², so 550,000 / 12,100 = 46 granules, and by twelve monthly revisits, and you have 550 granule-passes, 16 million tiles and about five hours of GPU per year of archive at that assumed rate.

The arithmetic is the deliverable here, not the seconds. Change the stride and you change the cost four-fold; change the tile size from 64 to 224 and you cut the tile count by a factor of twelve while making each tile more expensive. Do this on paper before you write the loop.

Reading windows instead of scenes

The one thing that must not happen is loading the scene into memory. A 11,000 × 11,000 raster in four bands at 16 bits is 11,000 × 11,000 × 4 × 2 bytes = 968 MB per granule before any float conversion, and float32 triples it. Cloud-optimised GeoTIFF exists precisely so you do not have to: internal tiling lets a reader fetch the bytes for one window and nothing else, over HTTP range requests if the file is in object storage.

Windowed reads are only fast if your windows line up with the file’s internal structure. A GeoTIFF written with 512-pixel internal tiles serves a 64-pixel window by reading the whole 512-pixel tile that contains it, so a 64-pixel window costs the same as a 512-pixel one. Reading eight aligned 64-pixel windows from one internal tile costs eight times what reading them in one pass would. Where the file is remote, the same arithmetic applies to HTTP range requests and the constant is much larger, because each request carries a round trip. Check src.block_shapes before choosing a tile size, and prefer a tile size that divides the internal block size rather than one that straddles it.

A striped GeoTIFF — one written in scanlines rather than tiles, which is still the GDAL default for some drivers — has no useful windowing at all: any window costs a full row of the raster. If the imagery you are handed is striped, re-tiling it once with gdal_translate -co TILED=YES before the run is cheaper than paying for it on every window.

The pipeline

1. Install the three things you need. pip install rasterio numpy torch. Rasterio wraps GDAL and gives you windowed reads and the affine transform; everything else here is standard.
2. Open the bands and confirm the geometry. Read src.transform, src.crs and src.shape before anything else. If the CRS is not what you expect, stop — a UTM zone mismatch between the imagery and your labels produces output that is entirely plausible and entirely in the wrong place.
3. Generate windows, do not slice arrays. rasterio.windows.Window(col_off, row_off, width, height) describes a read without performing it.
4. Drop cloudy tiles before the model sees them. Read the co-registered scene classification band and skip any tile whose cloud fraction exceeds your threshold. This is the cheapest speedup in the pipeline and it is also the difference between a land-cover map and a cloud map.
5. Batch, normalise, infer. Accumulate tiles into batches of 64 or 128 and run one forward pass per batch. Per-tile inference wastes most of the GPU on kernel launches.
6. Write predictions with their window origin so the class map can be assembled and georeferenced.

import numpy as np, rasterio, torch
from rasterio.windows import Window

TILE, STRIDE, BATCH = 64, 64, 128
CLOUD_CLASSES = (3, 8, 9, 10)   # Sentinel-2 SCL: shadow, cloud med/high, cirrus

model = torch.jit.load("landcover.pt").eval().cuda()

def tiles(width, height, tile=TILE, stride=STRIDE):
    for row in range(0, height - tile + 1, stride):
        for col in range(0, width - tile + 1, stride):
            yield Window(col, row, tile, tile)

def run(band_paths, scl_path, out_path, cloud_max=0.10):
    srcs = [rasterio.open(p) for p in band_paths]
    scl = rasterio.open(scl_path)
    ref = srcs[0]
    n_rows = (ref.height - TILE) // STRIDE + 1
    n_cols = (ref.width - TILE) // STRIDE + 1
    preds = np.full((n_rows, n_cols), 255, dtype=np.uint8)

    batch, coords = [], []
    for i, win in enumerate(tiles(ref.width, ref.height)):
        mask = scl.read(1, window=win, out_shape=(TILE, TILE))
        if np.isin(mask, CLOUD_CLASSES).mean() > cloud_max:
            continue
        arr = np.stack([s.read(1, window=win) for s in srcs]).astype("float32")
        batch.append(arr / 10000.0)          # S2 reflectance scale factor
        coords.append((i // n_cols, i % n_cols))

        if len(batch) == BATCH:
            preds = flush(model, batch, coords, preds)
            batch, coords = [], []
    if batch:
        preds = flush(model, batch, coords, preds)

    profile = ref.profile | dict(
        driver="GTiff", dtype="uint8", count=1, nodata=255,
        width=n_cols, height=n_rows,
        transform=ref.transform * rasterio.Affine.scale(STRIDE, STRIDE),
    )
    with rasterio.open(out_path, "w", **profile) as dst:
        dst.write(preds, 1)

def flush(model, batch, coords, preds):
    x = torch.from_numpy(np.stack(batch)).cuda()
    with torch.no_grad():
        y = model(x).argmax(dim=1).cpu().numpy().astype("uint8")
    for (r, c), label in zip(coords, y):
        preds[r, c] = label
    return preds

Writing the answer back to the map

The line that matters most in that listing is the transform:

transform = ref.transform * rasterio.Affine.scale(STRIDE, STRIDE)

The output raster has one pixel per tile, so its pixels are STRIDE times larger on the ground than the input’s. An affine transform is six numbers mapping pixel coordinates to projected coordinates, and scaling it by the stride is what keeps the class map aligned with the imagery. Get this wrong by writing the input transform unchanged and the class map is 64 times too small, anchored at the scene corner — which looks like a bug in the model and is a bug in one line of geometry.

Keep the CRS as the imagery’s UTM zone rather than reprojecting to WGS84 on write. Reprojection resamples, resampling a categorical raster with anything but nearest-neighbour invents classes that do not exist, and you will eventually want the metres that UTM gives you.

The four traps

• The reflectance scale factor. Sentinel-2 bands ship as integers scaled by 10,000. A model trained on 0–1 reflectance and fed raw integers sees inputs three to four orders of magnitude too large and will confidently predict one class everywhere. Divide, and check the product’s current scaling and any offset before assuming the constant.
• Band order. The model was trained on a specific band order. Nothing crashes if you supply red where the model expects near-infrared; the accuracy simply collapses. Pin the order in the same file as the weights.
• Tile-edge effects. A classifier fed a 64 × 64 window sees no context beyond it, so long thin features running across tile boundaries are classified inconsistently on each side. Overlapping tiles and averaging the class probabilities is the standard mitigation, and it costs the 4× computed above.
• Class imbalance in the training set. EuroSAT, from Helber and colleagues in 2019, holds 27,000 labelled Sentinel-2 patches across ten classes, which is a balanced research dataset. Real land cover is not balanced: a coastal scene may be 70% water, and a model tuned on balanced data over-predicts rare classes. Check the confusion behaviour on the classes you care about rather than the headline accuracy.

Before any of this, mask the clouds properly rather than with a class threshold, and understand what the mask is doing — see cloud masking before analysis. If the output you want is polygons rather than a class per tile, the pipeline shape changes to segmentation and vectorisation, described in building footprint detection.