Building Training Datasets from Satellite Imagery

A training dataset for a raster model is not a folder of images. It is a set of arrays, each one carrying the affine transform that locates it, paired with a label mask burned on exactly the same grid, tagged with the fold it belongs to and the scene it came from. Get that structure right and the modelling that follows is ordinary machine learning; get it wrong — a half-pixel offset here, a random fold assignment there — and no amount of architecture search will recover the signal. This topic sits inside Raster Machine Learning & Model Inference and feeds everything downstream of it.

The specific challenge is that the two halves of a training sample arrive in completely different coordinate systems. Imagery arrives as a pixel grid defined by a geotransform; labels arrive as vector geometries in whatever CRS the person digitising them happened to use. Bringing them together is a projection problem, a rasterization problem, and a sampling problem, in that order.

---

Prerequisites

pip install "rasterio>=1.3.0" "geopandas>=0.14" "shapely>=2.0" "numpy>=1.23" "pandas>=2.0"

Conceptually you need three things in place first. You need the label geometries and the imagery in a single CRS, which is the reprojection work described in mastering CRS transformations in rasterio. You need to know the scale factor and nodata value of the imagery, from extracting nodata and dtype from a GeoTIFF. And you need a cloud mask, because a chip that is half cloud teaches the model that clouds are whatever class the label says — the approaches are in cloud and shadow masking strategies.

---

Step-by-Step Workflow

Step 1 — Put labels and pixels in the same CRS

The first operation on any label layer is a reprojection into the raster’s CRS, followed by a validity check. Invalid geometries — self-intersecting polygons from a hand-digitised layer are the usual culprit — rasterize into surprising shapes rather than raising.

import geopandas as gpd
import rasterio

with rasterio.open("s2_stack_36NYF.tif") as src:
    raster_crs = src.crs
    raster_bounds = src.bounds

labels = gpd.read_file("field_labels.gpkg")
if labels.crs != raster_crs:
    labels = labels.to_crs(raster_crs)

# Shapely 2 make_valid fixes self-intersections without dropping features
labels["geometry"] = labels.geometry.make_valid()
labels = labels[labels.geometry.is_valid & ~labels.geometry.is_empty]

# Keep only labels that actually fall inside the scene
labels = labels.cx[raster_bounds.left:raster_bounds.right,
                   raster_bounds.bottom:raster_bounds.top]
print(f"{len(labels)} usable label features in {raster_crs}")

The .cx spatial slice is cheap and removes the single most wasteful failure mode in chip generation: iterating over a national label layer to build chips from one scene.

Step 2 — Choose chip size and stride from the objects, not the model

Chip size is a geography decision before it is a hyperparameter. The chip must be large enough that a typical object and enough of its surroundings fit inside it, because context is exactly what a convolutional model is for.

At Sentinel-2’s 10 m resolution a 256-pixel chip covers 2.56 km — comfortably more than a field, less than a catchment. For Landsat at 30 m the same pixel count covers 7.7 km, which is usually more context than a field-scale problem needs, so 128 is often the better choice there. The rule of thumb that holds across sensors: the chip should be roughly four to eight times the diameter of the objects you are labelling.

import math

import rasterio
from rasterio.windows import Window


def chip_windows(src: rasterio.DatasetReader, size: int = 256, stride: int | None = None):
    """Yield the windows of a regular chip grid over a scene."""
    stride = stride or size
    for row in range(0, src.height - 1, stride):
        for col in range(0, src.width - 1, stride):
            yield Window(col_off=col, row_off=row, width=size, height=size)


with rasterio.open("s2_stack_36NYF.tif") as src:
    total = sum(1 for _ in chip_windows(src, size=256, stride=128))
    gsd = abs(src.transform.a)
    print(f"{total} chips, each {256 * gsd:.0f} m across at {gsd:.0f} m GSD")

Step 3 — Read the image chip

The read itself is a windowed read with boundless=True, so a chip that runs off the edge of the scene is padded rather than truncated. A fixed chip shape is a hard requirement for any batched model, and discovering the exception at chip 40,000 of 50,000 is a bad way to spend an afternoon.

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


