Stationarity & Trend Analysis in Spatial Modeling

Stationarity & Trend Analysis is the diagnostic gateway between raw spatial data and reliable geostatistical inference. Within the broader Core Concepts of Spatial Statistics & Geostatistics framework, it answers a deceptively simple question: is this process homogeneous enough to model with a single variogram? When environmental gradients, urban heat islands, or systematic land-use transitions dominate a study area, ignoring non-stationarity inflates kriging variance, biases regression coefficients, and produces uncertainty bounds that are wide in the wrong places. The five-step workflow below gives you the Python tools to diagnose, quantify, and remove trend before fitting any covariance model.

Prerequisites #

• Python 3.10+
• geopandas >= 0.14, numpy >= 1.26, scipy >= 1.12, statsmodels >= 0.14, scikit-gstat >= 1.0, pyproj >= 3.6
• Data projected into a locally appropriate distance-preserving CRS (UTM, national equal-area, or State Plane) — not geographic lat/lon
• No exact-coordinate duplicates; missing values resolved via spatially aware imputation or complete-case filtering

Install in one step:

# pip install geopandas>=0.14 numpy>=1.26 scipy>=1.12 statsmodels>=0.14 scikit-gstat>=1.0 pyproj>=3.6

Mathematical Core #

Geostatistical stationarity operates at two nested levels.

First-order stationarity requires that the expected value of the spatial random function $Z(\mathbf{s})$ is constant everywhere:

$$\mathrm{E}[Z(\mathbf{s})] = \mu \quad \forall\,\mathbf{s} \in \mathcal{D}$$

Second-order (weak) stationarity additionally constrains the covariance to depend only on the separation vector $\mathbf{h}$, not on the absolute location:

$$\mathrm{Cov}[Z(\mathbf{s}),\, Z(\mathbf{s}+\mathbf{h})] = C(\mathbf{h}) \quad \forall\,\mathbf{s},\,\mathbf{h}$$

This implies constant variance $\sigma^2 = C(\mathbf{0})$ and a well-defined, bounded semivariogram:

$$\gamma(\mathbf{h}) = \frac{1}{2}\,\mathrm{Var}[Z(\mathbf{s}+\mathbf{h}) - Z(\mathbf{s})] = C(\mathbf{0}) - C(\mathbf{h})$$

When second-order stationarity fails but spatial increments remain stationary (unbounded variance), the intrinsic hypothesis holds and ordinary kriging remains valid. If even increments are non-stationary — recognisable by a variogram that grows parabolically with no sill — a deterministic drift $m(\mathbf{s})$ must be modelled explicitly:

$$Z(\mathbf{s}) = m(\mathbf{s}) + \varepsilon(\mathbf{s})$$

where $m(\mathbf{s}) = \sum_{k=0}^{K} a_k f_k(\mathbf{s})$ is typically a polynomial in the coordinates and $\varepsilon(\mathbf{s})$ is the residual random field assumed second-order stationary.

Every symbol in that decomposition has a physical reading. $Z(\mathbf{s})$ is the single realisation you actually measured — one draw from an ensemble that can never be resampled, which is precisely why a stationarity assumption is needed at all: it is what licenses pairs separated by the same $\mathbf{h}$ in different parts of the domain to stand in for repeated draws at one location. $m(\mathbf{s})$ is the portion of the field you are willing to call deterministic and explain with coordinates or covariates; $\varepsilon(\mathbf{s})$ is everything left over, and the variogram has to describe it unaided. The basis functions $f_k(\mathbf{s})$ are the drift terms — $1, x, y$ for a first-order surface, plus $x^2, xy, y^2$ for a second-order one — while the coefficients $a_k$ are what least squares estimates. The lag $\mathbf{h}$ is a vector, not a scalar: retaining its direction is what makes directional variograms possible, and collapsing it to its length $h = \|\mathbf{h}\|$ is already an isotropy assumption smuggled in before any model is fitted.

The three hypotheses nest. Strict stationarity requires the entire joint distribution to be translation-invariant; second-order stationarity asks only that the first two moments be, which is all that kriging actually consumes; the intrinsic hypothesis relaxes further still and asks only that increments $Z(\mathbf{s}+\mathbf{h}) - Z(\mathbf{s})$ have constant mean and variance. Brownian motion is the canonical process that is intrinsic but not second-order stationary: its variance grows without bound, so $C(\mathbf{h})$ does not exist, yet $\gamma(h) = h/2$ remains perfectly well defined. That asymmetry is why the semivariogram rather than the covariance function is the primitive object in geostatistics — it survives assumptions the covariance does not.

