Feature Engineering for Pixel-Based Models

A pixel-based model sees a vector of numbers and nothing else. Every bit of geography that matters — that this pixel sits on a north-facing slope, that its neighbours are all bright, that it was green in June and bare in October — has to be encoded into that vector or it is invisible to the model. Feature engineering for raster data is the work of doing that encoding deliberately, in an order that can be reproduced at inference time. This topic is part of Raster Machine Learning & Model Inference.

The organising principle is simple: a feature stack is a raster, so all the usual raster rules apply. Every feature band must be on the same grid, with the same nodata convention, in a fixed and documented order. The moment a feature is computed on a different grid — a DEM at 30 m joined to Sentinel-2 at 10 m without resampling, say — the whole matrix is quietly wrong.

---

Prerequisites

pip install "rasterio>=1.3.0" "xarray>=2023.12" "rioxarray>=0.15" "numpy>=1.23" "scipy>=1.11"

You will want the band-math grounding from band math operations with xarray, the alignment tooling from aligning two rasters with reproject_match, and a cloud mask from masking clouds with the Sentinel-2 SCL band.

---

Step-by-Step Workflow

Step 1 — Fix the base bands and the scaling

Everything downstream is derived from a small set of calibrated reflectance bands. Decide which ones, decide the scale factor, and write both into the manifest before computing anything.

import numpy as np
import rioxarray
import xarray as xr

BASE_BANDS = ["B02", "B03", "B04", "B08", "B11", "B12"]
REFLECTANCE_SCALE = 10_000.0


def load_base(paths: dict[str, str]) -> xr.Dataset:
    """Load the base reflectance bands onto one grid as float32."""
    ref = rioxarray.open_rasterio(paths[BASE_BANDS[0]], masked=True).squeeze("band")
    out = {}
    for name in BASE_BANDS:
        da = rioxarray.open_rasterio(paths[name], masked=True).squeeze("band")
        if da.shape != ref.shape:
            da = da.rio.reproject_match(ref)     # 20 m bands onto the 10 m grid
        out[name] = (da / REFLECTANCE_SCALE).astype("float32")
    return xr.Dataset(out)

The reproject_match call is doing quiet but critical work: Sentinel-2’s shortwave bands ship at 20 m and must land on the 10 m grid before they can share a feature vector with the visible bands. The resampling choice there is discussed in resampling Sentinel-2 20 m bands to 10 m.

Step 2 — Add spectral indices that separate your classes

Indices are ratios that cancel out illumination and emphasise a physical contrast. Adding every index in the literature is a waste; adding three that map onto the classes you care about is most of the signal.

import xarray as xr


def add_indices(ds: xr.Dataset) -> xr.Dataset:
    """Append the normalised-difference features a land cover model needs."""
    def nd(a: str, b: str) -> xr.DataArray:
        return ((ds[a] - ds[b]) / (ds[a] + ds[b])).astype("float32")

    ds["NDVI"] = nd("B08", "B04")     # vegetation vigour
    ds["NDWI"] = nd("B03", "B08")     # open water
    ds["NDMI"] = nd("B08", "B11")     # canopy moisture
    ds["NDBI"] = nd("B11", "B08")     # built-up surfaces
    ds["BSI"] = (((ds.B11 + ds.B04) - (ds.B08 + ds.B02))
                 / ((ds.B11 + ds.B04) + (ds.B08 + ds.B02))).astype("float32")
    return ds

Note that NDBI is exactly -NDMI. Including both is harmless for a tree model and actively misleading for a linear one, and it is the kind of redundancy that creeps into a feature list over months. The index definitions themselves are covered in depth in spectral index calculation pipelines.

Step 3 — Encode local context with focal statistics

This is the step that closes most of the gap between a pixel model and a small convolutional one. A focal mean tells the model what the neighbourhood looks like; a focal standard deviation tells it how textured the neighbourhood is. Together they distinguish a smooth field from a broken-up settlement at the same average brightness.

import numpy as np
from scipy import ndimage


def focal_stats(arr: np.ndarray, size: int = 5) -> tuple[np.ndarray, np.ndarray]:
    """NaN-aware focal mean and standard deviation over a square window."""
    valid = np.isfinite(arr).astype("float32")
    filled = np.where(np.isfinite(arr), arr, 0.0).astype("float32")

    count = ndimage.uniform_filter(valid, size=size) * size * size
    total = ndimage.uniform_filter(filled, size=size) * size * size
    mean = np.divide(total, count, out=np.full_like(total, np.nan), where=count > 0)

    sq = ndimage.uniform_filter(filled ** 2, size=size) * size * size
    var = np.divide(sq, count, out=np.full_like(sq, np.nan), where=count > 0) - mean ** 2
    return mean, np.sqrt(np.clip(var, 0, None))

