Spatial Features Explained: Turning Geometry into Columns
On this page
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 — Sentinel-2 greenness predicted from terrain over part of Snowdonia — the same random forest with different feature sets:
terrain only (slope, elevation, aspect, TPI) in-area R² 0.348
terrain + raw x, y in-area R² 0.509
A large gain from two columns. Under block cross-validation, which withholds whole areas, the same two columns gained far less — 0.258 to 0.283 — because a coordinate describes one place, not a condition.
Features that describe the conditions at a place should travel where coordinates cannot. Whether they do is a measurement: on the same problem, trained on one half of the area and tested on the other, no terrain feature set transferred at all.
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])Shape, distance, neighbourhood, then coordinates
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.
Shape columns, distances and multi-scale neighbourhoods
The data used on this page
The problem from spatial machine learning explained — Sentinel-2 greenness and the Copernicus DEM on a 20 m grid over part of Snowdonia — with OpenStreetMap's rivers, roads and buildings for the same area.
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 as a file, and OpenStreetMap rivers, roads and buildings 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, Polygon
QUERIES = {
"snowdonia_osm.json": """(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););""",
"snowdonia_buildings.json": """way["building"](53.03,-4.16,53.12,-4.00);""",
}
for path, body in QUERIES.items():
if not os.path.exists(path):
reply = requests.post("https://overpass-api.de/api/interpreter", timeout=180,
data={"data": f"[out:json][timeout:120];{body}out geom tags;"},
headers={"User-Agent": "spatialworkflow.io guide"})
reply.raise_for_status()
json.dump(reply.json(), open(path, "w"))
lines = 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)
buildings = gpd.GeoDataFrame(
[{"building": e["tags"].get("building"), "geometry": Polygon([(p["lon"], p["lat"]) for p in e["geometry"]])}
for e in json.load(open("snowdonia_buildings.json"))["elements"] if len(e["geometry"]) >= 4], crs=4326)
print(f"{len(lines):,} river and road ways, {len(buildings):,} buildings "
f"({(buildings['building'] == 'ruins').sum()} of them ruins)")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,518 river and road ways, 1,026 buildings (151 of them ruins)
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.nanshapes = geometry_features(buildings)
shapes["ruin"] = buildings["building"].eq("ruins")
print(shapes.groupby("ruin")[["area", "compactness", "elongation"]].median().round(2).to_string()) 1026 rows, ['area', 'perimeter', 'compactness', 'convexity', 'elongation', 'vertices']
area perimeter compactness convexity elongation vertices
mean 110.663 43.783 0.641 0.933 1.891 7.418
50% 70.443 37.333 0.711 1.000 1.508 5.000
std 157.736 31.563 0.179 0.143 1.594 5.155
area compactness elongation
ruin
False 77.74 0.71 1.52
True 12.82 0.63 1.38
On real buildings the shape columns already carry a signal: the 151 mapped ruins here have a median footprint of 12.8 m² against 77.7 m² for the rest, with no coordinate involved.
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, namesThe reference layers here are lines, so their vertices stand in for them — much closer than one representative point per line, and still an approximation where vertices are far apart:
rng = np.random.default_rng(1)
sample_coords = coords[rng.choice(len(coords), 4000, replace=False)]
vertices = {kind: np.vstack([np.asarray(g.coords) for g in group.geometry])
for kind, group in lines.groupby("kind")}
context, context_names = context_features(sample_coords, vertices, radii=(200, 1000)) 6 features from 2 layers
dist_river median 122.7 zero 0.0%
river_within_200m median 11.0 zero 32.5%
river_within_1000m median 553.5 zero 0.1%
dist_road median 624.0 zero 0.0%
road_within_200m median 0.0 zero 77.7%
road_within_1000m median 57.0 zero 33.0%
The vertex shortcut held up here: the median distance to a river came out at 122.7 m, against 121 m computed exactly to the lines on the feature-building page. A count feature that is zero for most rows carries little information at that radius — road_within_200m is zero for 77.7% of cells. 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, namesrelief, relief_names = raster_features("elev.tif", sample_coords)
print(np.round(np.nanpercentile(relief, [5, 50, 95], axis=0), 1)[:, :4], relief_names[:4]) 10 features at [60, 180, 540] m
[[126. 126.6 1.1 -1.4]
[383.9 383.6 5.1 -0. ]
[791.5 791.3 13.2 1.5]] ['value', 'mean_60m', 'std_60m', 'relative_60m']
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.
Why the feature set beats the algorithm
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 — in the measured example two coordinate columns took the in-area R² from 0.348 to 0.509.
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, on the same problem: coordinates gained 0.135 of R² under random cross-validation and 0.061 under 2 km blocks, while multi-scale terrain and distances to rivers and roads gained about 0.05 under both. That is the pattern the argument predicts. But trained on the western half and tested on the eastern, every terrain feature set fell to an R² near zero and the coordinates reached 0.215, because greenness there has a north–south trend that crosses the split. Transferable in principle is not transferable in fact; only the held-out test says which.
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.
Perimeters, relative features and vertex-count leakage
- 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.
Internal links
- How to build spatial features for a machine learning model — the implementation
- Spatial machine learning explained — why location is special
- Feature importance says coordinates are the best predictor — the coordinate trade-off measured
- How to do spatial cross-validation in Python — measuring transferability
- How to sample training points from rasters and polygons — where the rows come from
- How to extract raster values at point locations with rasterio — the sampling mechanics
- How to find the nearest point in GeoPandas — distance features
- Extrapolation in space explained: the area of applicability — where the features stop being valid
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.348 to 0.509 and gained far less under block CV — and whether anything transfers has to be tested, not assumed.
Which spatial features transfer best?
In principle, distances to reference features and relative neighbourhood measures: "500 m from a river" means the same thing anywhere; an easting does not. In practice, test it — on a measured Snowdonia example none of the terrain features transferred from one half of the area to the other.
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.
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.