Stationarity is also a property of the model at a chosen scale, not an intrinsic property of the ground. A soil-metal field that is plainly non-stationary across a 40 km catchment can often be treated as stationary inside a 2 km moving-window neighbourhood, which is exactly what a local search radius in ordinary kriging buys you. Before committing to a global drift surface, ask whether shrinking the search neighbourhood achieves the same result more cheaply and with fewer parameters to defend.

Step 1 — Coordinate System Validation & Topology Cleaning #

Stationarity tests are inherently coordinate-system-dependent. Geographic coordinates (latitude/longitude) introduce severe metric distortion at non-equatorial latitudes, compressing lag distances near the poles and stretching them at the equator. Always transform data to a locally appropriate projected CRS before computing any distances:

import geopandas as gpd
import numpy as np

# Load data and assert projected CRS
gdf = gpd.read_file("spatial_data.gpkg")
assert gdf.crs is not None and gdf.crs.is_projected, (
    "Data must be in a projected CRS for distance calculations. "
    f"Current CRS: {gdf.crs}"
)

# Remove exact-coordinate duplicates (masked by overlapping point symbols)
gdf = gdf.drop_duplicates(subset=["geometry"])

# Drop rows with missing target values
gdf = gdf.dropna(subset=["z"])

# Extract coordinates as a plain numpy array for downstream routines
coords = np.column_stack([gdf.geometry.x, gdf.geometry.y])  # shape (n, 2)
values = gdf["z"].to_numpy()
print(f"{len(gdf)} observations in CRS: {gdf.crs.to_epsg()}")

Uneven sampling density can masquerade as a spatial trend: dense high-value points concentrated in one corner of the domain will deflect moving averages even when the true process is stationary. Address this with sampling bias mitigation strategies — declustering or kernel density correction — before moving to variogram estimation.

Step 2 — Exploratory Gradient Detection #

Visual diagnostics should precede any formal test. Generate three complementary plots:

import matplotlib.pyplot as plt

fig, axes = plt.subplots(1, 2, figsize=(11, 4))

# Coordinate-vs-value scatter: systematic slope = first-order drift
axes[0].scatter(coords[:, 0], values, s=8, alpha=0.5)
axes[0].set(xlabel="Easting (m)", ylabel="z", title="z vs Easting")

axes[1].scatter(coords[:, 1], values, s=8, alpha=0.5)
axes[1].set(xlabel="Northing (m)", ylabel="z", title="z vs Northing")

plt.tight_layout()
plt.savefig("gradient_scatter.png", dpi=150)

A consistent positive or negative slope in either panel is strong evidence of first-order drift. Next, compute directional variograms (N–S and E–W) and compare them at short lags:

import skgstat as skg

# Omnidirectional variogram as baseline
V_omni = skg.Variogram(
    coords, values,
    model="spherical",
    n_lags=15,
    maxlag=0.6,  # 60% of max pair distance
    estimator="matheron",
)

# Directional variograms — azimuth 0° (N-S) and 90° (E-W)
V_ns = skg.DirectionalVariogram(coords, values, azimuth=0,   tolerance=22.5, n_lags=12)
V_ew = skg.DirectionalVariogram(coords, values, azimuth=90,  tolerance=22.5, n_lags=12)

If the N–S and E–W experimental variograms diverge strongly at short lags (rather than at medium lags where geometric anisotropy is expected), the divergence is likely caused by directional drift, not structural anisotropy. Always detrend before interpreting anisotropy ratios.

Step 3 — Formal Stationarity Testing #

Exploratory plots provide intuition; production pipelines need quantitative evidence. Two complementary tests cover the two stationarity levels:

Moran’s I on Raw Values (First-Order Check) #

Spatial autocorrelation metrics such as Moran’s I measure the degree to which nearby values are more similar than expected under spatial randomness. A large, positive Moran’s I on the raw field — especially when paired with a gradient in the coordinate-scatter plots — indicates unmodeled drift:

from libpysal.weights import DistanceBand
from esda.moran import Moran

# Build distance-band spatial weights (threshold = 10 km here)
w = DistanceBand.from_array(coords, threshold=10_000, binary=True)
w.transform = "r"   # row-standardise

