Tiling Large Rasters into Training Patches

Cut a labelled scene into training patches without leaking: block-level splits, stride and overlap choices, label-density filtering, and a manifest that makes the dataset reproducible.

Assign spatial blocks to train, validation and test first, then cut patches inside each block, dropping any patch whose footprint crosses a block boundary. Filter on label density, write the patches to a materialised store, and record every patch’s source window in a manifest. That order — blocks, then patches — is what separates a dataset that measures generalisation from one that measures memorisation.

This page covers the patch pipeline. For what the patches are fed into, see Deep Learning Segmentation for Satellite Imagery.

Blocks decide the split; patches are cut inside them A scene is divided into four large blocks, two assigned to training and one each to validation and test. Inside the training blocks, overlapping patches are drawn. Patches whose footprint would straddle two blocks are shown crossed out, because they are discarded rather than assigned to either split. Split at block scale, sample at patch scale, discard the straddlers block 1 · train block 2 · train block 3 · validation block 4 · test patches inside a block stride = 0.75 × patch straddler → discarded costs a few % of data, buys an honest score

Why This Fails in Geospatial ML Pipelines

Patch-level random splitting is the default in every image-classification tutorial and it is wrong here. Two patches cut 64 pixels apart share three quarters of their pixels. Put one in train and one in validation and the validation score measures how well the model reproduces pixels it has already seen. In practice this inflates intersection-over-union by 20–30 points, which is enough to make a broken model look production-ready — the image-space version of the leakage described in reducing spatial leakage in model training.

The second failure is class starvation. Rare classes occupy a small fraction of a scene, so a naive grid of patches is overwhelmingly negative. Keep them all and most batches carry no gradient for the rare class; drop them all and the model has never seen what “not the class” looks like, so it hallucinates the class across the whole scene at inference. The fix is a controlled negative fraction, chosen and recorded rather than emergent.

Third, irreproducibility. A patch dataset cut by an ad-hoc script is a pile of arrays with no provenance: nobody can say which scene window a given patch came from, whether the fifth run used the same random seed, or why the class balance shifted between two experiments. A manifest — one row per patch with its source window, split, and label histogram — turns the dataset into an artifact that can be versioned alongside the model, as described in spatial dataset versioning with DVC and lakeFS.

Core Principles

  • Blocks first, patches second. The split boundary must be much coarser than the patch.
  • Discard straddlers. A patch that crosses a split boundary belongs to neither side.
  • Overlap within a split only. Stride around 0.75 of the patch is a good default.
  • Control the negative fraction explicitly. Do not let the grid decide it.
  • Write a manifest. Source window, split, class histogram, seed.
  • Materialise the patches. Re-reading GeoTIFF windows every epoch starves the GPU.

Production-Ready Code

from __future__ import annotations

import json
import logging
from dataclasses import dataclass, asdict

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

logger = logging.getLogger(__name__)
IGNORE = 255


@dataclass(frozen=True)
class PatchSpec:
    patch: int = 256
    stride: int = 192            # 0.75 overlap within a split
    block_px: int = 2048         # split granularity
    min_label_frac: float = 0.10 # at least this share of pixels labelled
    negative_fraction: float = 0.2
    seed: int = 0


def assign_blocks(height: int, width: int, spec: PatchSpec,
                  fractions=(0.7, 0.15, 0.15)) -> dict:
    """Assign whole blocks to train/val/test. Deterministic given spec.seed."""
    rows = int(np.ceil(height / spec.block_px))
    cols = int(np.ceil(width / spec.block_px))
    ids = np.arange(rows * cols)
    np.random.default_rng(spec.seed).shuffle(ids)

    n_tr = int(round(fractions[0] * len(ids)))
    n_va = int(round(fractions[1] * len(ids)))
    return {"train": set(ids[:n_tr].tolist()),
            "val": set(ids[n_tr:n_tr + n_va].tolist()),
            "test": set(ids[n_tr + n_va:].tolist()),
            "rows": rows, "cols": cols}


