Raster Machine Learning & Model Inference

Fitting a model to satellite imagery is rarely the hard part. The hard parts are getting pixels and labels onto the same grid without silently shifting one against the other, splitting the data so the reported accuracy means something, and then running the trained model across scenes that are far larger than memory while keeping every output pixel georeferenced. This section covers that whole loop for remote sensing practitioners who already know their way around scikit-learn or PyTorch but need the raster half to be correct.

Everything here assumes the grid mechanics from Core Raster Fundamentals & STAC Mapping and the preprocessing from Satellite Processing Workflows & Index Pipelines. When a model has to run over a continent rather than a county, the execution patterns in Cloud Execution & Orchestration are what carry it.

---

The Shape of a Raster Modelling Pipeline

A raster modelling pipeline is a loop, not a line. Labels produce chips, chips produce a model, the model produces a raster, and the raster is scored against reference data that usually sends you back to the labels. Every arrow in that loop is a place where a geotransform can drift.

The loop has one non-negotiable invariant: every array that enters or leaves it carries a transform, a CRS, and a nodata policy. A chip is not a 256×256 array, it is a 256×256 array plus the affine transform of the window it came from. Drop that transform at any stage and the prediction raster cannot be written back onto the scene without guesswork.

---

Key Components

Two of these deserve emphasis. window_transform is the function that keeps georeferencing alive across the model boundary, and rasterize is the function that guarantees labels and pixels share a grid. Most of the failures in this section trace back to one of them being skipped.

---

Production Patterns

Pattern 1 — A chip reader that never loses georeferencing

Every training set in this section is built by the same primitive: take a geometry, compute the window that centres it, read the bands, and keep the window transform alongside the array. The details of label burning are covered in building training datasets from satellite imagery; this is the core read.

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


def read_chip(src: rasterio.DatasetReader, row: int, col: int, size: int = 256):
    """Read a square chip and return the array with the transform that locates it."""
    win = Window(col_off=col, row_off=row, width=size, height=size)
    # boundless=True pads with fill rather than raising when the window
    # runs off the edge of the scene — essential near scene borders.
    arr = src.read(window=win, boundless=True, fill_value=src.nodata or 0)
    return arr.astype("float32"), src.window_transform(win)


with rasterio.open("s2_stack_36NYF.tif") as src:
    chip, chip_transform = read_chip(src, row=4096, col=2048)
    print(chip.shape, chip_transform)

The boundless=True argument is what makes edge chips safe. Without it, a window that extends past the scene raises or silently clips, and a clipped chip that is 256×241 instead of 256×256 will crash a fixed-input model halfway through a long run.

Pattern 2 — Pixel-based prediction as a reshape

For non-convolutional models, inference over a block is three lines: read, reshape to (pixels, features), predict, reshape back. The band order used here must be the band order used at training time, which is why the feature list belongs in the model artifact rather than in a comment.

import numpy as np


def predict_block(model, block: np.ndarray, nodata_mask: np.ndarray) -> np.ndarray:
    """Apply a pixel-based sklearn model to a (bands, rows, cols) block."""
    bands, rows, cols = block.shape
    flat = block.reshape(bands, rows * cols).T          # (pixels, features)
    out = np.zeros(rows * cols, dtype="uint8")          # 0 reserved for nodata
    valid = ~nodata_mask.reshape(-1)
    if valid.any():
        out[valid] = model.predict(flat[valid]) + 1     # classes shifted off 0
    return out.reshape(rows, cols)

Reserving class 0 for “no prediction” costs one class value and saves an enormous amount of downstream confusion: a pixel that was cloud, off-scene, or missing a band is now distinguishable from a pixel the model genuinely assigned to the first class. The masking inputs usually come from cloud and shadow masking strategies.

Pattern 3 — Streaming inference over a whole scene

The full-scene loop writes block by block so peak memory stays at one window rather than one scene. It is the same windowed-read discipline described in optimizing rasterio window reads for memory efficiency, with a model call in the middle.

import rasterio
from rasterio.windows import Window


def infer_scene(model, src_path: str, dst_path: str, block: int = 1024) -> None:
    with rasterio.open(src_path) as src:
        profile = src.profile | {
            "count": 1, "dtype": "uint8", "nodata": 0,
            "compress": "deflate", "tiled": True,
            "blockxsize": 512, "blockysize": 512,
        }
        with rasterio.open(dst_path, "w", **profile) as dst:
            for row in range(0, src.height, block):
                for col in range(0, src.width, block):
                    win = Window(col, row,
                                 min(block, src.width - col),
                                 min(block, src.height - row))
                    arr = src.read(window=win).astype("float32")
                    mask = (arr == (src.nodata or 0)).any(axis=0)
                    dst.write(predict_block(model, arr, mask), 1, window=win)

