⚙️ How-To Spatial machine learning Raster & rasterio Spatial Machine Learning · step 4 of 14 Updated

How to Sample Training Points from Rasters and Polygons

Problem statement

Training rows come from somewhere, and where they come from determines what the model learns. Three decisions matter more than the algorithm:

  • how many, which matters less than people expect
  • where, which matters more
  • how they are grouped, which determines whether validation means anything

Measured on a real problem, training on randomly placed cells and predicting the rest:

training cells    true map R²
           100         0.413
           300         0.480
         1,000         0.498
         3,000         0.570
        10,000         0.599

A hundredfold increase in rows bought 0.19 R². That flatness is the signature of spatially correlated data: neighbouring cells repeat information, so extra rows add far less than the count suggests.

Quick answer 

import numpy as np
import rasterio


def sample_raster(paths, n=2000, seed=0, mask_nodata=True):
    """Random valid cells, with their coordinates kept for grouping later."""
    with rasterio.open(paths[0]) as src:
        shape, transform = src.shape, src.transform

    stack, names = [], []
    for name, path in paths.items():
        with rasterio.open(path) as src:
            band = src.read(1).astype("float32")
            if mask_nodata and src.nodata is not None:
                band[band == src.nodata] = np.nan
        stack.append(band)
        names.append(name)

    data = np.stack(stack)
    valid = np.isfinite(data).all(axis=0)
    rows, cols = np.nonzero(valid)

    rng = np.random.default_rng(seed)
    pick = rng.choice(len(rows), size=min(n, len(rows)), replace=False)
    rows, cols = rows[pick], cols[pick]

    x = transform.c + (cols + 0.5) * transform.a
    y = transform.f + (rows + 0.5) * transform.e
    return data[:, rows, cols].T, np.column_stack([x, y]), names

Keeping the coordinates is not optional — every spatial validation scheme needs them.

Random, stratified and clustered training samples over a study area, with the clustered design leaving large unsampled gaps.
The same budget spent three ways. Clustering is what real surveys drift into.

Step-by-step solution

1. Sample only complete cases

A row with a missing feature is dropped by most models, and the dropping is not random — it follows wherever that layer has gaps. Filtering to cells where every layer is valid, before sampling, makes the sample representative of what the model can actually use.

2. Keep the coordinates

Without them there is no spatial cross-validation, no applicability mask, and no way to check whether the sample is clustered. Carry them alongside the feature matrix from the start.

3. Prefer stratified sampling for classification

Random sampling gives class proportions matching the landscape, which means rare classes get very few rows. A class covering 1% of the area gets 20 rows out of 2,000 — not enough to learn.

per_class = min_count = 200
for label in classes:
    idx = np.nonzero(labels == label)[0]
    pick = rng.choice(idx, min(per_class, len(idx)), replace=False)

Stratifying changes the prior, so the model's predicted probabilities no longer reflect landscape frequencies. That is a trade to make deliberately — see How to handle class imbalance in a spatial classifier.

4. Avoid sampling adjacent cells

Two neighbouring 20 m cells are nearly identical. Sampling both adds a row and almost no information, and it puts near-duplicates in different cross-validation folds.

Enforce a minimum separation, or sample on a coarse grid with a random offset within each cell.

5. Record the group structure

Field data usually arrives in plots, transects or survey days. Rows from one plot are not independent, and they must stay together in every fold.

Carry a group column from the start; GroupKFold needs it and there is no way to reconstruct it later.

True map R-squared rising only from 0.413 to 0.599 as training cells rise from 100 to 10,000.
A hundredfold increase in rows bought 0.19 R². Correlated rows are not independent observations.

Code examples

Example 1 — sampling with a minimum separation 

import numpy as np
import rasterio
from scipy.spatial import cKDTree


def sample_separated(layers, transform, n=2000, min_distance=None, seed=0):
    """Random valid cells, no two closer than min_distance."""
    data = np.stack(list(layers.values()))
    valid = np.isfinite(data).all(axis=0)
    rows, cols = np.nonzero(valid)

    x = transform.c + (cols + 0.5) * transform.a
    y = transform.f + (rows + 0.5) * transform.e
    coords = np.column_stack([x, y])

    rng = np.random.default_rng(seed)
    order = rng.permutation(len(coords))

    if min_distance is None:
        pick = order[:n]
    else:
        chosen, chosen_coords = [], []
        for i in order:
            if len(chosen) >= n:
                break
            if chosen_coords:
                tree = cKDTree(np.array(chosen_coords))
                if tree.query(coords[i], k=1)[0] < min_distance:
                    continue
            chosen.append(i)
            chosen_coords.append(coords[i])
        pick = np.array(chosen)
        print(f"  {len(pick):,} of {n:,} requested after enforcing "
              f"{min_distance} m separation")

    X = data[:, rows[pick], cols[pick]].T
    print(f"  {X.shape[0]:,} rows x {X.shape[1]} features from "
          f"{valid.sum():,} valid cells ({X.shape[0] / valid.sum():.2%})")
    return X, coords[pick], list(layers)

Rebuilding the tree each iteration is fine for a few thousand points and slow beyond that. For large samples, use a grid-based rejection instead: divide the area into cells of min_distance and take at most one point per cell.

Example 2 — stratified sampling from a class raster 

import numpy as np


