How to Sample Training Points from Rasters and Polygons
On this page
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
The examples below use a real problem: predicting summer NDVI from terrain over a 10 km window of Snowdonia, from a Sentinel-2 scene and the Copernicus 30 m elevation model on a 20 m grid. Trained on randomly placed cells and scored on every other cell, a random forest's R² went from 0.194 with 100 training cells to 0.335 with 10,000 (Example 4). A hundredfold increase in rows bought 0.14 R², and the last 7,000 rows bought 0.036. That flattening is the signature of spatially correlated data: neighbouring cells repeat information, so extra rows add far less than the count suggests.
Quick answer
def sample_raster(paths, n=2000, seed=0, mask_nodata=True):
"""Random cells valid in every layer, with their coordinates kept for grouping later."""
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 and not np.isnan(src.nodata):
band[band == src.nodata] = np.nan
transform = src.transform
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
features, xy, names = sample_raster(LAYERS, n=2000)
print(features.shape, names)
print("first row:", {k: round(float(v), 3) for k, v in zip(names, features[0])}, "at", xy[0].round(0))(2000, 4) ['elevation', 'slope', 'tpi', 'ndvi']
first row: {'elevation': 290.461, 'slope': 12.78, 'tpi': -0.861, 'ndvi': 0.78} at [ 425590. 5883830.]
Keeping the coordinates is not optional — every spatial validation scheme needs them.
Complete cases, coordinates, strata, separation, groups
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:
for name, path in LAYERS.items():
with rasterio.open(path) as src:
band = src.read(1)
print(f"{name:10} missing in {np.isnan(band).mean():6.1%} of cells")
with rasterio.open(LAYERS["ndvi"]) as src:
ndvi_band = src.read(1)
missing = np.isnan(ndvi_band)
print(f"mean elevation where NDVI is missing {np.nanmean(elevation[missing]):.0f} m, "
f"where it is present {np.nanmean(elevation[~missing]):.0f} m")elevation missing in 0.0% of cells
slope missing in 0.0% of cells
tpi missing in 0.0% of cells
ndvi missing in 2.1% of cells
mean elevation where NDVI is missing 188 m, where it is present 408 m
Only 2.1% of cells lack an NDVI value — the ones the scene classification marked as water, cloud or shadow — but they are not scattered at random: they sit 220 m lower on average, because most of them are the lakes in the valley floors. A model trained on the complete cases has seen none of that ground, and knowing that is part of knowing what it can predict.
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. sample_raster returns 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:
classes, counts = np.unique(landcover[np.isfinite(landcover)], return_counts=True)
random_rows = np.random.default_rng(0).choice(landcover[np.isfinite(landcover)], 2000, replace=False)
print(f"{'class':22} {'share of area':>13} {'rows in a random 2,000':>23}")
for value, count in zip(classes, counts):
print(f"{WORLDCOVER[int(value)]:22} {count / counts.sum():13.2%} {int((random_rows == value).sum()):>23}")class share of area rows in a random 2,000
tree cover 11.81% 242
grassland 83.88% 1675
built-up 0.88% 16
bare / sparse 0.56% 12
water 2.86% 55
herbaceous wetland 0.00% 0
moss and lichen 0.00% 0
A class covering a fraction of a percent of the window gets a handful of rows out of 2,000 — not enough to learn. Stratifying fixes that and 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. Measure it on the target itself:
for lag in (1, 5, 25, 100): # cells: 20 m, 100 m, 500 m, 2 km
a, b = ndvi_band[:, :-lag], ndvi_band[:, lag:]
both = np.isfinite(a) & np.isfinite(b)
print(f"NDVI correlation with the cell {lag * 20:>5} m to the east: {np.corrcoef(a[both], b[both])[0, 1]:.2f}")NDVI correlation with the cell 20 m to the east: 0.91
NDVI correlation with the cell 100 m to the east: 0.67
NDVI correlation with the cell 500 m to the east: 0.38
NDVI correlation with the cell 2000 m to the east: 0.12
Sampling both of two adjacent cells 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. Example 3 carries the polygon a point came from as its group.
Separated, stratified and in-polygon samples over Snowdonia
The data used on this page
A Sentinel-2 L2A scene of Snowdonia from summer 2025 (the least cloudy on Element 84's Earth Search catalogue between June and September; Copernicus Sentinel data), the Copernicus GLO-30 elevation model, ESA WorldCover 2021 (CC BY 4.0), all on the scene's 20 m grid; and HM Land Registry's INSPIRE parcels for Rutland (© Crown copyright and database rights; OGL) for sampling inside polygons. The setup writes each layer to a GeoTIFF, which is what a sampling script usually starts from.
from pathlib import Path
import geopandas as gpd
import numpy as np
import rasterio
import rasterio.merge
import requests
from pystac_client import Client
from rasterio.enums import Resampling
from rasterio.warp import reproject, transform_bounds
from scipy import ndimage
from scipy.spatial import cKDTree
from shapely.geometry import Point
Path("data").mkdir(exist_ok=True)
AOI = (-4.16, 53.03, -4.00, 53.12)
WORLDCOVER = {10: "tree cover", 20: "shrubland", 30: "grassland", 40: "cropland", 50: "built-up",
60: "bare / sparse", 70: "snow and ice", 80: "water", 90: "herbaceous wetland",
95: "mangroves", 100: "moss and lichen"}
LAYERS = {name: f"data/{name}.tif" for name in ("elevation", "slope", "tpi", "ndvi")}
if not Path("data/ndvi.tif").exists():
items = Client.open("https://earth-search.aws.element84.com/v1").search(
collections=["sentinel-2-l2a"], bbox=list(AOI), datetime="2025-06-01/2025-10-01",
query={"eo:cloud_cover": {"lt": 60}}).items()
item = min(items, key=lambda i: i.properties["eo:cloud_cover"])
with rasterio.open(item.assets["scl"].href) as ds: # the 20 m grid
win = rasterio.windows.from_bounds(*transform_bounds("EPSG:4326", ds.crs, *AOI),
transform=ds.transform).round_offsets().round_lengths()
scl, transform, crs = ds.read(1, window=win), ds.window_transform(win), ds.crs
bounds = rasterio.windows.bounds(win, ds.transform)
shape = scl.shape
band = {}
for name in ("red", "nir"):
with rasterio.open(item.assets[name].href) as ds:
band[name] = ds.read(1, window=rasterio.windows.from_bounds(*bounds, transform=ds.transform),
out_shape=shape, resampling=Resampling.average).astype("float32") / 1e4
ndvi = (band["nir"] - band["red"]) / (band["nir"] + band["red"])
ndvi[~np.isin(scl, (2, 4, 5, 7))] = np.nan # keep clear land only
tiles = [rasterio.open(f"https://copernicus-dem-30m.s3.amazonaws.com/Copernicus_DSM_COG_10_N53_00_{t}_00_DEM/"
f"Copernicus_DSM_COG_10_N53_00_{t}_00_DEM.tif") for t in ("W005", "W004")]
mosaic, mosaic_transform = rasterio.merge.merge(tiles, bounds=(-4.18, 53.02, -3.98, 53.13))
elev = np.full(shape, np.nan, dtype="float32")
reproject(mosaic[0], elev, src_transform=mosaic_transform, src_crs=tiles[0].crs,
dst_transform=transform, dst_crs=crs, resampling=Resampling.bilinear)
dz_dy, dz_dx = np.gradient(elev, -20, 20)
layers = {"elevation": elev, "slope": np.degrees(np.arctan(np.hypot(dz_dx, dz_dy))).astype("float32"),
"tpi": (elev - ndimage.uniform_filter(elev, size=15)).astype("float32"), "ndvi": ndvi}
cover = np.zeros(shape, dtype="uint8")
with rasterio.open("https://esa-worldcover.s3.eu-central-1.amazonaws.com/v200/2021/map/"
"ESA_WorldCover_10m_2021_v200_N51W006_Map.tif") as ds:
reproject(rasterio.band(ds, 1), cover, dst_transform=transform, dst_crs=crs, resampling=Resampling.mode)
layers["landcover"] = np.where(cover > 0, cover, np.nan).astype("float32")
profile = dict(driver="GTiff", height=shape[0], width=shape[1], count=1, dtype="float32",
crs=crs, transform=transform, nodata=np.nan, compress="deflate")
for name, values in layers.items():
with rasterio.open(f"data/{name}.tif", "w", **profile) as out:
out.write(values.astype("float32"), 1)
Path("data/scene.txt").write_text(item.id)
if not Path("data/Rutland.zip").exists():
Path("data/Rutland.zip").write_bytes(requests.get(
"https://use-land-property-data.service.gov.uk/datasets/inspire/download/Rutland.zip", timeout=300).content)
with rasterio.open("data/elevation.tif") as src:
elevation, grid = src.read(1), src.transform
with rasterio.open("data/landcover.tif") as src:
landcover = src.read(1)
print(Path("data/scene.txt").read_text(), "|", elevation.shape[1], "x", elevation.shape[0], "cells of",
grid.a, "m | layers:", list(LAYERS) + ["landcover"])S2C_30UVD_20250712_0_L2A | 543 x 509 cells of 20.0 m | layers: ['elevation', 'slope', 'tpi', 'ndvi', 'landcover']
Example 1 — sampling with a minimum separation
def sample_separated(layers, n=2000, min_distance=None, seed=0):
"""Random valid cells, no two closer than min_distance (grid-based rejection)."""
data = np.stack([rasterio.open(path).read(1) for path in layers.values()])
valid = np.isfinite(data).all(axis=0)
rows, cols = np.nonzero(valid)
x = grid.c + (cols + 0.5) * grid.a
y = grid.f + (rows + 0.5) * grid.e
coords = np.column_stack([x, y])
order = np.random.default_rng(seed).permutation(len(coords))
if min_distance is None:
pick = order[:n]
else:
chosen, taken = [], {}
for i in order: # one bucket per min_distance square
if len(chosen) >= n:
break
key = (int(coords[i, 0] // min_distance), int(coords[i, 1] // min_distance))
near = [taken[(key[0] + dx, key[1] + dy)] for dx in (-1, 0, 1) for dy in (-1, 0, 1)
if (key[0] + dx, key[1] + dy) in taken]
if any(np.hypot(*(coords[i] - coords[j])) < min_distance for j in near):
continue
if key in taken:
continue
taken[key] = i
chosen.append(i)
pick = np.array(chosen)
X = data[:, rows[pick], cols[pick]].T
print(f" min distance {min_distance or 0:>4} m: {len(pick):,} of {n:,} requested "
f"({len(pick) / valid.sum():.2%} of {valid.sum():,} valid cells)")
return X, coords[pick]
for distance in (None, 200, 500, 1000):
X_sep, xy_sep = sample_separated(LAYERS, n=2000, min_distance=distance)
nearest = cKDTree(xy_sep).query(xy_sep, k=2)[0][:, 1]
print(f"at 1,000 m the closest pair is {nearest.min():,.0f} m apart") min distance 0 m: 2,000 of 2,000 requested (0.74% of 270,536 valid cells)
min distance 200 m: 1,908 of 2,000 requested (0.71% of 270,536 valid cells)
min distance 500 m: 322 of 2,000 requested (0.12% of 270,536 valid cells)
min distance 1000 m: 84 of 2,000 requested (0.03% of 270,536 valid cells)
at 1,000 m the closest pair is 1,001 m apart
Taking at most one point per min_distance square and checking the neighbouring squares keeps this fast; rebuilding a KD-tree after every accepted point is fine for a few hundred points and slow beyond that. The answer is honest about what the window holds: at 500 m separation it cannot supply 2,000 points, because it does not contain 2,000 places that far apart.
Example 2 — stratified sampling from a class raster
def stratified_sample(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)
flat = labels.ravel()
ok = np.isfinite(flat)
classes, counts = np.unique(flat[ok], return_counts=True)
picked, dropped = [], []
for label, available in zip(classes, counts):
if available < min_class:
dropped.append((label, available))
continue
idx = np.nonzero(flat == label)[0]
take = min(per_class, available)
picked.append(rng.choice(idx, take, replace=False))
print(f" {WORLDCOVER[int(label)]:20} {available:8,} available, {take:4,} taken ({take / available:7.2%})")
for label, available in dropped:
print(f" ! {WORLDCOVER[int(label)]} has only {available} cells — dropped (minimum {min_class})")
selection = np.concatenate(picked)
print(f" {len(selection):,} rows across {len(picked)} classes")
return selection
selection = stratified_sample(landcover) tree cover 32,644 available, 200 taken ( 0.61%)
grassland 231,845 available, 200 taken ( 0.09%)
built-up 2,441 available, 200 taken ( 8.19%)
bare / sparse 1,548 available, 200 taken ( 12.92%)
water 7,901 available, 200 taken ( 2.53%)
! herbaceous wetland has only 5 cells — dropped (minimum 30)
! moss and lichen has only 3 cells — dropped (minimum 30)
1,000 rows across 5 classes
A stratified classifier's probabilities are conditional on the sampling design: this sample says every class is equally common, which the landscape table above says is false. Using those probabilities as landscape probabilities is a common and consequential mistake.
Example 3 — sampling within polygons
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, short = [], []
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({"group": idx, "label": row[label_column] if label_column else None, "geometry": point})
found += 1
if found < per_polygon:
short.append((idx, found, geometry.area / (east - west) / (north - south)))
out = gpd.GeoDataFrame(records, crs=gdf.crs)
print(f" {len(out):,} points from {out['group'].nunique()} polygons; {len(short)} polygons came up short")
for idx, found, fill in short[:3]:
print(f" parcel {idx}: {found}/{per_polygon} points, fills {fill:.1%} of its bounding box")
return out
parcels = gpd.read_file("/vsizip/data/Rutland.zip/Land_Registry_Cadastral_Parcels.gml",
columns=["INSPIREID"], rows=2000, DOWNLOAD_SCHEMA="NO")
points = sample_in_polygons(parcels, per_polygon=10) 19,990 points from 1999 polygons; 1 polygons came up short
parcel 518: 0/10 points, fills 0.0% of its bounding box
The group column is the important output. Ten points from one parcel are ten views of one place, and they must stay in the same cross-validation fold or the score is meaningless.
The one parcel that came up short is the case rejection sampling struggles with: a sliver so thin that it fills almost none of its bounding box, so a hundred random tries per point never landed inside it. Long thin strips — roads, verges, field margins — behave the same way on a smaller scale. For those, sample along the geometry rather than within its box, or use geometry.sample_points() in GeoPandas 1.x, which does the rejection for you.
Example 4 — how much rows are worth, measured
from sklearn.ensemble import RandomForestRegressor
from sklearn.metrics import r2_score
stack = np.stack([rasterio.open(path).read(1) for path in LAYERS.values()])
ok = np.isfinite(stack).all(axis=0)
X_all, y_all = stack[:3, ok].T, stack[3, ok] # elevation, slope, TPI -> NDVI
forest = lambda: RandomForestRegressor(n_estimators=200, min_samples_leaf=2, n_jobs=-1, random_state=0)
print("training cells -> R2 over every other cell")
for n in (100, 300, 1000, 3000, 10000):
pick = np.random.default_rng(0).choice(len(y_all), n, replace=False)
others = np.setdiff1d(np.arange(len(y_all)), pick)
print(f" {n:6,} {r2_score(y_all[others], forest().fit(X_all[pick], y_all[pick]).predict(X_all[others])):+.3f}")training cells -> R2 over every other cell
100 +0.194
300 +0.219
1,000 +0.292
3,000 +0.299
10,000 +0.335
Why more rows help so little
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 correlation table in step 4 shows how slowly it decays.
The measured curve makes this concrete. 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 class teaches it that all of them 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 — Example 1 could not place 2,000 points 500 m apart in this window — which is honest, because those rows were never independent observations.
Class priors, thin polygons and recording the design
- 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.14 R² here.
- Record the sampling design with the model; it is part of what the model means.
Internal links
- Spatial machine learning explained — why rows are not independent
- How to build spatial features for a machine learning model — what to sample
- How to do spatial cross-validation in Python — what the groups are for
- How to handle class imbalance in a spatial classifier — stratification and its consequences
- Sample design explained: where to measure — the same question for interpolation
- How to extract raster values at point locations with rasterio — the sampling mechanics
- Spatial leakage explained: why random cross-validation lies — what adjacent samples do to a score
- How to find and remove spatial outliers in a point dataset — cleaning the sample
FAQ
How many training points do I need?
Fewer than you think, and spread further. Going from 100 to 10,000 cells improved the R² over every other cell from 0.194 to 0.335 — a hundredfold increase for 0.14.
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.