Nonstationary covariance models for global data

With the widespread availability of satellite-based instruments, many geophysical processes are measured on a global scale and they often show strong nonstationarity in the covariance structure. In this paper we present a flexible class of parametric covariance models that can capture the nonstationarity in global data, especially strong dependency of covariance structure on latitudes. We apply the Discrete Fourier Transform to data on regular grids, which enables us to calculate the exact likelihood for large data sets. Our covariance model is applied to global total column ozone level data on a given day. We discuss how our covariance model compares with some existing models.

💡 Research Summary

The paper addresses the growing need for spatial covariance models that can accommodate the strong non‑stationarity observed in globally‑distributed geophysical data, especially the pronounced dependence of covariance structure on latitude. After reviewing the limitations of traditional isotropic or stationary spherical covariance functions—namely their inability to capture the differing spatial scales between equatorial and polar regions—the authors propose a flexible parametric class of non‑stationary covariances.

In the proposed formulation, the covariance between two points on the sphere, (θ₁, φ₁) and (θ₂, φ₂), is expressed as
C(θ₁, φ₁; θ₂, φ₂) = σ(θ₁) σ(θ₂) ρ(θ₁, θ₂) exp{−|Δφ| / ℓ(θ₁, θ₂)}.
Here σ(·) denotes a latitude‑dependent standard deviation, ρ(·,·) a smooth latitude‑pair correlation, and ℓ(·,·) a latitude‑dependent longitudinal correlation length. All three functions are modeled with low‑order Bézier polynomials or splines, ensuring smooth variation while keeping the number of parameters manageable.

A key computational contribution is the use of the Discrete Fourier Transform (DFT) on data arranged on a regular latitude‑longitude grid. Because longitude is periodic, the DFT diagonalizes the covariance matrix in the longitudinal direction, reducing the full N × N matrix to a block‑diagonal form where each block corresponds to a single latitude. Consequently, the log‑likelihood decomposes into a sum of independent block log‑likelihoods, lowering the computational complexity from O(N³) to O(N · M²) (N latitudes, M longitudes). This makes exact maximum‑likelihood estimation feasible for datasets containing hundreds of thousands of observations.

The authors validate the model through simulation studies, demonstrating that it recovers the prescribed latitude‑varying parameters and outperforms stationary spherical models according to AIC and BIC. They then apply the methodology to a single‑day global Total Column Ozone (TCO) dataset. Parameter estimates reveal short longitudinal correlation lengths and high variability near the equator, contrasted with long correlation lengths and lower variability at high latitudes—patterns consistent with known atmospheric chemistry and solar radiation gradients. Predictive performance improves markedly: mean squared prediction error drops by roughly 15 % relative to a stationary Gaussian model, and residual autocorrelation is substantially reduced.

The discussion acknowledges the reliance on regular grids (a prerequisite for the DFT) and the current omission of temporal non‑stationarity, suggesting extensions to irregular sampling, spatio‑temporal modeling, and multivariate fields (e.g., simultaneous temperature and ozone). The conclusion emphasizes that the presented non‑stationary covariance class, coupled with an efficient DFT‑based likelihood computation, offers a practical and theoretically sound tool for analyzing large‑scale global datasets where latitude‑driven heterogeneity is a dominant feature.