How to Build Spatial Features for a Machine Learning Model

Problem statement

Feature construction decides more about a spatial model than the algorithm does. Measured on a real problem — Sentinel-2 greenness predicted from terrain in Snowdonia — the same random forest with two extra columns:

terrain only (slope, elevation, northness, eastness, TPI)   R² 0.348
the same plus raw x and y                                   R² 0.509

Two columns, a large improvement in-area. And the same two columns gained far less under block cross-validation — 0.258 to 0.283 — because a coordinate learned here does not describe anywhere else.

The aim is to build features that carry the same information and travel: distances, neighbourhood summaries and terrain derivatives. Whether they do is a measurement, not an assumption — and on this page's own data, trained on the western half of the area and tested on the eastern, none of the terrain features did (Example 4).

Quick answer

import numpy as np
from scipy.spatial import cKDTree
from scipy.ndimage import uniform_filter


def spatial_features(coords, dem, transform, cell, reference_layers,
                     windows=(3, 9, 27)):
    """Terrain derivatives, neighbourhood summaries and distances."""
    features, names = [], []

    gy, gx = np.gradient(dem, cell)
    for name, array in (("elevation", dem),
                        ("slope", np.degrees(np.arctan(np.hypot(gx, gy)))),
                        ("northness", np.cos(np.arctan2(-gx, gy))),
                        ("eastness", np.sin(np.arctan2(-gx, gy)))):
        features.append(sample(array, coords, transform))
        names.append(name)

    for window in windows:
        smoothed = uniform_filter(np.nan_to_num(dem), size=window)
        features.append(sample(dem - smoothed, coords, transform))
        names.append(f"tpi_{int(window * cell)}m")

    for name, targets in reference_layers.items():
        distances, _ = cKDTree(targets).query(coords, k=1)
        features.append(distances)
        names.append(f"dist_{name}")

    return np.column_stack(features), names

median slope: 17.1° in UTM, 10.4° in Web Mercator; cos(53°) = 0.60

Raster sampling, terrain in metres, scales, relative features

1. Sample rasters correctly

def sample(array, coords, transform):
    """The raster value under each coordinate, in the raster's own CRS."""
    h, w = array.shape
    rows, cols = rasterio.transform.rowcol(transform, coords[:, 0], coords[:, 1])
    return array[np.clip(rows, 0, h - 1), np.clip(cols, 0, w - 1)]

Two details matter. rowcol expects the raster's own CRS, so the coordinates must already be in it. And clipping the indices handles points exactly on the far edge, which would otherwise index one past the end.

Where a point falls on nodata, the feature should be NaN rather than the fill value — a fill of −9999 sampled as an elevation will dominate every split in the model.

2. Compute terrain derivatives in a metric CRS

Slope is a ratio of vertical to horizontal distance. In a geographic CRS the horizontal distance is in degrees, so the slope is wrong by roughly the number of metres in a degree.

The same applies to any Web Mercator raster: distances there are inflated by 1/cos(latitude), so a slope computed from one is understated by that factor. Measured on this page's DEM, at 53° N:

from rasterio.warp import calculate_default_transform

mercator_transform, mw, mh = calculate_default_transform(crs, "EPSG:3857", shape[1], shape[0],
                                                         *rasterio.transform.array_bounds(*shape, transform))
mercator = np.full((mh, mw), np.nan, dtype="float32")
reproject(elev, mercator, src_transform=transform, src_crs=crs, src_nodata=np.nan,
          dst_transform=mercator_transform, dst_crs="EPSG:3857", resampling=Resampling.bilinear)
my, mx = np.gradient(mercator, -mercator_transform.e, mercator_transform.a)
mercator_slope = np.degrees(np.arctan(np.hypot(mx, my)))
print(f"median slope: {np.nanmedian(slope):.1f}° in UTM, {np.nanmedian(mercator_slope):.1f}° in Web Mercator; "
      f"cos(53°) = {np.cos(np.radians(53)):.2f}")

median slope: 17.1° in UTM, 10.4° in Web Mercator; cos(53°) = 0.60

3. Build neighbourhood features at several scales

for window in (3, 9, 27):
    smoothed = uniform_filter(data, size=window)

At 20 m cells those windows are 60 m, 180 m and 540 m — local, hillslope and catchment scale. They measure different things, and which one matters is usually the question the model can answer.

Two or three scales spanning an order of magnitude is enough. Beyond that the columns are strongly correlated and add little.

4. Prefer relative features

tpi = elevation - neighbourhood_mean

"Higher or lower than around here" transfers between regions; "742 metres" does not. Topographic position index, local relative slope and spectral values relative to a local mean are all more transferable than their absolute versions.

