Removing Spatial Trend with Geographic Detrending

Remove large-scale spatial trend with trend-surface analysis: fit a polynomial on projected X/Y coordinates, model detrended residuals, then add the surface back.

TL;DR Answer

Fit a low-order polynomial regression on the projected X/Y coordinates of your samples, subtract that fitted surface from the target to get detrended residuals, train your main model on the residuals, then add the trend surface back at prediction time. Removing the broad regional gradient before modelling frees the estimator to learn local relationships instead of memorising geography, and it is one of the core techniques for handling spatial autocorrelation in a training pipeline.

import numpy as np
from sklearn.linear_model import LinearRegression
from sklearn.preprocessing import PolynomialFeatures

# coords: (n, 2) projected easting/northing in metres; z: target values
poly = PolynomialFeatures(degree=1, include_bias=True)
C = poly.fit_transform(coords)
trend = LinearRegression().fit(C, z)

residual = z - trend.predict(C)   # detrended target — model this

Part of: Handling Spatial Autocorrelation

Why a Broad Spatial Trend Hurts Your Model

Most environmental targets carry a large-scale gradient that has little to do with the local predictors you actually care about. Soil organic carbon rises toward wetter northern latitudes; house prices climb toward a city centre; air temperature falls with a regional elevation ramp. That gradient is real, but it is also the easiest signal in the dataset, and a flexible learner will spend its capacity encoding position rather than the causal features you engineered.

The symptom shows up as strong spatial autocorrelation in the residuals. When you plot where the model over- and under-predicts, the errors form smooth patches — whole regions biased high, others biased low — instead of the salt-and-pepper pattern you want. That structure means the model has failed to explain a systematic component of the target, and it corrupts everything downstream: feature importances get inflated for any covariate that happens to correlate with location, and error bars computed under an independence assumption are far too narrow.

Trend-surface analysis attacks the problem head-on. Instead of hoping a tree ensemble untangles the regional gradient from local effects, you model the gradient explicitly with a smooth function of coordinates, remove it, and let the main model work on what remains. This is complementary to the leakage controls in reducing spatial leakage in model training: detrending removes the signal that causes autocorrelation, while spatial splitting stops that autocorrelation from inflating your scores.

How Trend-Surface Detrending Works

A trend surface is an ordinary polynomial regression whose only predictors are the coordinates themselves. A first-order surface is a tilted plane,

z = b0 + b1·x + b2·y,

which captures a single directional ramp. A second-order surface adds three curvature terms,

z = b0 + b1·x + b2·y + b3·x² + b4·y² + b5·x·y,

so it can bend into domes, basins, and saddles. The fitted surface is the low-frequency backbone of the field; the residual z − ẑ is the high-frequency signal your main model should learn.

The workflow is a clean three-step sandwich:

  1. Fit the polynomial on (x, y) to obtain the trend surface.
  2. Transform by subtracting the surface, producing detrended residuals that you feed to the downstream estimator alongside your real features.
  3. Inverse-transform at inference by evaluating the same surface at the prediction coordinates and adding it back to the model’s residual prediction.

Because the trend is additive, predictions decompose cleanly into regional trend + local residual, which also makes results easy to explain to stakeholders. The diagram below shows the full round trip.

Geographic detrending round trip: fit, model residuals, add trend back The raw target is decomposed into a smooth trend surface fitted on coordinates and a residual field. The main model learns the residuals; at prediction the trend surface is evaluated and added back to reconstruct the target. Raw target z Decompose Model & recombine z(x, y) gradient + local signal Trend surface poly(x, y) → ẑ Residual r = z − ẑ Main model learns r̂ Prediction ŷ = ẑ + r̂ trend added back

Core Principles for Reliable Detrending

  • Project first. The polynomial measures distance through its coordinate terms, so fit it in a metric CRS. Angular latitude/longitude warps the surface; see CRS Alignment and Projection Handling for the reprojection mechanics.
  • Centre and scale the coordinates. Raw UTM eastings run into the hundreds of thousands of metres. Squaring them for a second-order term produces values near 10¹¹ and wrecks the conditioning of the design matrix. Subtract the mean and divide by a length scale before building polynomial features.
  • Stay low-order. First order is the default, second order when a plane leaves visible curvature. Fourth order and beyond oscillate near the data boundary and start eating the local signal you are trying to preserve.
  • Fit the trend on training data only. Estimate surface coefficients inside each cross-validation fold, never on the full dataset, or trend information leaks from validation points into the fit.
  • Verify with Moran’s I. Detrending has succeeded when global Moran’s I on the residuals drops toward zero. This is the objective stopping rule for choosing polynomial order.
  • Keep the surface for inference. Persist the fitted coefficients and the coordinate scaler; you must evaluate the identical surface at prediction time to add the trend back.

