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.

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

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

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