Building a scikit-learn Pipeline for Raster Features

Package raster feature engineering into a scikit-learn Pipeline: custom transformers, ColumnTransformer, spatial CV in GridSearchCV, and one artifact that serves exactly what it trained.

Put every fitted step — imputation, scaling, encoding, reduction — inside a Pipeline, tune it with a group-aware splitter whose groups are spatial blocks, and serialise the whole object. The result is a single artifact that applies at inference exactly the transformations it learned during training, which removes an entire family of training-serving bugs by construction.

This page is the packaging pattern. For the estimator choices and evaluation around it, see Training with scikit-learn for Geospatial Data.

Anything fitted outside the pipeline sees the validation fold Two arrangements are compared. In the first, scaling and imputation are applied to the whole dataset before splitting, so the validation fold contributes to the statistics and the reported score is optimistic. In the second, the same steps live inside the pipeline, so cross-validation refits them on each training fold and the validation fold is never seen. fitted outside fitted inside scaler.fit(ALL rows) train fold validation fold the validation fold contributed to the mean and variance reported score is optimistic Pipeline([scaler, model]) train fold validation fold the transformer is refitted per fold, on training rows only score reflects a deployable process The same argument applies to imputers, encoders and PCA — anything with a .fit().

Why This Fails in Geospatial ML Pipelines

Preprocessing outside the pipeline is the classic leak, and geospatial feature sets make it worse than usual because so many steps have fitted state: a scaler, an imputer whose medians come from the data, a target encoder for administrative regions, a PCA rotation. Fit any of them on the full table and every cross-validation score is optimistic — not by a rounding error, but by however much the held-out fold contributed to the statistics.

The second failure is a split that ignores geography. GridSearchCV defaults to KFold, which shuffles rows. With spatially autocorrelated data that puts a pixel’s neighbours in the training fold while the pixel itself is validated, so hyperparameter search selects whatever memorises best — usually maximum depth and minimum regularisation. Substituting GroupKFold with spatial-block groups changes both the chosen hyperparameters and the reported score, in the direction of reality. The reasoning is set out in spatial cross-validation strategies.

Third, the serving mismatch. When preprocessing lives in a notebook cell and the model is pickled alone, inference has to re-implement the transformations, and the re-implementation drifts. A single serialised Pipeline makes that impossible: the artifact cannot be loaded without its transformers, and it applies them in the order it learned them.

Core Principles

  • Everything with a .fit() goes inside the pipeline.
  • Split by space, not at random. GroupKFold with spatial block groups.
  • Use ColumnTransformer so continuous, categorical and cyclic features get their own treatment.
  • Write custom transformers as BaseEstimator subclasses so they clone correctly during search.
  • Serialise the whole pipeline, and record the library versions with it.
  • Assert the feature contract at load. Column names and order are part of the artifact.

Production-Ready Code

from __future__ import annotations

import logging
import numpy as np
import pandas as pd
from sklearn.base import BaseEstimator, TransformerMixin
from sklearn.compose import ColumnTransformer
from sklearn.impute import SimpleImputer
from sklearn.pipeline import Pipeline
from sklearn.preprocessing import OneHotEncoder, StandardScaler

logger = logging.getLogger(__name__)


class CyclicAspect(BaseEstimator, TransformerMixin):
    """Turn a compass bearing column into northness and eastness.

    Written as a proper estimator so GridSearchCV can clone it per fold.
    """

    def fit(self, X, y=None):
        self.n_features_in_ = X.shape[1]
        return self

    def transform(self, X):
        rad = np.radians(np.asarray(X, dtype="float64"))
        out = np.hstack([np.cos(rad), np.sin(rad)])
        return np.nan_to_num(out, nan=0.0)

    def get_feature_names_out(self, input_features=None):
        base = list(input_features or [f"x{i}" for i in range(self.n_features_in_)])
        return np.array([f"{b}_north" for b in base] + [f"{b}_east" for b in base])