Because the destination profile is derived from the source profile, the prediction inherits the source CRS and transform exactly — the output is aligned with the input by construction rather than by a later reproject_match call. Convolutional models need an overlap on top of this loop, which is the subject of running model inference over large rasters.

---

Choosing Between Pixel-Based and Convolutional Models

The choice of model class is really a choice about where the signal lives. If a pixel’s own spectrum tells you the answer — water is dark in the near infrared, bare soil is bright in the shortwave — a pixel-based model is not a compromise, it is the correct tool. If the answer depends on shape, texture or context — this bright rectangle is a greenhouse, that dark strip is a road rather than a shadow — no amount of feature engineering on a single pixel will recover it, and a convolutional model earns its complexity.

The speckle in the last row is worth calling out because it is so easy to fix and so often left in. A pixel-based classifier has no reason to prefer spatially coherent output, so single misclassified pixels are scattered across otherwise clean fields. A majority filter over a 3×3 or 5×5 neighbourhood removes most of them without touching real boundaries, and the technique is set out in smoothing and post-processing classification rasters.

There is a middle path that is under-used: a pixel-based model on features that encode context. Adding a focal mean, a focal standard deviation, and a couple of terrain derivatives gives a gradient boosting model most of what a shallow convolution would extract, at a fraction of the labelling and training cost. That construction is the subject of feature engineering for pixel-based models, and for many operational land-cover problems it is the right answer.

---

The Model Artifact as a Data Contract

A trained model that ships as a bare pickle file is a liability. Six months later nobody remembers which bands went in, in what order, whether reflectance was scaled, what the class codes meant, or which nodata value was masked. Every one of those is required to reproduce a prediction, and none of them are recoverable from the weights.

Treat the artifact as a contract: the estimator plus the exact preprocessing it assumes, written down in a form the inference script reads rather than a form a human reads.

import json
from pathlib import Path

import joblib


def save_model_artifact(model, path: str, *, bands: list[str],
                        scale: float, nodata: float,
                        classes: dict[int, str]) -> None:
    """Persist an estimator together with the contract its inputs must satisfy."""
    out = Path(path)
    out.mkdir(parents=True, exist_ok=True)
    joblib.dump(model, out / "estimator.joblib")
    (out / "contract.json").write_text(json.dumps({
        "band_order": bands,          # names, not indices — indices drift
        "reflectance_scale": scale,   # divide raw DN by this before predicting
        "nodata": nodata,             # pixels equal to this are never predicted
        "classes": classes,           # code -> human label, code 0 = no prediction
        "chip_size": 1,               # 1 for pixel models, N for chip models
    }, indent=2))


def load_model_artifact(path: str):
    out = Path(path)
    contract = json.loads((out / "contract.json").read_text())
    return joblib.load(out / "estimator.joblib"), contract

The inference script then asserts the contract against the file it has been handed rather than trusting the caller:

import rasterio

model, contract = load_model_artifact("models/landcover_v3")

with rasterio.open("scene_stack.tif") as src:
    names = [src.descriptions[i] for i in range(src.count)]
    if names != contract["band_order"]:
        raise ValueError(
            f"band order mismatch: file has {names}, "
            f"model expects {contract['band_order']}"
        )

That assertion has caught more silent accuracy losses than any amount of cross-validation, because a band-order swap produces a prediction raster that looks entirely plausible — coherent regions, sensible boundaries — and is wrong everywhere. Writing band descriptions into the stack in the first place is covered in reading and writing GDAL tags and band descriptions.

The same contract travels with the output. A prediction raster that records the model version, the contract hash, and the input scene identifiers in its GDAL metadata can be traced back to the run that produced it; one that does not is a picture, not a data product. Attaching model metadata to STAC items takes that further and makes the output discoverable.

---

How This Section Fits Together

The five topics below are ordered the way a project moves, and each one assumes the previous one produced something well formed.

Building training datasets from satellite imagery is where labels and pixels are forced onto a common grid. It is the least glamorous stage and the one that determines the ceiling on everything after it.

Feature engineering for pixel-based models turns raw bands into a matrix worth fitting — spectral indices, terrain derivatives, and temporal summaries stacked in a fixed, documented order.

Running model inference over large rasters is the engineering stage: overlapping windows, batched feeds, and a write loop that keeps the output aligned with the input.

Validating raster predictions and accuracy assessment produces the numbers that other people will rely on, including area estimates with honest confidence intervals rather than a pixel count.