5. Handle NaN deliberately

Neighbourhood filters propagate NaN across the whole window, so one nodata cell can void a 27-cell neighbourhood. Filling before filtering avoids that and biases the result towards the fill value.

The defensible approach is to fill with the local mean, compute, and then mask the result back to the original valid area.

Terrain derivatives and vector distances that travel

The data used on this page

The problem from spatial machine learning explained: Sentinel-2 greenness (NDVI) over part of Snowdonia on a 20 m grid, the Copernicus 30 m DEM resampled onto it, and — new here — OpenStreetMap's rivers, streams and roads for the same area as vector layers.

import numpy as np
import rasterio
import rasterio.merge
from pystac_client import Client
from rasterio.enums import Resampling
from rasterio.warp import reproject, transform_bounds
from scipy import ndimage

AOI = (-4.16, 53.03, -4.00, 53.12)
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"):                                         # 10 m bands averaged onto it
    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"])

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)                           # rows run south
slope = np.degrees(np.arctan(np.hypot(dz_dx, dz_dy)))
aspect = np.arctan2(-dz_dx, -dz_dy)                                 # direction the slope faces, from north
tpi = elev - ndimage.uniform_filter(elev, size=15)                  # against the 300 m neighbourhood

valid = np.isin(scl, (2, 4, 5, 7)) & np.isfinite(ndvi) & np.isfinite(elev)
rows, cols = np.nonzero(valid)
xs, ys = rasterio.transform.xy(transform, rows, cols)
coords = np.column_stack([xs, ys])
X = np.column_stack([elev[valid], slope[valid], np.cos(aspect[valid]), np.sin(aspect[valid]), tpi[valid]])
y = ndvi[valid]
FEATURES = ["elevation", "slope", "northness", "eastness", "tpi"]
west = coords[:, 0] < np.median(coords[:, 0])
pick = np.random.default_rng(0).choice(np.nonzero(west)[0], 3000, replace=False)
X_train, y_train, X_new, y_new = X[pick], y[pick], X[~west], y[~west]
feature_names = FEATURES
print(f"{item.id}: {shape[1]} x {shape[0]} cells of 20 m, {valid.sum():,} usable "
      f"({valid.mean():.1%}); {len(X_train):,} training cells in the west, {len(X_new):,} to predict in the east")

# the DEM on the 20 m grid, as a file, and OSM rivers and roads for the same area
with rasterio.open("elev.tif", "w", driver="GTiff", height=shape[0], width=shape[1], count=1,
                   dtype="float32", crs=crs, transform=transform, nodata=np.nan) as out:
    out.write(elev, 1)

import json
import os

import geopandas as gpd
import requests
from shapely.geometry import LineString

if not os.path.exists("snowdonia_osm.json"):
    QUERY = """[out:json][timeout:120];
(way["waterway"~"^(river|stream)$"](53.03,-4.16,53.12,-4.00);
 way["highway"~"^(primary|secondary|tertiary|unclassified|residential|track)$"](53.03,-4.16,53.12,-4.00););
out geom tags;"""
    reply = requests.post("https://overpass-api.de/api/interpreter", data={"data": QUERY}, timeout=180,
                          headers={"User-Agent": "spatialworkflow.io guide"})
    reply.raise_for_status()
    json.dump(reply.json(), open("snowdonia_osm.json", "w"))
osm = gpd.GeoDataFrame(
    [{"kind": "river" if "waterway" in e["tags"] else "road",
      "geometry": LineString([(p["lon"], p["lat"]) for p in e["geometry"]])}
     for e in json.load(open("snowdonia_osm.json"))["elements"]], crs=4326).to_crs(crs)
rivers, roads = osm[osm["kind"] == "river"], osm[osm["kind"] == "road"]
print(f"{len(rivers):,} river and stream ways, {len(roads):,} road ways")

S2C_30UVD_20250712_0_L2A: 543 x 509 cells of 20 m, 270,536 usable (97.9%); 3,000 training cells in the west, 135,292 to predict in the east
1,157 river and stream ways, 361 road ways

Example 1 — terrain derivatives with the CRS checked

import numpy as np
import rasterio
from scipy.ndimage import uniform_filter