class RatioFeatures(BaseEstimator, TransformerMixin):
    """Normalised differences between band pairs, computed inside the pipeline."""

    def __init__(self, pairs: list[tuple[str, str]] | None = None):
        self.pairs = pairs or [("nir", "red"), ("green", "nir")]

    def fit(self, X: pd.DataFrame, y=None):
        missing = [c for p in self.pairs for c in p if c not in X.columns]
        if missing:
            raise ValueError(f"RatioFeatures needs columns {sorted(set(missing))}")
        return self

    def transform(self, X: pd.DataFrame) -> np.ndarray:
        cols = []
        for a, b in self.pairs:
            denom = X[a].to_numpy("float64") + X[b].to_numpy("float64")
            num = X[a].to_numpy("float64") - X[b].to_numpy("float64")
            with np.errstate(invalid="ignore", divide="ignore"):
                cols.append(np.where(np.abs(denom) > 1e-6, num / denom, np.nan))
        return np.column_stack(cols)

    def get_feature_names_out(self, input_features=None):
        return np.array([f"nd_{a}_{b}" for a, b in self.pairs])


CONTINUOUS = ["blue", "green", "red", "nir", "swir1", "elevation", "slope_deg", "tpi_11"]
CATEGORICAL = ["soil_class"]
CYCLIC = ["aspect_deg"]
RATIO_SOURCE = ["nir", "red", "green"]


def build_pipeline(estimator) -> Pipeline:
    """One artifact: preprocessing plus estimator, in the order they are learned."""
    pre = ColumnTransformer(
        transformers=[
            ("cont", Pipeline([("impute", SimpleImputer(strategy="median")),
                               ("scale", StandardScaler())]), CONTINUOUS),
            ("cat", Pipeline([("impute", SimpleImputer(strategy="most_frequent")),
                              ("onehot", OneHotEncoder(handle_unknown="ignore",
                                                       min_frequency=20))]), CATEGORICAL),
            ("cyc", CyclicAspect(), CYCLIC),
            ("ratio", RatioFeatures([("nir", "red"), ("green", "nir")]), RATIO_SOURCE),
        ],
        remainder="drop",
        verbose_feature_names_out=False,
    )
    return Pipeline([("pre", pre), ("model", estimator)])

Step-by-Step Walkthrough

Step 1 — Build spatial groups alongside the frame

import geopandas as gpd
import numpy as np

samples = gpd.read_file("samples.gpkg").to_crs(32633)
coords = np.c_[samples.geometry.x, samples.geometry.y]

BLOCK_M = 5_000
groups = (np.floor(coords[:, 0] / BLOCK_M).astype("int64") * 100_000
          + np.floor(coords[:, 1] / BLOCK_M).astype("int64"))
print(f"{len(np.unique(groups))} spatial blocks over {len(samples)} samples")

Step 2 — Assemble and smoke-test the pipeline

from sklearn.ensemble import HistGradientBoostingClassifier

X = samples[CONTINUOUS + CATEGORICAL + CYCLIC]
y = samples["land_cover"].to_numpy()

pipe = build_pipeline(HistGradientBoostingClassifier(random_state=0))
pipe.fit(X.head(500), y[:500])
names = pipe.named_steps["pre"].get_feature_names_out()
print(f"{len(names)} engineered features: {list(names)[:8]} ...")

Fitting on 500 rows first catches column-name mistakes in seconds rather than after a full search.

Step 3 — Tune with a group-aware splitter

from sklearn.model_selection import GridSearchCV, GroupKFold

grid = {
    "model__max_depth": [4, 6, None],
    "model__learning_rate": [0.05, 0.1],
    "model__min_samples_leaf": [20, 50],
}
search = GridSearchCV(
    pipe, grid, cv=GroupKFold(n_splits=5), scoring="f1_macro",
    n_jobs=-1, refit=True, verbose=1,
)
search.fit(X, y, groups=groups)

print("best spatial CV f1_macro:", round(search.best_score_, 4))
print("best params:", search.best_params_)

Run the same search once with plain KFold for comparison. The gap between the two best scores is the size of the leak you would have shipped.

Step 4 — Serialise the whole artifact with its versions

import json
import sklearn
import joblib

joblib.dump(search.best_estimator_, "artifacts/landcover_pipeline.joblib")
with open("artifacts/landcover_pipeline.json", "w") as fh:
    json.dump({
        "sklearn": sklearn.__version__,
        "input_columns": list(X.columns),
        "cv": "GroupKFold(5) on 5 km spatial blocks",
        "best_score_f1_macro": float(search.best_score_),
        "best_params": {k: str(v) for k, v in search.best_params_.items()},
    }, fh, indent=2)