Exporting and serving model outputs closes the loop by writing a Cloud-Optimized GeoTIFF that the rest of the stack — tile servers, catalogs, downstream pipelines — can consume without a conversion step.

---

Common Pitfalls & Failure Modes

Label–pixel misalignment. Labels digitised on a basemap in WGS84 and burned onto a UTM scene without reprojection land half a field away. Reproject the geometries into the raster CRS before rasterizing, using the patterns in mastering CRS transformations in rasterio. The symptom is a model that trains to 0.99 on spectrally meaningless features.

Optimistic validation from random splits. The single most common defect in published remote sensing models. Spatial autocorrelation means a random 80/20 split leaves near-duplicate pixels on both sides. Split by block, tile, or region instead — the mechanics are in designing spatial cross-validation for raster models.

Nodata entering the feature matrix. A fill value of −9999 or 0 that reaches the model is treated as a real reflectance. Tree ensembles will happily learn “if band 4 is exactly 0, class is water”. Mask before flattening, always.

Band order drift between training and inference. The training script read bands [B04, B08, B11]; the inference script reads src.read() and gets whatever order the file happens to store. Pin the band list by name in the model artifact and assert it at load time.

Scale factor applied once, or twice. Sentinel-2 L2A reflectance is stored as scaled integers. Apply the factor at training and forget it at inference and every feature is 10,000× off. The rules are set out in handling nodata and scale factors in band math.

Tile seams in convolutional output. A segmentation model applied to abutting tiles produces visible discontinuities at every tile edge, because pixels near a tile border have no context. Predict with overlap and keep only the centre, exactly as a mosaic keeps only the interior of each scene in seamless mosaicking and edge blending.

---

Performance and Scale

Three numbers govern how a raster model scales.

Read throughput. Inference over a scene reads every pixel exactly once, so the floor on runtime is set by how fast bytes arrive. A well-formed Cloud-Optimized GeoTIFF read with a matching block size is several times faster than a striped GeoTIFF read block by block; the measurements are in benchmarking COG read throughput from object storage.

Block size versus model batch. The natural read unit is the file’s internal tile — usually 512×512. The natural model unit is a batch of chips. Making the read block an exact multiple of the chip size means no chip ever straddles two reads, which removes an entire class of buffering code.

Parallelism model. Pixel-based models parallelise perfectly by block, so a process pool over windows scales linearly until the storage saturates. Convolutional models on a single GPU do not: the correct shape there is one reader thread feeding a queue while the GPU drains it. When the scene count grows past what one machine can finish overnight, move the block loop onto Dask or a batch queue rather than making the single-machine loop cleverer.

Memory. Peak memory is bands × block² × 4 bytes × concurrency, plus the model. Six bands at 1024² float32 is 25 MB per worker — comfortable. The same at 4096² is 400 MB per worker, which is how a sixteen-way pool exhausts a 32 GB machine.

---

Frequently Asked Questions

Q: Do I need a GPU to run models over satellite rasters? Not for pixel-based models. Gradient boosting and random forests on a handful of bands run comfortably on CPU and are usually I/O bound rather than compute bound. A GPU only becomes necessary for convolutional segmentation models where each tile costs tens of milliseconds of matrix work; even then the raster read is often the bottleneck, so overlap reads with inference before buying hardware.

Q: Why do my validation scores collapse when the model is applied to a new region? Almost always because the validation split was random rather than spatial. Neighbouring pixels are highly autocorrelated, so a random split puts near-identical samples in both train and validation, and the score measures memorisation. Splitting by block, tile, or administrative unit gives a number that survives contact with a new scene.

Q: Should predictions be written as classes or probabilities? Write both when storage allows. A uint8 class raster is what most downstream consumers want, but the probability or logit band is what lets you re-threshold later without re-running the model, and it is essential for any uncertainty reporting. Store probabilities as a scaled integer band to keep the file small.

Q: How much training data does a pixel-based land cover model need? Far less than a segmentation model, but more spread than people expect: a few thousand pixels per class drawn from many separate locations beats a hundred thousand pixels drawn from six polygons. Diversity of location dominates raw sample count, because the second polygon in the same field adds almost no information.

---

Related

• Building Training Datasets from Satellite Imagery — chip extraction, label rasterization, and spatially disjoint folds.
• Running Model Inference over Large Rasters — overlapping windows, seam-free stitching, and batched GPU feeds.
• Feature Engineering for Pixel-Based Models — turning bands, indices, terrain and time into a stable feature matrix.
• Validating Raster Predictions and Accuracy Assessment — confusion matrices, spatial folds, and defensible area estimates.
• Exporting and Serving Model Outputs — writing prediction COGs, vectorizing masks, and recording provenance.