def _block_of(r: int, c: int, spec: PatchSpec, cols: int) -> int:
    return (r // spec.block_px) * cols + (c // spec.block_px)


def candidate_windows(height: int, width: int, spec: PatchSpec, blocks: dict):
    """Yield (row, col, split) for patches lying entirely within one block."""
    for r in range(0, height - spec.patch + 1, spec.stride):
        for c in range(0, width - spec.patch + 1, spec.stride):
            b0 = _block_of(r, c, spec, blocks["cols"])
            b1 = _block_of(r + spec.patch - 1, c + spec.patch - 1, spec, blocks["cols"])
            if b0 != b1:
                continue                                  # straddles a boundary
            for split in ("train", "val", "test"):
                if b0 in blocks[split]:
                    yield r, c, split
                    break


def build_patch_dataset(feature_path: str, label_path: str, out_prefix: str,
                        spec: PatchSpec = PatchSpec()) -> list[dict]:
    """Cut, filter and materialise patches; return the manifest rows."""
    with rasterio.open(feature_path) as fsrc, rasterio.open(label_path) as lsrc:
        if (fsrc.width, fsrc.height) != (lsrc.width, lsrc.height) or fsrc.crs != lsrc.crs:
            raise ValueError("feature and label rasters are not on the same grid")
        blocks = assign_blocks(fsrc.height, fsrc.width, spec)
        rng = np.random.default_rng(spec.seed + 1)

        manifest, kept_x, kept_y = [], [], []
        n_seen = n_empty = n_dropped = 0

        for r, c, split in candidate_windows(fsrc.height, fsrc.width, spec, blocks):
            n_seen += 1
            win = Window(c, r, spec.patch, spec.patch)
            y = lsrc.read(1, window=win)

            labelled = float((y != IGNORE).mean())
            if labelled < spec.min_label_frac:
                n_dropped += 1
                continue

            classes, counts = np.unique(y[y != IGNORE], return_counts=True)
            is_negative = len(classes) == 1 and classes[0] == 0
            if is_negative:
                n_empty += 1
                if rng.random() > spec.negative_fraction:
                    continue

            x = fsrc.read(window=win).astype("float32")
            kept_x.append(x)
            kept_y.append(y)
            manifest.append({
                "index": len(manifest), "split": split, "row": int(r), "col": int(c),
                "labelled_frac": round(labelled, 4), "negative": bool(is_negative),
                "classes": {int(k): int(v) for k, v in zip(classes, counts)},
            })

    np.save(f"{out_prefix}_x.npy", np.stack(kept_x))
    np.save(f"{out_prefix}_y.npy", np.stack(kept_y))
    with open(f"{out_prefix}_manifest.json", "w") as fh:
        json.dump({"spec": asdict(spec), "patches": manifest}, fh)

    logger.info("%d candidate windows -> %d patches (%d sparse dropped, %d negatives seen)",
                n_seen, len(manifest), n_dropped, n_empty)
    return manifest

Step-by-Step Walkthrough

Step 1 — Verify the grids before cutting anything

import rasterio

with rasterio.open("features.tif") as f, rasterio.open("labels.tif") as l:
    assert f.crs == l.crs and (f.width, f.height) == (l.width, l.height)
    assert f.transform.almost_equals(l.transform, precision=1e-6)

An offset here becomes an offset in every patch; the diagnosis is in aligning rasterized labels to an existing feature grid.

Step 2 — Cut, and read the summary line

spec = PatchSpec(patch=256, stride=192, block_px=2048,
                 min_label_frac=0.1, negative_fraction=0.2, seed=0)
manifest = build_patch_dataset("features.tif", "labels.tif", "data/patches", spec)

The log line is the health check. If “candidate windows” is far larger than “patches”, the label coverage is thin and the scene may not be worth tiling at this patch size.

Step 3 — Check the split-level class balance

import collections
import pandas as pd

rows = []
for p in manifest:
    for cls, n in p["classes"].items():
        rows.append({"split": p["split"], "class": int(cls), "pixels": n})
balance = pd.DataFrame(rows).groupby(["split", "class"])["pixels"].sum().unstack(fill_value=0)
print((balance.T / balance.sum(axis=1)).T.round(4))

The three splits should have similar class proportions. If the test split has three times the water of the training split, the block assignment happened to isolate a lake — regenerate with a different seed, or stratify block assignment on class presence.

Step 4 — Confirm no patch pair crosses a split

import itertools

by_split = collections.defaultdict(list)
for p in manifest:
    by_split[p["split"]].append((p["row"], p["col"]))

for a, b in itertools.combinations(by_split, 2):
    for (r1, c1) in by_split[a]:
        for (r2, c2) in by_split[b]:
            assert abs(r1 - r2) >= spec.patch or abs(c1 - c2) >= spec.patch, (
                f"patches from {a} and {b} overlap at ({r1},{c1}) and ({r2},{c2})")

Run this once on a small scene rather than every build — it is quadratic — but run it, because it is the only direct proof that the split is clean.

Keep some negatives, not all and not none A grid of candidate patches is shown. A small number contain the rare class and are all retained. The majority contain only background; a controlled fraction of these is retained and the rest discarded, producing a training set where the rare class is well represented but background is still learned. Candidate patches → retained training set 20 candidates · 3 contain the rare class sample 3 positives + 4 of 17 negatives negative_fraction = 0.2, recorded in the manifest Keep none and the model paints the class everywhere; keep all and most batches teach it nothing.

Verification

import numpy as np


def test_no_patch_straddles_a_block():
    spec = PatchSpec(patch=256, stride=192, block_px=2048, seed=0)
    blocks = assign_blocks(8192, 8192, spec)
    for r, c, _ in candidate_windows(8192, 8192, spec, blocks):
        assert (r // spec.block_px) == ((r + spec.patch - 1) // spec.block_px)
        assert (c // spec.block_px) == ((c + spec.patch - 1) // spec.block_px)


def test_same_seed_gives_the_same_dataset(tmp_path):
    a = build_patch_dataset("features.tif", "labels.tif", str(tmp_path / "a"), PatchSpec(seed=7))
    b = build_patch_dataset("features.tif", "labels.tif", str(tmp_path / "b"), PatchSpec(seed=7))
    assert [(p["row"], p["col"], p["split"]) for p in a] == \
           [(p["row"], p["col"], p["split"]) for p in b]


def test_negative_fraction_is_respected(tmp_path):
    spec = PatchSpec(negative_fraction=0.2, seed=3)
    man = build_patch_dataset("features.tif", "labels.tif", str(tmp_path / "c"), spec)
    neg = sum(p["negative"] for p in man)
    assert 0.05 <= neg / max(1, len(man)) <= 0.45


def test_sparse_patches_are_dropped(tmp_path):
    spec = PatchSpec(min_label_frac=0.5)
    man = build_patch_dataset("features.tif", "labels.tif", str(tmp_path / "d"), spec)
    assert all(p["labelled_frac"] >= 0.5 for p in man)

Beyond the tests, render a handful of patches with their labels overlaid before every training run. Five minutes of looking catches offsets, band-order mistakes and inverted masks that no assertion will.

The manifest is what makes a patch dataset an artifact A table shows one row per patch with columns for the patch index, the split it belongs to, its source row and column in the scene, the fraction of its pixels that are labelled, and whether it is a negative patch. This record lets any patch be traced back to its source window. One row per patch — traceable back to the scene indexsplitrow collabellednegative 0train2 3041 1520.86false 1train2 3041 3440.41true 2val6 1444 0960.93false Version this file with the arrays. Two runs with the same seed must produce the same manifest. A patch you cannot trace is a patch you cannot debug.

FAQ

Should training patches overlap?

Within a split, yes — a stride of about 0.75 of the patch adds augmentation and ensures every object appears somewhere other than at an edge. Across splits, never: overlapping patches in different splits are the leak.

How do I stop empty patches from dominating the dataset?

Filter on label density and keep a controlled fraction of negatives — 10–30% is a sensible start. Dropping them all teaches the model the rare class is everywhere; keeping them all starves it of gradient. The same trade-off appears in handling class imbalance in land cover classification.

Should patches be written to disk or cut on the fly?

Written. Reading compressed GeoTIFF windows every epoch usually makes the loader the bottleneck; a materialised store plus a manifest is faster and reproducible.

Part of: Deep Learning Segmentation for Satellite Imagery Part of: Training Geospatial Predictive Models in Python