Production-Ready SpatialDetrender

The class below implements the fit/transform/inverse_transform contract so the detrender drops into an sklearn pipeline or a TransformedTargetRegressor. It validates inputs, centres and scales the coordinates internally, and stores everything needed to reconstruct the surface at inference.

from __future__ import annotations

import numpy as np
from numpy.typing import NDArray
from sklearn.base import BaseEstimator, RegressorMixin
from sklearn.linear_model import LinearRegression
from sklearn.preprocessing import PolynomialFeatures


class SpatialDetrender(BaseEstimator, RegressorMixin):
    """Remove a low-order polynomial spatial trend from a target field.

    Fits a trend surface z ~ poly(x, y) on projected coordinates. ``transform``
    returns detrended residuals; ``inverse_transform`` adds the surface back.

    Parameters
    ----------
    degree : int, default=1
        Polynomial order of the trend surface. Use 1 (plane) or 2 (curved).
        Higher orders overfit the regional gradient and extrapolate poorly.
    """

    def __init__(self, degree: int = 1) -> None:
        self.degree = degree

    def _check_coords(self, coords: NDArray[np.float64]) -> NDArray[np.float64]:
        coords = np.asarray(coords, dtype=float)
        if coords.ndim != 2 or coords.shape[1] != 2:
            raise ValueError(
                "coords must have shape (n_samples, 2): projected easting/northing."
            )
        if not np.isfinite(coords).all():
            raise ValueError("coords contains NaN or infinite values.")
        return coords

    def fit(
        self,
        coords: NDArray[np.float64],
        z: NDArray[np.float64],
    ) -> "SpatialDetrender":
        coords = self._check_coords(coords)
        z = np.asarray(z, dtype=float).ravel()
        if z.shape[0] != coords.shape[0]:
            raise ValueError("coords and z must have the same number of rows.")
        if self.degree < 1:
            raise ValueError("degree must be >= 1.")

        # Centre and scale coordinates to keep the design matrix well conditioned.
        self.center_ = coords.mean(axis=0)
        self.scale_ = coords.std(axis=0)
        self.scale_[self.scale_ == 0] = 1.0
        scaled = (coords - self.center_) / self.scale_

        self.poly_ = PolynomialFeatures(degree=self.degree, include_bias=True)
        design = self.poly_.fit_transform(scaled)
        self.surface_ = LinearRegression(fit_intercept=False).fit(design, z)
        return self

    def _surface(self, coords: NDArray[np.float64]) -> NDArray[np.float64]:
        scaled = (coords - self.center_) / self.scale_
        return self.surface_.predict(self.poly_.transform(scaled))

    def transform(
        self,
        coords: NDArray[np.float64],
        z: NDArray[np.float64],
    ) -> NDArray[np.float64]:
        """Return detrended residuals z - trend(coords)."""
        coords = self._check_coords(coords)
        z = np.asarray(z, dtype=float).ravel()
        return z - self._surface(coords)

    def inverse_transform(
        self,
        coords: NDArray[np.float64],
        residual: NDArray[np.float64],
    ) -> NDArray[np.float64]:
        """Add the trend surface back to residual predictions."""
        coords = self._check_coords(coords)
        residual = np.asarray(residual, dtype=float).ravel()
        return residual + self._surface(coords)

Step-by-Step Walkthrough

1. Load samples and extract projected coordinates

import geopandas as gpd
import numpy as np

gdf = gpd.read_file("soil_samples.geojson").to_crs("EPSG:32633")  # UTM 33N
coords = np.column_stack([gdf.geometry.x, gdf.geometry.y])
z = gdf["organic_carbon"].to_numpy()

feature_cols = ["ndvi", "slope", "clay_pct", "twi"]
X_features = gdf[feature_cols].to_numpy()