def stratified_sample(features, labels, per_class=200, min_class=30, seed=0):
    """Equal-ish rows per class, reporting what each class could supply."""
    rng = np.random.default_rng(seed)
    classes, counts = np.unique(labels, return_counts=True)

    picked, dropped = [], []
    for label, available in zip(classes, counts):
        idx = np.nonzero(labels == label)[0]
        if available < min_class:
            dropped.append((label, available))
            continue
        take = min(per_class, available)
        picked.append(rng.choice(idx, take, replace=False))
        print(f"  class {label}: {available:8,} available, {take:5,} taken "
              f"({take / available:6.2%})")

    for label, available in dropped:
        print(f"  ! class {label} has only {available} cells — dropped "
              f"(minimum {min_class})")

    selection = np.concatenate(picked)
    print(f"  {len(selection):,} rows across {len(picked)} classes")
    print("  note: stratifying changes the class prior — predicted "
          "probabilities no longer reflect landscape frequencies")
    return features[selection], labels[selection], selection

The note about the prior is not decoration. A stratified classifier's probabilities are conditional on the sampling design, and using them as landscape probabilities is a common and consequential mistake.

Example 3 — sampling within polygons 

import numpy as np
import geopandas as gpd
from shapely.geometry import Point


def sample_in_polygons(gdf, per_polygon=10, label_column=None, seed=0,
                       max_attempts=100):
    """Random points inside each polygon, carrying the polygon id as a group."""
    rng = np.random.default_rng(seed)
    records = []

    for idx, row in gdf.iterrows():
        geometry = row.geometry
        if geometry is None or geometry.is_empty:
            continue
        west, south, east, north = geometry.bounds
        found, attempts = 0, 0
        while found < per_polygon and attempts < per_polygon * max_attempts:
            attempts += 1
            point = Point(rng.uniform(west, east), rng.uniform(south, north))
            if geometry.contains(point):
                records.append({
                    "x": point.x, "y": point.y, "group": idx,
                    "label": row[label_column] if label_column else None,
                })
                found += 1
        if found < per_polygon:
            print(f"  polygon {idx}: only {found}/{per_polygon} points "
                  f"({attempts} attempts) — very thin or narrow geometry")

    out = gpd.GeoDataFrame(records, crs=gdf.crs,
                           geometry=[Point(r["x"], r["y"]) for r in records])
    print(f"  {len(out):,} points from {out['group'].nunique()} polygons")
    return out

The group column is the important output. Ten points from one polygon are ten views of one place, and they must stay in the same cross-validation fold or the score is meaningless.

Rejection sampling struggles with long thin polygons — rivers, roads, field margins — where the bounding box is mostly outside the geometry. For those, sample along the geometry rather than within its box.

Explanation

Why more rows help so little

Neighbouring cells are strongly correlated: similar slope, similar elevation, similar target. The tenth cell in a neighbourhood adds a fraction of what the first added.

The measured curve makes this concrete: 100 to 10,000 training cells improved the true map R² from 0.413 to 0.599. A hundredfold increase in rows for a 0.19 gain.

Two consequences. Sampling more densely in the same places is nearly free of benefit, and any confidence interval computed from the row count is far too narrow — the effective sample size is closer to the number of distinct neighbourhoods.

Why the group column cannot be reconstructed

Rows from one field plot, one survey day or one polygon are not independent. If they land in different folds, the validation is measuring memorisation.

That grouping is knowable only at sampling time. Once the feature matrix is built, two rows from the same plot look like two rows from anywhere.

Carrying a group column from the start costs one column and makes correct validation possible. Adding it later usually means re-running the sampling.

Why stratified sampling changes the prior

A classifier learns the class frequencies in its training data. Sampling 200 rows from each of five classes teaches it that all five are equally common, whichever way the landscape is distributed.

The learned decision boundaries are usually improved by this — rare classes get enough examples to be learnable. The predicted probabilities are not landscape probabilities any more.

Where calibrated probabilities matter, either weight by class frequency at fit time or correct the priors afterwards. Where only the predicted class matters, stratifying is straightforward.

Why minimum separation is worth enforcing

Two adjacent cells are near-duplicates. Under random k-fold they land in different folds and the model predicts one from the other, which is memorisation scored as skill.

Enforcing a minimum separation of roughly the autocorrelation range makes the rows closer to independent. It also reduces the achievable sample size, sometimes drastically — which is honest, because those rows were never independent observations.

A group column carried from sampling through to validation, without which rows from one plot split across folds.
Adding it later usually means re-running the sampling. It costs one column now.

Edge cases or notes

  • Sample only complete cases, or incomplete rows are dropped non-randomly.
  • Keep coordinates — every spatial validation scheme needs them.
  • Carry a group column from sampling; it cannot be reconstructed later.
  • Enforce a minimum separation of roughly the autocorrelation range.
  • Stratifying changes the class prior, so probabilities need correcting.
  • Rejection sampling fails on long thin polygons. Sample along them instead.
  • More rows help far less than expected — 100× for 0.19 R² here.
  • Record the sampling design with the model; it is part of what the model means.

FAQ

How many training points do I need?

Fewer than you think, and spread further. Going from 100 to 10,000 cells improved the true map R² from 0.413 to 0.599 — a hundredfold increase for 0.19.

Should I sample randomly or on a grid?

A grid or stratified-random design gives better coverage for the same budget. Purely random sampling leaves gaps and clusters by chance.

Why keep the coordinates?

Spatial cross-validation, applicability masks and clustering diagnostics all need them. Without coordinates the sample cannot be validated spatially.

What is the group column for?

Rows from one plot, transect or polygon are not independent and must stay in the same fold. That structure is knowable only at sampling time.

Should I stratify by class?

For classification, usually yes — rare classes otherwise get too few rows. It changes the class prior, so predicted probabilities need correcting.

How far apart should samples be?

Roughly the autocorrelation range of the target. Closer than that and the rows are near-duplicates that leak across folds.

Why does my sample have fewer rows than I asked for?

Either the minimum separation could not be satisfied, or complete-case filtering removed cells where a layer had gaps.