mi_raw = Moran(values, w, permutations=999)
print(f"Moran's I (raw): {mi_raw.I:.4f}  p-value: {mi_raw.p_sim:.4f}")
# Expected: large positive I with p << 0.05 when drift is present

For the implementation details of Moran’s I, see the dedicated guide on how to calculate Moran’s I in PySAL.

Levene’s Test for Variance Homogeneity (Second-Order Check) #

Split the domain into distance bins and test whether variance is constant across bins. Significantly different bin variances indicate second-order non-stationarity:

from scipy.stats import levene

# Assign observations to 4 quadrants by coordinate median
x_med = np.median(coords[:, 0])
y_med = np.median(coords[:, 1])
mask_nw = (coords[:, 0] < x_med) & (coords[:, 1] >= y_med)
mask_ne = (coords[:, 0] >= x_med) & (coords[:, 1] >= y_med)
mask_sw = (coords[:, 0] < x_med) & (coords[:, 1] < y_med)
mask_se = (coords[:, 0] >= x_med) & (coords[:, 1] < y_med)

stat, p = levene(
    values[mask_nw], values[mask_ne],
    values[mask_sw], values[mask_se],
    center="median",  # robust to skew
)
print(f"Levene's test: W={stat:.3f}  p={p:.4f}")
# p < 0.05 → heterogeneous variance → second-order stationarity suspect

Both tests degrade in ways worth knowing before you trust a p-value. Moran’s I is a function of the weights matrix as much as of the data: widening threshold from 10 km to 30 km averages the drift over more neighbours and pulls $I$ toward zero, so a threshold tuned until the test passes is not a test at all. Report $I$ at two or three thresholds bracketing the expected correlation range and confirm that the sign and the significance are stable across them. Levene’s test on quadrants has the opposite weakness. It is a variance test, but quadrant means differ under first-order drift, and center="median" removes only the within-group centre — so a strong linear gradient inflates each quadrant’s spread by a different amount and Levene’s rejects for a reason that has nothing to do with second-order stationarity. Run it on detrended residuals as well as on raw values: if it is significant on the raw field and non-significant on the residuals, the apparent heteroscedasticity was drift in disguise and no further action is needed. Genuine second-order non-stationarity — a field whose variance really does grow toward one edge, as concentration data often does where the mean is high — survives detrending, and it calls for a variance-stabilising transform or a locally varying sill rather than for more drift terms.

Step 4 — Trend Detrending & Residual Extraction #

Once non-stationarity is confirmed, decompose the field into deterministic drift and stochastic residual. The right detrending strategy depends on the complexity of the gradient:

The polynomial approach is transparent and easy to validate:

from numpy.linalg import lstsq

# Build design matrix: intercept + x + y + x² + xy + y²
x = coords[:, 0]
y = coords[:, 1]
x_sc = (x - x.mean()) / x.std()   # centre & scale for numerical stability
y_sc = (y - y.mean()) / y.std()

A = np.column_stack([
    np.ones(len(x_sc)),    # β₀ — intercept
    x_sc,                  # β₁ — linear east-west gradient
    y_sc,                  # β₂ — linear north-south gradient
    x_sc**2,               # β₃ — quadratic easting
    x_sc * y_sc,           # β₄ — interaction term
    y_sc**2,               # β₅ — quadratic northing
])

beta, residuals, rank, sv = lstsq(A, values, rcond=None)
trend = A @ beta
residuals_z = values - trend

print(f"Trend R²: {1 - np.var(residuals_z)/np.var(values):.4f}")
# Aim for R² that captures the large-scale drift without over-fitting

Use AIC/BIC to select polynomial order: add terms only while the information criterion decreases. Over-fitting with high-order polynomials absorbs genuine spatial correlation into the drift estimate, leaving residuals with artificially low variance and inflated kriging uncertainty.

Step 5 — Residual Covariance Modelling & Validation #

After detrending, verify that residuals satisfy second-order stationarity before fitting the variogram:

# Re-check Moran's I on residuals — should be near zero
mi_res = Moran(residuals_z, w, permutations=999)
print(f"Moran's I (residuals): {mi_res.I:.4f}  p-value: {mi_res.p_sim:.4f}")
# Target: I close to 0, p >> 0.05