def read_image_chip(src: rasterio.DatasetReader, win: Window,
                    scale: float = 10_000.0) -> tuple[np.ndarray, object]:
    """Read one chip as float32 reflectance plus the transform that locates it."""
    fill = src.nodata if src.nodata is not None else 0
    arr = src.read(window=win, boundless=True, fill_value=fill)
    valid = arr[0] != fill
    chip = arr.astype("float32") / scale        # scaled ints -> reflectance
    chip[:, ~valid] = np.nan                    # NaN marks "no observation"
    return chip, src.window_transform(win)

Dividing by the scale factor here, once, at read time, is the discipline that prevents the double-scaling bug described in handling nodata and scale factors in band math. Record the scale in the dataset manifest so inference applies exactly the same transformation.

Step 4 — Burn the label mask onto the chip’s own grid

This is the step that makes or breaks alignment. rasterio.features.rasterize accepts a transform argument, and passing the chip’s transform — not the scene’s — is what guarantees the mask lines up pixel for pixel with the array just read.

import numpy as np
from rasterio.features import rasterize


def burn_label_mask(labels_gdf, transform, shape: tuple[int, int],
                    class_field: str = "class_id") -> np.ndarray:
    """Rasterize label polygons onto exactly the chip grid."""
    shapes = ((geom, int(value))
              for geom, value in zip(labels_gdf.geometry, labels_gdf[class_field])
              if geom is not None and not geom.is_empty)
    return rasterize(
        shapes=shapes,
        out_shape=shape,
        transform=transform,     # the CHIP transform, never the scene transform
        fill=0,                  # 0 = unlabelled, matching the prediction convention
        all_touched=False,       # only pixels whose centre falls inside the polygon
        dtype="uint8",
    )

all_touched is the parameter people get wrong. Leaving it False uses the pixel-centre rule, which matches how area is computed everywhere else in the stack and keeps thin polygons from inflating their class. Setting it True marks every pixel the polygon touches at all, which is right for linear features such as roads and rivers that would otherwise fall between pixel centres and vanish. The trade-off is the same one discussed in rasterizing vector polygons onto a raster grid.

Step 5 — Filter, tag, and persist

Not every chip is worth keeping. A chip that is 90% nodata, or 70% cloud, or contains no labelled pixel at all, costs training time and teaches nothing. Filter aggressively, and record why each surviving chip survived.

import numpy as np
import pandas as pd


def keep_chip(chip: np.ndarray, mask: np.ndarray, cloud: np.ndarray,
              *, min_valid: float = 0.7, min_labelled: int = 64) -> bool:
    valid_frac = float(np.isfinite(chip[0]).mean())
    cloud_frac = float(cloud.mean())
    labelled = int((mask > 0).sum())
    return (valid_frac >= min_valid
            and cloud_frac <= 0.2
            and labelled >= min_labelled)