Step 5 — Load and enforce the contract at inference

def load_pipeline(model_path: str, meta_path: str):
    meta = json.load(open(meta_path))
    if sklearn.__version__ != meta["sklearn"]:
        logger.warning("scikit-learn %s at serving vs %s at training",
                       sklearn.__version__, meta["sklearn"])
    pipe = joblib.load(model_path)

    def predict(frame: pd.DataFrame):
        missing = [c for c in meta["input_columns"] if c not in frame.columns]
        if missing:
            raise ValueError(f"inference frame is missing columns: {missing}")
        return pipe.predict(frame[meta["input_columns"]])

    return predict
ColumnTransformer routes each column type to the right treatment An input frame is split into four branches. Continuous columns are imputed with the median and standardised. Categorical columns are imputed with the mode and one-hot encoded. The aspect column is converted to northness and eastness. Selected bands feed a ratio transformer. The four outputs are concatenated into a single design matrix that reaches the estimator. One frame in, four treatments, one design matrix out GeoDataFrame 12 columns continuous median impute → standardise categorical mode impute → one-hot cyclic aspect → northness, eastness band pairs → normalised differences hstack 27 features estimator

Verification

import numpy as np
import pandas as pd
import pytest
from sklearn.base import clone
from sklearn.model_selection import GroupKFold, cross_val_score


def test_transformers_are_clonable():
    """A transformer that does not clone silently shares state across folds."""
    t = CyclicAspect()
    assert clone(t).get_params() == t.get_params()


def test_cyclic_transform_is_continuous_across_north():
    t = CyclicAspect().fit(np.array([[0.0]]))
    a = t.transform(np.array([[359.0]]))
    b = t.transform(np.array([[1.0]]))
    assert np.linalg.norm(a - b) < 0.05


def test_spatial_cv_scores_below_random_cv():
    """If they match, the groups are not actually separating anything."""
    pipe = build_pipeline(HistGradientBoostingClassifier(random_state=0))
    spatial = cross_val_score(pipe, X, y, cv=GroupKFold(5), groups=groups,
                              scoring="f1_macro").mean()
    random = cross_val_score(pipe, X, y, cv=5, scoring="f1_macro").mean()
    assert spatial <= random + 1e-6
    print(f"leak size: {random - spatial:.3f} f1_macro")


def test_unknown_category_does_not_crash_inference():
    pipe = build_pipeline(HistGradientBoostingClassifier()).fit(X, y)
    new = X.head(1).copy()
    new["soil_class"] = "a_class_never_seen"
    assert pipe.predict(new).shape == (1,)


def test_missing_column_is_reported_clearly():
    predict = load_pipeline("artifacts/landcover_pipeline.joblib",
                            "artifacts/landcover_pipeline.json")
    with pytest.raises(ValueError, match="missing columns"):
        predict(X.drop(columns=["elevation"]).head(5))

The third test is the one to keep visible: it prints the size of the leak, and that number is the single most useful figure to put in front of anyone asking why the offline score dropped after the pipeline was fixed.

The spatial score is the one that matches production Three bars compare a random k-fold cross-validation score, a spatially grouped cross-validation score, and the observed production performance. The random score is much higher than both, while the spatial score sits very close to the production value. Same pipeline, three numbers 0.89 random KFold 0.71 GroupKFold (blocks) 0.69 observed in production Fixing the split does not make the model worse — it stops the score from lying.

FAQ

Why put feature engineering inside the Pipeline instead of before it?

Because anything fitted outside sees the validation fold, so the reported score describes a process you cannot deploy. Inside, cross-validation refits every transformer on each training fold alone.

Can a Pipeline handle geometry columns?

Not directly — estimators need numeric arrays. Either convert geometry to numeric columns in a transformer at the head of the pipeline, or compute spatial features beforehand and pass the block labels separately as groups, which is what most production pipelines do.

How do I use spatial cross-validation with GridSearchCV?

Pass cv=GroupKFold(n_splits=...) and supply groups of spatial block ids to fit. Whole blocks then stay together in every fold, so hyperparameters are tuned against geographically held-out data. See implementing SpatialKFold in Python for a custom splitter.

Part of: Training with scikit-learn for Geospatial Data Part of: Training Geospatial Predictive Models in Python