def terrain_features(dem_path, windows=(3, 9, 27)):
    """Slope, aspect components and multi-scale TPI, as aligned arrays."""
    with rasterio.open(dem_path) as src:
        if not src.crs.is_projected:
            raise ValueError(f"{dem_path} is in {src.crs}; slope needs a "
                             "projected CRS")
        if src.crs.to_epsg() == 3857:
            print("  ! Web Mercator: slope will be understated by "
                  "1/cos(latitude)")
        dem = src.read(1).astype("float32")
        if src.nodata is not None:
            dem[dem == src.nodata] = np.nan
        cell = abs(src.transform.a)
        transform, crs = src.transform, src.crs

    gy, gx = np.gradient(dem, cell)
    aspect = np.arctan2(-gx, gy)

    layers = {
        "elevation": dem,
        "slope": np.degrees(np.arctan(np.hypot(gx, gy))),
        "northness": np.cos(aspect),
        "eastness": np.sin(aspect),
    }

    valid = np.isfinite(dem)
    filled = np.where(valid, dem, np.nanmean(dem))
    for window in windows:
        mean = uniform_filter(filled, size=window)
        squared = uniform_filter(filled ** 2, size=window)
        layers[f"tpi_{int(window * cell)}m"] = np.where(valid, dem - mean, np.nan)
        layers[f"roughness_{int(window * cell)}m"] = np.where(
            valid, np.sqrt(np.maximum(squared - mean ** 2, 0)), np.nan)

    for name, array in layers.items():
        print(f"  {name:18} {np.nanmin(array):9.2f} .. {np.nanmax(array):9.2f}  "
              f"{np.isnan(array).mean():5.2%} nan")
    return layers, transform, crs, cell

layers, dem_transform, dem_crs, cell = terrain_features("elev.tif")

  elevation              53.68 ..   1073.06  0.00% nan
  slope                   0.00 ..     67.00  0.00% nan
  northness              -1.00 ..      1.00  0.00% nan
  eastness               -1.00 ..      1.00  0.00% nan
  tpi_60m               -11.01 ..      9.99  0.00% nan
  roughness_60m           0.00 ..     37.65  0.00% nan
  tpi_180m              -34.99 ..     44.60  0.00% nan
  roughness_180m          0.00 ..     81.58  0.00% nan
  tpi_540m             -105.57 ..    127.03  0.00% nan
  roughness_540m          1.31 ..    136.54  0.00% nan

Raising on a geographic CRS and warning on Web Mercator catches the two errors that make every terrain derivative wrong, and neither raises an exception on its own.

Example 2 — sampling all the layers at once

import numpy as np
import rasterio


def sample_layers(layers, coords, transform, mask_nodata=True):
    """One row per coordinate, one column per layer."""
    height, width = next(iter(layers.values())).shape
    rows, cols = rasterio.transform.rowcol(transform, coords[:, 0], coords[:, 1])
    rows = np.asarray(rows)
    cols = np.asarray(cols)

    inside = (rows >= 0) & (rows < height) & (cols >= 0) & (cols < width)
    if not inside.all():
        print(f"  ! {int((~inside).sum()):,} of {len(coords):,} points fall "
              "outside the raster")
    rows = np.clip(rows, 0, height - 1)
    cols = np.clip(cols, 0, width - 1)

    names = list(layers)
    matrix = np.column_stack([layers[name][rows, cols] for name in names])
    matrix[~inside] = np.nan

    complete = np.isfinite(matrix).all(axis=1)
    print(f"  {matrix.shape[0]:,} rows x {matrix.shape[1]} features, "
          f"{complete.mean():.1%} complete")
    for i, name in enumerate(names):
        missing = float(np.isnan(matrix[:, i]).mean())
        if missing > 0.01:
            print(f"    {name:18} {missing:5.1%} missing")
    return matrix, names, complete

rng = np.random.default_rng(1)
pick = rng.choice(len(coords), 4000, replace=False)
points = coords[pick]
terrain, terrain_names, complete = sample_layers(layers, points, dem_transform)

  4,000 rows x 10 features, 100.0% complete

Reporting the complete-row fraction matters because most models drop incomplete rows. A feature that is 20% missing silently removes a fifth of the training data, and the removal is not random — it follows wherever that layer has gaps.

Example 3 — distance and density features from vector layers

import numpy as np
import geopandas as gpd
from scipy.spatial import cKDTree


def vector_features(coords, layers, radii=(250, 1000), log_distance=True):
    """Distance to, and count within, each reference layer."""
    features, names = [], []

    for name, gdf in layers.items():
        if gdf.geom_type.isin(["Point", "MultiPoint"]).all():
            targets = np.column_stack([gdf.geometry.x, gdf.geometry.y])
        else:
            points = gdf.geometry.representative_point()
            targets = np.column_stack([points.x, points.y])

        tree = cKDTree(targets)
        distances, _ = tree.query(coords, k=1)
        features.append(np.log1p(distances) if log_distance else distances)
        names.append(f"logdist_{name}" if log_distance else f"dist_{name}")

        for radius in radii:
            counts = tree.query_ball_point(coords, r=radius, return_length=True)
            features.append(np.asarray(counts, dtype=float))
            names.append(f"{name}_within_{radius}m")

    matrix = np.column_stack(features)
    for i, name in enumerate(names):
        column = matrix[:, i]
        print(f"  {name:24} median {np.median(column):9.2f}  "
              f"zero {np.mean(column == 0):5.1%}")
    return matrix, names

