🧭 Concept Spatial machine learning Spatial Machine Learning · step 3 of 14 Updated

Spatial Features Explained: Turning Geometry into Columns

Problem statement

A machine learning model takes a table of numbers. A geometry is not a number, so every spatial model needs geometry turned into columns — and that translation is where the modelling actually happens.

The choice matters more than the algorithm. Measured on a real problem, the same random forest with different feature sets:

terrain only (slope, elevation, aspect, TPI)   true map R²  0.559
terrain + raw x, y                             true map R²  0.721

A 29% improvement from two columns. And under block cross-validation, which simulates applying the model elsewhere, the same two columns gained far less — 0.221 to 0.324 — because a coordinate means nothing outside the area it was learned in.

Features that encode where a thing is beat features that encode which place it is, wherever the model has to move.

Quick answer

Four families, in increasing order of transferability:

# 1. raw position — powerful here, meaningless elsewhere
X = np.column_stack([x, y])

# 2. geometry properties — travel well
X = np.column_stack([gdf.geometry.area, gdf.geometry.length,
                     4 * np.pi * gdf.geometry.area / gdf.geometry.length ** 2])

# 3. distance to things — travel well and are interpretable
X = np.column_stack([distance_to_coast, distance_to_road, distance_to_water])

# 4. neighbourhood summaries — travel well, and need a scale
X = np.column_stack([mean_elevation_500m, density_1km, majority_class_200m])
Four families of spatial features: raw coordinates, geometry properties, distances to features, and neighbourhood summaries.
Only the first is a lookup table for one place. The other three describe conditions that recur.

Step-by-step solution

1. Geometry properties: what the shape is 

gdf["area"] = gdf.geometry.area                        # projected CRS
gdf["perimeter"] = gdf.geometry.length
gdf["compactness"] = 4 * np.pi * gdf["area"] / gdf["perimeter"] ** 2
gdf["elongation"] = gdf.geometry.minimum_rotated_rectangle.apply(elongation)
gdf["vertices"] = gdf.geometry.apply(lambda g: len(g.exterior.coords))

Compactness — 1 for a circle, lower for anything irregular — distinguishes a warehouse from a terrace row without any coordinate. Vertex count distinguishes hand-digitised outlines from automatically extracted ones, which is sometimes exactly the signal you want and sometimes leakage about the data source.

All of these require a projected CRS. Area in square degrees varies with latitude.

2. Distance to things: what the setting is 

from scipy.spatial import cKDTree

distances, _ = cKDTree(coast_points).query(sample_coords, k=1)

Distance to the coast, to a river, to a road, to the nearest town. These are the most transferable spatial features there are: 500 m from a river means something similar everywhere.

Use cKDTree for point targets and a spatial index with .distance() for lines and polygons. Both need projected coordinates.

3. Neighbourhood summaries: what is around 

mean_elevation = uniform_filter(dem, size=window)
neighbour_count = tree.query_ball_point(coords, r=500, return_length=True)

The scale is the parameter, and it is a modelling decision. Elevation averaged over 100 m and over 5 km are different variables — one is local relief, the other is regional setting.

Compute several scales and let the model choose. Two or three windows spanning an order of magnitude is usually enough.

4. Raw coordinates: last, and knowingly

They work, and they do not travel. Include them only when the model will stay inside the area they were fitted on — and record that decision, because it constrains every future use of the model.

5. Check that the features exist where you will predict

A feature computed from a road layer that covers only the study area is unavailable elsewhere. That is a deployment blocker, not a modelling detail, and it is worth checking before building anything.

Elevation averaged over a 100 metre window against a 5 km window, giving local relief and regional setting respectively.
The window size is the feature. Two scales of the same variable are two different predictors.

Code examples

Example 1 — geometry properties in one pass 

import numpy as np
import geopandas as gpd


def geometry_features(gdf, crs=None):
    """Shape descriptors that do not depend on where the feature is."""
    utm = gdf.to_crs(crs or gdf.estimate_utm_crs())
    geometry = utm.geometry

    out = gpd.GeoDataFrame(index=gdf.index)
    out["area"] = geometry.area
    out["perimeter"] = geometry.length
    out["compactness"] = np.where(
        out["perimeter"] > 0,
        4 * np.pi * out["area"] / out["perimeter"] ** 2, np.nan)

    hull = geometry.convex_hull
    out["convexity"] = np.where(hull.area > 0, out["area"] / hull.area, np.nan)

    boxes = geometry.minimum_rotated_rectangle()
    out["elongation"] = boxes.apply(_elongation)
    out["vertices"] = geometry.apply(
        lambda g: 0 if g is None else len(g.exterior.coords)
        if g.geom_type == "Polygon" else np.nan)

    print(f"  {len(out)} rows, {out.columns.tolist()}")
    print(out.describe().loc[["mean", "50%", "std"]].round(3).to_string())
    return out.drop(columns="geometry", errors="ignore")


def _elongation(box):
    if box is None or box.is_empty:
        return np.nan
    coords = np.array(box.exterior.coords[:-1])
    sides = np.linalg.norm(np.diff(np.vstack([coords, coords[:1]]), axis=0),
                           axis=1)
    long_side, short_side = max(sides[:2]), min(sides[:2])
    return float(long_side / short_side) if short_side else np.nan

Printing the summary statistics is worth the line. A compactness column that is entirely NaN means zero-length perimeters, which means empty or degenerate geometry that will break the model with a confusing error much later.

Example 2 — distances and neighbourhood counts 

import numpy as np
from scipy.spatial import cKDTree