2. Fit the detrender and inspect the residuals

detrender = SpatialDetrender(degree=1).fit(coords, z)
residual = detrender.transform(coords, z)

trend = detrender.inverse_transform(coords, np.zeros_like(z))  # surface alone
trend_r2 = 1.0 - residual.var() / z.var()

print(f"raw target std: {z.std():.3f}")
print(f"residual   std: {residual.std():.3f}")   # smaller — gradient removed
print(f"trend R²       : {trend_r2:.3f}")

A first-order trend that explains 20–40 % of the target’s variance is typical for a field with a genuine regional gradient. If it explains almost nothing, detrending will not help and you should leave the target as-is.

3. Train the main model on residuals

from xgboost import XGBRegressor

model = XGBRegressor(n_estimators=400, max_depth=5, learning_rate=0.05)
model.fit(X_features, residual)   # learns local structure, not geography

4. Predict and add the trend back

residual_pred = model.predict(X_features)
z_pred = detrender.inverse_transform(coords, residual_pred)

Because the surface is additive and deterministic, the reconstruction at step 4 is exact for the trend component — only the residual carries model error. When you evaluate, always score z_pred against the original z, never the residual against itself, and do it under a spatially aware split from Spatial Cross-Validation Strategies so the residual model is not rewarded for leftover geography.

Verification: Moran’s I Should Drop

The objective test of a good detrend is that spatial autocorrelation in the residuals collapses. Compute global Moran’s I on the target before detrending and on the residuals after; a large positive value falling toward zero confirms the surface absorbed the regional gradient.

import numpy as np
from libpysal.weights import KNN
from esda.moran import Moran

w = KNN.from_array(coords, k=8)
w.transform = "r"

moran_raw = Moran(z, w, permutations=999)
moran_res = Moran(residual, w, permutations=999)

print(f"Moran's I (raw target): {moran_raw.I:.3f}  p={moran_raw.p_sim:.3f}")
print(f"Moran's I (residual)  : {moran_res.I:.3f}  p={moran_res.p_sim:.3f}")

A typical result moves from I ≈ 0.62 on the raw target to I ≈ 0.11 on the residuals — still slightly positive, because a smooth surface cannot capture short-range structure, but no longer a dominant regional ramp. If Moran’s I on the residuals is still high, step up to degree=2; if bumping the degree barely moves it, the remaining autocorrelation is local and belongs to your main model, not the trend. The full interpretation of these statistics is covered in Diagnosing Spatial Autocorrelation with Global Moran’s I, and the neighbourhood weights you build here are the same ones used for Spatial Lag and Neighborhood Statistics.

Watch the failure mode in the opposite direction: if raising the polynomial order keeps driving Moran’s I down but your cross-validated error stops improving or gets worse, the surface has started chasing noise rather than the true gradient. Stop at the lowest order that removes the broad ramp and leaves a residual field with no large coherent patches.

FAQ

What is the difference between a first-order and second-order trend surface?

A first-order trend surface is a plane, z = b0 + b1·x + b2·y. It captures a single directional gradient, such as values rising steadily from south to north. A second-order surface adds the quadratic and interaction terms x², y², and x·y, letting the surface bend to model doming, basins, or saddle-shaped regional patterns. First-order is the safe default; move to second-order only when a linear plane leaves obvious curvature in the residuals.

Why should I use projected coordinates instead of latitude and longitude for detrending?

A polynomial trend surface measures distance implicitly through its coordinate terms. Raw latitude and longitude are angular, so one unit of longitude represents a different ground distance at every latitude. Fitting the polynomial on distorted axes warps the surface and biases the residuals. Reproject to a metric CRS such as UTM or an equal-area projection first so both axes are in metres and the fitted gradient is geometrically honest.

How high an order polynomial should I use for the trend surface?

Keep the order as low as possible, usually first or second order. A trend surface is meant to strip broad regional gradients, not local structure — that is the job of your main model. High-order polynomials oscillate wildly near the edges of the data, absorb local signal you want to keep, and extrapolate catastrophically outside the sampled extent. Validate the choice by checking that Moran’s I on the residuals drops toward zero without the surface starting to chase noise.

Part of: Handling Spatial Autocorrelation — Training Geospatial Predictive Models in Python