The NaN handling is not optional decoration. A plain uniform_filter over an array containing NaN propagates NaN across the whole window, so a single missing pixel erases a 5×5 block of features. The count-and-divide formulation above degrades gracefully instead, returning the mean of whatever valid pixels the window contains.

Two focal windows — a tight 3×3 and a broad 11×11 — usually beat one, because they encode texture at two scales. Three is rarely worth the read cost.

Step 4 — Join terrain and temporal features

Terrain features come from a DEM that almost never matches the imagery grid, so the join is a resampling operation. Temporal features come from collapsing a time series into statistics that survive missing observations.

import rioxarray
import xarray as xr


def add_terrain(ds: xr.Dataset, dem_path: str) -> xr.Dataset:
    """Resample a DEM onto the reference grid and derive slope and aspect."""
    ref = ds[BASE_BANDS[0]]
    dem = rioxarray.open_rasterio(dem_path, masked=True).squeeze("band")
    dem = dem.rio.reproject_match(ref).astype("float32")

    gy, gx = np.gradient(dem.values, abs(ref.rio.resolution()[0]))
    ds["elevation"] = dem
    ds["slope"] = (("y", "x"), np.degrees(np.arctan(np.hypot(gx, gy))).astype("float32"))
    ds["aspect_sin"] = (("y", "x"), np.sin(np.arctan2(-gx, gy)).astype("float32"))
    ds["aspect_cos"] = (("y", "x"), np.cos(np.arctan2(-gx, gy)).astype("float32"))
    return ds

Aspect is encoded as a sine and cosine pair rather than as degrees, because raw aspect is circular: 359° and 1° are neighbours on the ground and maximally distant as numbers. Every tree that splits on raw aspect wastes a split recovering that. The same transformation applies to day-of-year in a temporal stack.

For the temporal half, collapse the series into percentiles rather than means. The median and the 10th/90th percentiles of NDVI across a season describe the phenology far better than an average, and they are robust to the handful of cloudy observations that always survive masking:

import xarray as xr


def add_temporal(ds: xr.Dataset, ndvi_series: xr.DataArray) -> xr.Dataset:
    """Summarise a masked NDVI time series into robust seasonal features."""
    ds["ndvi_p10"] = ndvi_series.quantile(0.10, dim="time").astype("float32")
    ds["ndvi_p50"] = ndvi_series.quantile(0.50, dim="time").astype("float32")
    ds["ndvi_p90"] = ndvi_series.quantile(0.90, dim="time").astype("float32")
    ds["ndvi_range"] = (ds.ndvi_p90 - ds.ndvi_p10).astype("float32")
    ds["ndvi_count"] = ndvi_series.notnull().sum(dim="time").astype("float32")
    return ds

ndvi_count is a feature in its own right, and an unusual one: it tells the model how much evidence each pixel’s summaries rest on. Pixels with three valid observations are systematically noisier than pixels with twenty, and a tree model will happily learn to discount them.

Step 5 — Write a named stack

The stack leaves this stage as a raster with band descriptions, because a feature order that lives only in a Python list is an outage waiting to happen.

import numpy as np
import rasterio
import xarray as xr


def write_feature_stack(ds: xr.Dataset, path: str, names: list[str]) -> None:
    ref = ds[names[0]]
    profile = {
        "driver": "GTiff", "height": ref.sizes["y"], "width": ref.sizes["x"],
        "count": len(names), "dtype": "float32", "crs": ds.rio.crs,
        "transform": ds.rio.transform(), "nodata": np.nan,
        "compress": "zstd", "zstd_level": 3, "tiled": True,
        "blockxsize": 512, "blockysize": 512,
    }
    with rasterio.open(path, "w", **profile) as dst:
        for i, name in enumerate(names, start=1):
            dst.write(ds[name].values.astype("float32"), i)
            dst.set_band_description(i, name)

---

Which Features Actually Earn Their Place

Feature lists grow monotonically unless someone prunes them, and every addition is paid for on every pixel of every scene forever. The useful discipline is to rank features by permutation importance on a held-out spatial fold — not by the model’s built-in impurity importance, which flatters high-cardinality continuous features — and then to delete from the bottom until accuracy moves.