def context_features(sample_coords, layers, radii=(200, 500, 1000)):
    """Distance to, and count within, each of several reference layers.

    layers: {"road": coords_array, "water": coords_array, ...}
    """
    features, names = [], []

    for name, targets in layers.items():
        tree = cKDTree(np.asarray(targets))

        distances, _ = tree.query(sample_coords, k=1)
        features.append(distances)
        names.append(f"dist_{name}")

        for radius in radii:
            counts = tree.query_ball_point(sample_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)
    print(f"  {matrix.shape[1]} features from {len(layers)} layers")
    for i, name in enumerate(names):
        column = matrix[:, i]
        print(f"    {name:22} median {np.median(column):10.1f}  "
              f"zero {np.mean(column == 0):5.1%}")
    return matrix, names

A count feature that is zero for most rows carries little information at that radius. The zero fraction printed here tells you immediately which radii are useful for your data.

Example 3 — raster neighbourhood features at several scales 

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


def raster_features(raster_path, coords, windows=(3, 9, 27), band=1):
    """Sample a raster at several neighbourhood scales."""
    with rasterio.open(raster_path) as src:
        data = src.read(band).astype("float32")
        if src.nodata is not None:
            data[data == src.nodata] = np.nan
        transform = src.transform
        cell = abs(transform.a)

        rows, cols = rasterio.transform.rowcol(transform, coords[:, 0],
                                               coords[:, 1])
        rows = np.clip(np.asarray(rows), 0, src.height - 1)
        cols = np.clip(np.asarray(cols), 0, src.width - 1)

    filled = np.nan_to_num(data, nan=float(np.nanmean(data)))
    features, names = [data[rows, cols]], ["value"]

    for window in windows:
        mean = uniform_filter(filled, size=window)
        squared = uniform_filter(filled ** 2, size=window)
        std = np.sqrt(np.maximum(squared - mean ** 2, 0))

        features.append(mean[rows, cols])
        names.append(f"mean_{int(window * cell)}m")
        features.append(std[rows, cols])
        names.append(f"std_{int(window * cell)}m")
        features.append(data[rows, cols] - mean[rows, cols])
        names.append(f"relative_{int(window * cell)}m")

    matrix = np.column_stack(features)
    print(f"  {matrix.shape[1]} features at {[int(w * cell) for w in windows]} m")
    return matrix, names

The relative_ features — the value minus its neighbourhood mean — are often the most useful of the set. They express "higher or lower than around here", which is scale-explicit and transfers well, whereas the raw value carries the regional level.

Explanation

Why the feature set matters more than the algorithm

Swapping a random forest for gradient boosting typically moves a score by a few percent. Adding a well-chosen covariate can move it far more — the measured coordinate example was 29%.

That is because the algorithm can only combine what it is given. A model with no information about proximity to water cannot learn a relationship with proximity to water, however sophisticated it is.

Time spent on features is therefore usually better spent than time on hyperparameters, and it is also where domain knowledge enters a model at all.

Why distance features transfer and coordinates do not

Both encode position. The difference is whether the encoding means the same thing elsewhere.

"500 m from a river" implies similar drainage, soil moisture and flood risk anywhere. "Easting 422,000" implies nothing outside its own UTM zone, and something completely different in another.

The measured consequence: coordinates gained 0.162 R² in-area and 0.103 under block CV. Distance features gain similarly in both, because the relationship they encode is real rather than positional.

Why the neighbourhood scale is a feature, not a parameter

Elevation averaged over 100 m and over 5 km measure different things: local terrain and regional setting. A model given both can learn that vegetation responds to one and not the other.

Choosing a single "best" window discards that. Computing several and letting the model select is both easier and more informative — and the importance ranking then tells you which scale the process operates at, which is a genuine finding.

Two or three windows spanning an order of magnitude is enough. More adds correlated columns without adding information.

Why to check feature availability before modelling

A model is only deployable where its features can be computed. A feature derived from a detailed road network available for one country cannot be computed for another.

This is worth checking first, because it is a hard constraint rather than an accuracy trade-off. A model that scores 0.9 and cannot be run where it is needed is worth less than one that scores 0.7 and can.

Absolute elevation failing to transfer to a higher range against topographic position, which is centred near zero everywhere.
Anything with a regional level plus local variation should usually be entered as the local part.

Edge cases or notes

  • Project before computing area, length or distance.
  • Compactness needs a non-zero perimeter. Guard against degenerate geometry.
  • Compute several neighbourhood scales and let the model choose.
  • Relative features — value minus neighbourhood mean — often transfer best.
  • Vertex count can be leakage about the digitising source rather than the object.
  • Check feature availability at deployment, not just at training.
  • Distance features are unbounded; consider a log transform for skewed ones.
  • Raw coordinates last, and knowingly.

FAQ

What are spatial features?

Columns derived from geometry: shape properties, distances to reference features, neighbourhood summaries, and raw coordinates. They are how geometry enters a model that takes a table.

Should I include x and y as features?

Only if the model will stay inside the area it was fitted on. Measured, they raised the in-area R² from 0.559 to 0.721 and gained far less under block CV.

Which spatial features transfer best?

Distances to reference features and relative neighbourhood measures. "500 m from a river" means the same thing anywhere; an easting does not.

What neighbourhood size should I use?

Several. Compute two or three windows spanning an order of magnitude and let the model choose; the importance ranking then tells you the scale of the process.

Do I need a projected CRS?

Yes, for anything measuring area, length or distance. Square degrees vary with latitude.

Is vertex count a useful feature?

Sometimes, and it can be leakage — it often reflects how the data was digitised rather than what the object is.

How many features should I build?

Enough to express the domain knowledge you have. Feature choice usually moves a score more than algorithm choice does.