vectors, vector_names = vector_features(points, {"river": rivers, "road": roads})

# how far off is the representative-point shortcut for lines?
from shapely import STRtree, points as to_points

exact = {}
for name, lines in (("river", rivers), ("road", roads)):
    tree = STRtree(lines.geometry.values)
    probes = to_points(points)
    nearest = tree.nearest(probes)
    exact[name] = np.array([p.distance(lines.geometry.values[i]) for p, i in zip(probes, nearest)])
    shortcut = np.expm1(vectors[:, vector_names.index(f"logdist_{name}")])
    print(f"  {name}: exact median {np.median(exact[name]):,.0f} m, representative-point median "
          f"{np.median(shortcut):,.0f} m, overstated {np.median(shortcut / np.maximum(exact[name], 1)):.1f}x")

  logdist_river            median      5.44  zero  0.0%
  river_within_250m        median      1.00  zero 45.9%
  river_within_1000m       median     27.00  zero  0.2%
  logdist_road             median      6.56  zero  0.0%
  road_within_250m         median      0.00  zero 80.1%
  road_within_1000m        median      4.00  zero 36.4%
  river: exact median 121 m, representative-point median 229 m, overstated 1.6x
  road: exact median 624 m, representative-point median 705 m, overstated 1.1x

Log-transforming distances is usually right: the difference between 10 m and 100 m from a river matters far more than between 10 km and 10.1 km, and a raw distance column is heavily skewed.

For line and polygon layers, using representative points is an approximation, and for long lines a poor one: a river's representative point is one spot somewhere along it, so the distance to it is the distance to that spot, not to the river. Where distance to a line matters, use the geometry's .distance() with a spatial index, as the check above does — slower, and correct.

Example 4 — do the features travel?

The test that decides whether any of this worked: the same random forest with four feature sets, scored by ordinary random cross-validation and by blocked cross-validation on 2 km squares, where each fold predicts squares it has never seen.

from sklearn.ensemble import RandomForestRegressor
from sklearn.model_selection import GroupKFold, KFold, cross_val_score

target = y[pick]
blocks = (np.floor(points[:, 0] / 2000) * 1000 + np.floor(points[:, 1] / 2000)).astype(int)
sets = {
    "terrain (5 columns)": X[pick],
    "terrain + raw x, y": np.column_stack([X[pick], points]),
    "multi-scale terrain": terrain,
    "multi-scale + exact distances": np.column_stack([terrain, np.log1p(exact["river"]),
                                                      np.log1p(exact["road"])]),
}
model = RandomForestRegressor(n_estimators=200, min_samples_leaf=5, n_jobs=4, random_state=0)
for name, features in sets.items():
    ok = np.isfinite(features).all(axis=1)
    random_r2 = cross_val_score(model, features[ok], target[ok], cv=KFold(5, shuffle=True, random_state=0)).mean()
    block_r2 = cross_val_score(model, features[ok], target[ok], cv=GroupKFold(5), groups=blocks[ok]).mean()
    west_half = points[ok, 0] < np.median(points[:, 0])
    fitted = model.fit(features[ok][west_half], target[ok][west_half])
    east_r2 = fitted.score(features[ok][~west_half], target[ok][~west_half])
    print(f"{name:30} random CV R² {random_r2:.3f}   block CV {block_r2:.3f}   "
          f"trained west, tested east {east_r2:.3f}")

terrain (5 columns)            random CV R² 0.372   block CV 0.241   trained west, tested east 0.009
terrain + raw x, y             random CV R² 0.507   block CV 0.302   trained west, tested east 0.215
multi-scale terrain            random CV R² 0.387   block CV 0.256   trained west, tested east -0.014
multi-scale + exact distances  random CV R² 0.422   block CV 0.290   trained west, tested east -0.001

Read the last column first. Trained on the western half of the area and tested on the eastern, every terrain feature set — five columns, ten multi-scale columns, or those plus exact distances to rivers and roads — scored an R² near zero: whatever relates greenness to terrain in the west does not hold in the east. Raw coordinates scored best there, because greenness has a north–south trend that continues across the split and a forest can carry the edge of the west's pattern a little way east. Within the area, the multi-scale features and distances did help, by 0.05 of R² under random CV and a little under blocked CV.