Three patterns show up again and again in these rankings. Temporal summaries beat single-date reflectance for anything that grows, because phenology is the strongest signal in the data and a single date throws it away. Shortwave infrared punches above its weight for soil, moisture and burn discrimination, which is why a six-band base that omits B11 and B12 handicaps a model before it starts. And raw visible bands are usually the first things to go once indices are present, because the index has already extracted the contrast that mattered and discarded the illumination that did not.

The counter-intuitive entry is ndvi_count. It ranks low here, which is the healthy outcome; when it ranks high it is a warning, as the troubleshooting section below explains.

Pruning has a second benefit that only shows up at scale. A twelve-feature stack read from object storage is twelve byte-range requests per window; a twenty-eight-feature stack is twenty-eight. On a continental run that difference is hours, and it lands squarely on the egress bill described in reducing S3 egress costs in raster pipelines.

---

Parameter Reference

---

Verification & Testing

Check three properties before a feature stack is used for anything.

import numpy as np
import rasterio

with rasterio.open("features.tif") as src:
    names = list(src.descriptions)
    assert all(names), "every band must carry a description"
    assert len(set(names)) == len(names), "duplicate feature names"

    sample = src.read(window=((0, 512), (0, 512)))
    finite = np.isfinite(sample).reshape(len(names), -1).mean(axis=1)
    for name, frac in zip(names, finite):
        print(f"{name:<14} {frac:6.1%} finite")

A feature that is 5% finite is a bug, not a feature: it usually means a source raster covered only part of the scene. A feature that is 100% finite where every other band has gaps is equally suspicious — something filled the nodata silently.

The second test is a correlation sweep. Compute pairwise correlations on a sample of pixels and look for pairs above 0.98; those are duplicated information, like the NDBI/NDMI pair above, and dropping one costs nothing.

The third test is the one that catches drift: recompute a single feature for one window using the independent, naive implementation and compare. If the fast focal filter and a plain loop over a 5×5 neighbourhood disagree by more than floating-point noise, the fast path has a boundary bug.

---

Troubleshooting

Every focal feature is NaN near the scene edge

ndimage pads with zeros by default, and the count-based correction then divides by a partial count — correct, but only if mode="constant" with cval=0 is used consistently for both numerator and denominator. Mixing mode="reflect" for one and the default for the other produces exactly this artefact.

Feature values differ between training and inference

The band order changed. This is why descriptions are written into the file and asserted on load. If the assertion is missing, add it before debugging anything else.

MemoryError building the stack for a full tile

A 10,980 × 10,980 scene with 28 float32 features is 13 GB. Build the stack per window rather than per scene, or write bands incrementally as the code above does and never hold more than one feature at a time.

Terrain features look blocky compared to the imagery

A 30 m DEM resampled to 10 m with nearest-neighbour produces 3×3 plateaus. Use bilinear for continuous terrain, and consider computing slope at the DEM’s native resolution before resampling, which is more faithful than computing it from the upsampled surface.

The model ranks ndvi_count as the most important feature

That is usually real and usually a problem: the model has learned which regions are cloudy rather than what is on the ground. Check whether cloud frequency correlates with your classes geographically, and if it does, drop the count feature and rebalance the training sample.

---

Frequently Asked Questions

Q: Do tree-based models need features to be scaled? No. Decision trees split on thresholds, so any monotonic rescaling leaves the model unchanged. Scaling matters for linear models, neural networks and anything distance-based, and it matters for storage regardless, because a float32 stack is four times the size of a scaled int16 one.

Q: How many features is too many for a raster model? The practical ceiling is set by inference cost, not by the model. Every feature is a band that must be read, computed and held in memory for every pixel of every scene. Thirty well-chosen features usually beat a hundred and twenty, and they run four times faster.

Q: Should spectral indices be computed before or after cloud masking? After. A cloudy pixel has a physically meaningless index value, and if it survives into a temporal summary it will drag the median or the percentile with it. Mask first, then compute, then aggregate.

Q: Is it worth adding features from a second sensor? Often, but only after harmonisation. Landsat and Sentinel-2 reflectance differ enough that an unharmonised join teaches the model to detect which sensor a pixel came from. The correction is covered in harmonising Landsat and Sentinel-2 reflectance.

---

Related

• Stacking Spectral Indices as Model Features — which indices earn their place and how to name them.
• Adding Terrain Derivatives as Predictors — slope, aspect and curvature joined to an imagery grid.
• Building Temporal Feature Stacks from Image Time Series — percentiles, phenology metrics and gap handling.
• Harmonising Features across Sensors for Transfer — making a model trained on one sensor work on another.
• Band Math Operations with xarray — the array mechanics every feature here is built on.