# Fit variogram to residuals
V_res = skg.Variogram(
    coords, residuals_z,
    model="spherical",
    n_lags=15,
    maxlag=0.5,
    estimator="cressie",   # robust Cressie-Hawkins estimator
)
print(V_res)
# Inspect: sill should stabilise; nugget/sill ratio < 0.3 is typical

One caveat undercuts the whole two-stage recipe and is worth stating plainly: the variogram of OLS residuals is biased low, and increasingly so at long lags. Least squares picks $\hat{a}_k$ to minimise the residual sum of squares under an implicit independence assumption, so it absorbs exactly the low-frequency variation the residual variogram is supposed to detect. The symptom is a residual sill and range that both come out smaller than the truth — typically by 10–30% for a second-order drift fitted over a domain only a few correlation ranges wide, and worse as the ratio of range to domain size grows. The consequence for prediction is mild, because kriging weights are fairly insensitive to a proportional error in the sill; the consequence for uncertainty is not, because a shrunken sill produces prediction intervals that are too narrow everywhere, including in the sparsely sampled corners where they matter most. The fix is to stop treating detrending and covariance estimation as independent steps. Either estimate both jointly by restricted maximum likelihood — wrapping scipy.optimize.minimize around the profile REML objective, since scikit-gstat fits the variogram alone — or iterate: fit the drift by OLS, estimate the residual variogram, refit the drift by generalised least squares using that covariance, re-estimate, and stop when the fitted range moves less than a few percent. Two or three passes is usually enough, and the range almost always grows on the second pass, which is the bias correcting itself.

For the full Python implementation including likelihood ratio tests and reproducible notebook output, see the step-by-step guide on Testing for Second-Order Stationarity in Python.

Output Interpretation #

A successful detrending run produces four observable signals:

1. Residual scatter vs coordinates: flat, structureless cloud centred on zero — no remaining slope.
2. Residual variogram: sill stabilises before $h = 0.5 \times$ maximum pair distance; nugget-to-sill ratio below 0.35.
3. Moran’s I on residuals: $I < 0.05$ with $p > 0.05$ (non-significant).
4. LOOCV standardised errors: mean $\approx 0$, standard deviation $\approx 1$ — correctly calibrated uncertainty.

Warning signs:

• Residual variogram still growing parabolically: increase polynomial order or switch to a GAM.
• Moran’s I significant after detrending: the trend model misses a regional sub-domain; consider local detrending or moving-window kriging.
• LOOCV errors systematically positive in one quadrant: the drift function has a sign error or the CRS transformation introduced a coordinate-axis swap.

Production Considerations #

Computational scaling: Polynomial surface fitting via least squares is $O(n \cdot K^2)$ — fast for any realistic sample size. Variogram estimation scales roughly $O(n^2)$ for exhaustive pair enumeration; use scikit-gstat’s random subsampling option (subsample=5000) for $n > 20{,}000$.

Chunked processing for large grids: When the prediction grid is large ($> 10^6$ cells), rasterise the fitted trend surface separately from the kriging residuals and add them back after prediction. This avoids loading the full prediction matrix into memory.

Reproducibility: Cache fitted polynomial coefficients (beta) and the residual variogram parameters (nugget, sill, range) as JSON alongside the model artefacts. Re-applying the same trend to a new batch of observations requires only the stored beta vector and the coordinate scaling parameters (x.mean(), x.std(), etc.).

Parallel variogram estimation: For directional variograms on large datasets, use dask to parallelise pair computation across azimuth bands:

# Pseudo-code for dask-parallelised directional variogram
import dask
tasks = [
    dask.delayed(skg.DirectionalVariogram)(coords, residuals_z, azimuth=az, n_lags=12)
    for az in [0, 45, 90, 135]
]
variograms = dask.compute(*tasks)

Troubleshooting #

Next Steps #

For a full implementation with likelihood ratio model comparison and reproducible notebook output, follow the step-by-step guide on Testing for Second-Order Stationarity in Python. Once residuals are confirmed stationary, proceed to spatial weight matrices to encode the neighbourhood structure for regression, or move directly to ordinary and universal kriging to generate interpolated prediction surfaces from the detrended residuals.

← Back to Core Concepts of Spatial Statistics & Geostatistics

Related #

• Testing for Second-Order Stationarity in Python
• Spatial Autocorrelation Metrics
• Spatial Weight Matrices
• Ordinary & Universal Kriging
• Sampling Bias Mitigation