⚙️ How-To Spatial machine learning Spatial statistics Spatial Machine Learning · step 5 of 14 Updated

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, the same random forest with two extra columns:

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

Two columns, a 29% improvement. And the same two columns gained far less under block cross-validation — 0.221 to 0.324 — because a coordinate does not travel.

The job is to build features that carry the same information and travel: distances, neighbourhood summaries, shape descriptors and terrain derivatives.

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
Sample locations feeding terrain derivatives, multi-scale neighbourhood summaries and distance features into one matrix.
Every column is a modelling decision. The scale of each neighbourhood feature is part of it.

Step-by-step solution

1. Sample rasters correctly 

rows, cols = rasterio.transform.rowcol(transform, coords[:, 0], coords[:, 1])
values = 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.

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.

Topographic position index at 60 m, 180 m and 540 m windows, describing local roughness, hillslope position and catchment position.
The same variable at three scales is three predictors. The importance ranking then reveals the scale of the process.

Code examples

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

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

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

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. Where exact distance to a boundary matters, use the geometry's .distance() with a spatial index instead — slower, and correct.

Explanation

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 transfer better

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.

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.

A neighbourhood filter propagating NaN across its whole window against filling before filtering and masking afterwards.
With three scales, the union of voided windows can remove a large share of the training rows.

Edge cases or notes

  • 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.

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.