So the features are worth building, and not a guarantee. The transfer test is the only thing that says whether a model can be used somewhere it was not trained, and on this problem the honest answer is that it cannot.

The neighbourhood scale is the feature

Why the scale of a neighbourhood feature is the feature

Elevation minus a 60 m mean measures local roughness. Elevation minus a 540 m mean measures hillslope position. They are different variables that happen to share a formula.

Computing several and letting the model choose means the importance ranking tells you which scale the process operates at — a genuine finding, and one you cannot get by picking a single window in advance.

The cost is correlated columns. Tree models handle that adequately; linear models do not, and need either regularisation or a single scale.

Why relative features should transfer better — and when they do not

A model fitted on elevations from 50 to 1,000 m has never seen 2,000 m. Applied to a higher range, a tree saturates and a linear model extrapolates.

Topographic position — the difference from the local mean — is centred near zero everywhere. A value of +30 m means the same thing in the Alps and the Pennines, so a model using it can transfer where one using raw elevation cannot. It does only if the relationship itself is the same in both places. In Example 4 it was not: adding multi-scale position and roughness took the east-half R² from 0.009 to −0.014. A relative feature removes one reason a model fails to transfer; it cannot remove a change in what the terrain means for the vegetation.

The same reasoning applies to reflectance relative to a scene median, and to any variable with a regional level plus local variation.

Why nodata handling is not a detail

Two failure modes, both silent.

A nodata fill value sampled as a real value — an elevation of −9999 — becomes the most extreme value in the dataset and dominates every split.

A neighbourhood filter propagates NaN across its whole window, so one bad cell voids a 27 × 27 block. With several scales, the union of those voids can remove a large fraction of the training rows.

Filling before filtering and masking afterwards is the standard compromise, and it should be recorded because it introduces a small bias towards the fill value near gaps.

Why to build features once and reuse them

Feature construction over a large raster is expensive and deterministic. Computing it inside a cross-validation loop repeats the work per fold, and computing it separately for training and prediction risks the two diverging.

Build the feature stack once, as a multi-band raster with the band names recorded. Sample it for training, and read it whole for prediction. Then the training and prediction features are the same by construction.

Web Mercator slope, NaN propagation and log distances

• Terrain derivatives need a projected CRS. Web Mercator understates slope by 1/cos(latitude).
• Convert nodata to NaN before sampling, or the fill value becomes a feature value.
• Neighbourhood filters propagate NaN. Fill, filter, then mask back.
• Log-transform distance features; they are heavily skewed.
• Report the complete-row fraction — models drop incomplete rows non-randomly.
• Two or three neighbourhood scales is enough; more adds correlated columns.
• Representative points approximate distance to polygons. Use .distance() where it matters.
• Build the feature stack once and reuse it for training and prediction.

Internal links

• Spatial features explained: turning geometry into columns — which features and why
• How to sample training points from rasters and polygons — where the rows come from
• Spatial machine learning explained — the wider context
• Feature importance says coordinates are the best predictor — the coordinate trade-off
• How to calculate slope and aspect from a DEM in Python — the terrain derivatives
• How to extract raster values at point locations with rasterio — the sampling mechanics
• How to turn model predictions back into a raster — using the same stack for prediction
• Extrapolation in space explained: the area of applicability — where the features stop being valid

FAQ

What spatial features should I build?

Terrain derivatives, multi-scale neighbourhood summaries, distances to reference features, and shape descriptors. Prefer relative measures over absolute ones.

What neighbourhood sizes should I use?

Two or three spanning an order of magnitude — at 20 m cells, something like 60 m, 180 m and 540 m. The importance ranking then reveals the scale of the process.

Why does my slope look wrong?

Almost certainly a geographic or Web Mercator CRS. Slope is vertical over horizontal distance, and both must be in metres.

How do I handle nodata in neighbourhood features?

Fill before filtering and mask afterwards. Filters propagate NaN across the whole window, so one bad cell voids a large neighbourhood.

Should I log-transform distance features?

Usually. Raw distances are heavily skewed, and the difference between 10 m and 100 m matters far more than between 10 km and 10.1 km.

Why report the complete-row fraction?

Because models drop incomplete rows, and the dropped rows follow wherever a layer has gaps — so the removal is systematic, not random.

Should I compute features inside the cross-validation loop?

No. Build the stack once and sample it. Recomputing per fold repeats deterministic work and risks training and prediction features diverging.