def block_fold(transform, block_km: float = 20.0, n_folds: int = 5) -> int:
    """Fold id from the chip's map coordinates, so nearby chips share a fold."""
    x, y = transform.c, transform.f          # chip origin in projected metres
    bx, by = int(x // (block_km * 1000)), int(y // (block_km * 1000))
    return (bx * 31 + by * 17) % n_folds     # deterministic, spatially blocky

The block_fold function is small and does most of the work of keeping validation honest: chips that are physically close land in the same fold, so a held-out fold is a held-out region. A random np.random.randint in its place would produce the inflated scores shown on the section overview. The more careful variants — buffered blocks, leave-one-region-out — are in designing spatial cross-validation for raster models.

Write the chips into one array file per fold and keep a sidecar index table:

import numpy as np
import pandas as pd

images, masks, index = [], [], []
# ... loop appending chip arrays, mask arrays and metadata rows ...

np.savez_compressed("chips_fold0.npz",
                    images=np.stack(images), masks=np.stack(masks))
pd.DataFrame(index).to_parquet("chips_fold0_index.parquet")

---

Parameter Reference

---

Verification & Testing

Three assertions catch almost every alignment bug before a model ever sees the data.

import numpy as np


def verify_chip(chip: np.ndarray, mask: np.ndarray, transform) -> None:
    # 1. Image and mask share a grid
    assert chip.shape[-2:] == mask.shape, (chip.shape, mask.shape)

    # 2. The transform is square and north-up, as every chip read should be
    assert abs(transform.a) == abs(transform.e), "non-square pixels"
    assert transform.b == 0 and transform.d == 0, "rotated transform"

    # 3. Labelled pixels are observed pixels
    labelled_but_missing = (mask > 0) & ~np.isfinite(chip[0])
    assert labelled_but_missing.sum() == 0, "labels over nodata"

The third assertion is the useful one. A label sitting over a nodata pixel means either the label layer extends past the scene or the two are misaligned; both are worth failing the build over.

For a visual check, write a handful of chips back out as GeoTIFFs with their transforms and open them alongside the label layer in any desktop GIS. If the mask edges sit exactly on the field boundaries in the imagery, the pipeline is correct. If they sit consistently one field to the north-west, the label reprojection in step 1 was skipped.

A statistical check worth running once per dataset: the class histogram per fold. Wildly different class proportions between folds mean the blocks are too large relative to the scene, and the cross-validation will report variance that has nothing to do with the model.

---

Troubleshooting

ValueError: operands could not be broadcast together when combining chip and mask

The chip was read with boundless=False near the scene edge and came back smaller than out_shape. Pass boundless=True and a fill_value to src.read, and keep out_shape in rasterize pinned to the nominal chip size rather than chip.shape.

Every mask comes back empty

Either the labels are in a different CRS from the transform used to burn them, or the class_field values are strings. rasterize silently produces an empty array when no geometry intersects the window; add an assertion that at least one chip in the first hundred has a non-zero mask and fail fast.

TopologyException or invalid geometry during rasterization

A self-intersecting polygon survived the validity filter. Call .make_valid() on the whole layer, then drop empty geometries. On Shapely versions below 2.0, geom.buffer(0) is the equivalent fix.

Chips look correct but the model plateaus at the majority class

The dataset is dominated by background. Check the labelled-pixel fraction across chips; if 95% of pixels are class 0 the loss is optimised by predicting background everywhere. Raise min_labelled, sample chips centred on labels rather than on a regular grid, or weight the loss — the options are compared in balancing class frequencies in chip datasets.

Memory grows until the process is killed

Chips are being accumulated in a Python list before writing. Write in batches of a few thousand and clear the list, or stream directly into a chunked array format. A 50,000-chip dataset at 256×256×6 float32 is 75 GB — it was never going to fit.

---

Frequently Asked Questions

Q: What chip size should I use for Sentinel-2 segmentation? 256 by 256 at 10 m covers 2.56 km, which is large enough to contain whole fields and their surroundings and small enough to batch efficiently. Use 512 when the objects of interest are larger than about a kilometre, and 128 only when labels are sparse and you need more chips from the same scene.

Q: Should I store chips as individual GeoTIFFs or one array file? One array file per fold, in a chunked format, is faster to train from and far kinder to a filesystem than a hundred thousand small GeoTIFFs. Keep a sidecar table of chip transforms and source scene ids so any chip can still be traced back to its location.

Q: How do I stop chips from overlapping between training and validation? Assign folds by a spatial block id computed from the chip centre, not by a random shuffle, and make the block edge at least one chip wider than the chip itself so no training chip touches a validation chip.

Q: Can I build chips directly from a STAC search instead of a local file? Yes, and it is usually the better route for multi-scene datasets. Search the catalog, load the items into an array with the patterns in querying STAC catalogs programmatically, then take windows from that array exactly as you would from a local dataset.

---

Related

• Extracting Image Chips Around Labelled Points — centring windows on point labels rather than a regular grid.
• Rasterizing Labels into Segmentation Masks — multi-class burn order, priority rules, and boundary handling.
• Splitting Chips into Spatially Disjoint Folds — block assignment, buffers, and leave-one-region-out designs.
• Balancing Class Frequencies in Chip Datasets — sampling, weighting, and when to do neither.
• Rasterizing Vector Polygons onto a Raster Grid — the general-purpose burn used outside the modelling context.