Ensemble Spatial Interpolation Methods

• Ensemble spatial interpolation is a method that constructs multiple weak predictors from randomized spatial partitions and local interpolators to form an aggregated estimate.
• It integrates techniques like inverse distance weighting, kriging, and deep neural networks to compute point estimates and quantify prediction uncertainty.
• The approach is applied in meteorology, mining, and epidemic tracking, offering improved accuracy and efficiency, especially for large-scale spatial datasets.

Ensemble spatial interpolation is a formulation in which one observes a spatial process $z(x)$ at locations $\mathcal P=\{x_j\}_{j=1}^{N_s}$ with values $\mathcal M=\{z_j\}$, constructs multiple spatial predictors, and aggregates them into a final estimate. In Spatialize v1.0, this is written as $\widehat z(x)=G(f_1(x),\dots,f_m(x))$, where the $f_k$ are “weak voters” produced from stochastic spatial partitions and local interpolators; the ensemble outputs can also be interpreted as an empirical posterior predictive distribution (Egaña et al., 23 Jul 2025). Related constructions aggregate bootstrap deep neural networks for bivariate wind fields, weighted off-site neural submodels for frost prediction, and quantile-regression base learners for probabilistic precipitation correction (Nag et al., 2023, Zhou et al., 2022, Papacharalampous et al., 2024).

1. Core estimator and mathematical structure

In the partition-based ESI formulation, the domain is partitioned repeatedly, and each partition induces a local interpolation problem. Once a tree $T_k$ is fixed, a target location $x$ is associated with the subset $\mathcal L_k\subset\mathcal P$ lying in the same leaf cell, and the $k$-th weak voter is

$f_k(x)\equiv S_{\mathcal L_k}(x),$

where $\mathcal P=\{x_j\}_{j=1}^{N_s}$0 is a local interpolator such as IDW or ordinary kriging. The final estimate is then obtained by an aggregation operator $\mathcal P=\{x_j\}_{j=1}^{N_s}$1; the simplest case is the ensemble mean,

$\mathcal P=\{x_j\}_{j=1}^{N_s}$2

Posterior-style uncertainty is derived from the empirical ensemble, with posterior variance

$\mathcal P=\{x_j\}_{j=1}^{N_s}$3

and posterior $\mathcal P=\{x_j\}_{j=1}^{N_s}$4-quantile $\mathcal P=\{x_j\}_{j=1}^{N_s}$5 (Egaña et al., 23 Jul 2025).

The local interpolator can be purely distance-based or geostatistical. For IDW,

$\mathcal P=\{x_j\}_{j=1}^{N_s}$6

with user-supplied exponent $\mathcal P=\{x_j\}_{j=1}^{N_s}$7. For ordinary kriging, a local variogram model $\mathcal P=\{x_j\}_{j=1}^{N_s}$8 is assumed, with nugget $\mathcal P=\{x_j\}_{j=1}^{N_s}$9, range $\mathcal M=\{z_j\}$0, sill $\mathcal M=\{z_j\}$1, and a model such as spherical or exponential (Egaña et al., 23 Jul 2025).

A distinct neural formulation appears in bivariate DeepKriging. There, a spatial location $\mathcal M=\{z_j\}$2 is embedded through radial basis features

$\mathcal M=\{z_j\}$3

and a deep neural network maps $\mathcal M=\{z_j\}$4 to the two components of the wind field,

$\mathcal M=\{z_j\}$5

In practice, $\mathcal M=\{z_j\}$6 is a fully connected feed-forward net with hidden layers

$\mathcal M=\{z_j\}$7

followed by a linear output layer (Nag et al., 2023).

2. Modes of ensemble construction

The literature contains several distinct mechanisms for generating the ensemble members that are ultimately aggregated. Some methods randomize spatial partitions, some randomize training data, and some combine heterogeneous predictive algorithms.

Paradigm	Ensemble elements	Aggregation / output
Partition-based ESI (Egaña et al., 23 Jul 2025)	Mondrian Forest or Voronoi Forest partitions with local IDW or kriging	mean, median, percentiles, empirical posterior
Bootstrap Bivariate DeepKriging (Nag et al., 2023)	$\mathcal M=\{z_j\}$8 bootstrap replicates; master fit $\mathcal M=\{z_j\}$9; re-train last $\widehat z(x)=G(f_1(x),\dots,f_m(x))$0 layers	ensemble mean $\widehat z(x)=G(f_1(x),\dots,f_m(x))$1 and variance $\widehat z(x)=G(f_1(x),\dots,f_m(x))$2
Off-site ANN frost interpolation (Zhou et al., 2022)	60 “weak learner” ANNs built from source stations	simple average, weighted average, weighted voting
Quantile-regression ensembles (Papacharalampous et al., 2024)	QR, QRF, GRF, GBM, LightGBM, QRNN	mean, median, best combiner, stacking

In Spatialize, two partition generators are explicit. A Mondrian Forest samples random axis-aligned trees over the bounding box $\widehat z(x)=G(f_1(x),\dots,f_m(x))$3, while a Voronoi Forest samples Voronoi partitions by drawing $\widehat z(x)=G(f_1(x),\dots,f_m(x))$4 nuclei and assigning each $\widehat z(x)=G(f_1(x),\dots,f_m(x))$5 to its nearest nucleus. When conditioning on data, each sub-box $\widehat z(x)=G(f_1(x),\dots,f_m(x))$6 is replaced by the minimal box $\widehat z(x)=G(f_1(x),\dots,f_m(x))$7 containing all sample points in $\widehat z(x)=G(f_1(x),\dots,f_m(x))$8, which guarantees each leaf has at least one sample (Egaña et al., 23 Jul 2025).

In bootstrap Bivariate DeepKriging, the original training set $\widehat z(x)=G(f_1(x),\dots,f_m(x))$9 is resampled $f_k$0 times. For each bootstrap replicate, the network is initialized at the weights $f_k$1 from a single master fit, the first $f_k$2 layers are frozen, and only the last $f_k$3 layers are re-trained. The ensemble predictions at location $f_k$4 are $f_k$5, with ensemble mean

$f_k$6

and ensemble variance

$f_k$7

(Nag et al., 2023).

The frost-prediction framework uses a different decomposition. Instead of a single interpolator, it trains one off-site ANN per source station. Each of the 60 trained submodels yields a temperature estimate $f_k$8 for a target site $f_k$9, and the outputs are fused either by the simple average, a weighted average based on normalized geographic distance, DEM difference, and NDVI difference, or weighted voting for frost-event classification (Zhou et al., 2022).

Quantile-regression ensembles replace identical weak learners by heterogeneous base learners. With $T_k$0 base quantile predictors—QR, QRF, GRF, GBM, LightGBM, and QRNN—the ensemble can be the mean combiner, median combiner, best combiner, or a stacked second-level quantile regressor $T_k$1 trained on the base predictions (Papacharalampous et al., 2024).

An adjacent but structurally different interpolation strategy appears in clustering-based EMOS. Stations are clustered in a feature space using standard $T_k$2-means, one EMOS model is estimated per cluster, and an ungauged location is assigned to the nearest cluster center in Euclidean distance. Here $T_k$3 recovers regional EMOS and $T_k$4 recovers local EMOS (Baran et al., 2024).

3. Predictive distributions and uncertainty quantification

A defining feature of ensemble spatial interpolation is that the ensemble can serve not only as a point predictor but also as a distributional object. In Spatialize, the set $T_k$5 is interpreted as draws from the posterior predictive distribution $T_k$6, from which one computes empirical means, variances, credible intervals, quantiles, or custom loss-based precision measures

$T_k$7

This places uncertainty quantification directly inside the interpolation procedure rather than in a separate post hoc stage (Egaña et al., 23 Jul 2025).

Bivariate DeepKriging develops a distribution-free uncertainty quantification method based on bootstrap and ensemble DNN. For a single component, empirical quantile intervals use the empirical $T_k$8-th quantiles of the bootstrap predictions. A second option is a $T_k$9-based bootstrap interval,

$x$0

where $x$1 is the sample standard deviation of the component-wise bootstrap predictions. A third construction yields uniform confidence bands over a region $x$2 by computing the sup-deviations

$x$3

and the empirical quantile $x$4 of $x$5. For real wind data, prediction-interval coverage for 95% intervals is near nominal for DeepKriging, at 94–95% (Nag et al., 2023).

Quantile-regression ensembles produce predictive distributions through conditional quantiles rather than through an explicit parametric law. For quantile level $x$6, a forecast is evaluated by the strictly proper quantile scoring rule

$x$7

and the mean test score $x$8 yields a quantile skill score against a reference model. Empirical coverage frequency is also used as an evaluation criterion (Papacharalampous et al., 2024).

Operational weather post-processing provides two further distributional templates. In the clustering-based interpolation of EMOS, wind speed is modeled by a left-truncated normal predictive law,

$x$9

with affine link functions for $\mathcal L_k\subset\mathcal P$0 and $\mathcal L_k\subset\mathcal P$1 estimated by minimizing the mean CRPS. At ungauged locations, forecast-based features are used for cluster assignment and parameter transfer, and coverage and average width of central 96.23% prediction intervals are verified (Baran et al., 2024). In spatially adaptive post-processing for temperature, both the short-term temperature average $\mathcal L_k\subset\mathcal P$2 and the local uncertainty predictor $\mathcal L_k\subset\mathcal P$3 are interpolated through intrinsic Gaussian random fields, and the interpolation variance $\mathcal L_k\subset\mathcal P$4 is added to the post-processed forecast variance,

$\mathcal L_k\subset\mathcal P$5

so that uncertainty in the interpolated mean is propagated explicitly (Scheuerer et al., 2013).

4. Computational regimes and implementation

The computational motivation for ensemble spatial interpolation is especially clear when classical covariance-based methods become prohibitive. In bivariate wind interpolation, classical cokriging for $\mathcal L_k\subset\mathcal P$6 training sites requires $\mathcal L_k\subset\mathcal P$7 cost to decompose the $\mathcal L_k\subset\mathcal P$8 covariance matrix and $\mathcal L_k\subset\mathcal P$9 to predict at $k$0 new sites, which becomes prohibitive once $k$1 exceeds a few thousands. DeepKriging training by stochastic gradient descent costs roughly $k$2, scales nearly linearly in $k$3, can be easily parallelized on GPUs, and has prediction cost $k$4 at a single new location, independent of $k$5. Empirically, on real wind data with $k$6, Bivariate DeepKriging is $k$7 faster than exact cokriging, and in smaller simulated experiments up to $k$8 it begins to outperform cokriging in wall-clock time around $k$9 (Nag et al., 2023).

Spatialize emphasizes stochastic partitioning and parallel evaluation. Its core algorithms are written in C++ and compiled as a Python extension via pybind11, while the high-level interface is provided in Python 3.x through esi_griddata(), esi_nongriddata(), and esi_hparams_search(). Each partition $f_k(x)\equiv S_{\mathcal L_k}(x),$0 can be processed independently; the C++ core uses OpenMP or std::thread to evaluate the $f_k(x)\equiv S_{\mathcal L_k}(x),$1 trees in parallel. The framework is described as designed to handle $f_k(x)\equiv S_{\mathcal L_k}(x),$2 up to $f_k(x)\equiv S_{\mathcal L_k}(x),$3–$f_k(x)\equiv S_{\mathcal L_k}(x),$4 on modern multicore machines. Default recommendations include ensemble size $f_k(x)\equiv S_{\mathcal L_k}(x),$5, Mondrian $f_k(x)\equiv S_{\mathcal L_k}(x),$6, Voronoi $f_k(x)\equiv S_{\mathcal L_k}(x),$7, and IDW exponent $f_k(x)\equiv S_{\mathcal L_k}(x),$8 (Egaña et al., 23 Jul 2025).

The post-processing literature supplies additional operational benchmarks. In adaptive EMOS for temperature, rolling-window training requires solving a CRPS minimization in $f_k(x)\equiv S_{\mathcal L_k}(x),$9 time; with $\mathcal P=\{x_j\}_{j=1}^{N_s}$00 and $\mathcal P=\{x_j\}_{j=1}^{N_s}$01, this is described as negligible. Kriging for $\mathcal P=\{x_j\}_{j=1}^{N_s}$02 stations requires solving one linear system twice, once for $\mathcal P=\{x_j\}_{j=1}^{N_s}$03 and once for $\mathcal P=\{x_j\}_{j=1}^{N_s}$04, and in practice $\mathcal P=\{x_j\}_{j=1}^{N_s}$05 leads to an inversion in under a second on a standard CPU (Scheuerer et al., 2013). In EnOSI, the full covariance matrix is never stored explicitly. Instead, the covariance kernel is precomputed or applied by FFT-based convolution, which avoids the cost of a dense $\mathcal P=\{x_j\}_{j=1}^{N_s}$06 covariance in high-dimensional spatial epidemic models; by contrast, a standard EnKF would recompute the sample covariance at each assimilation, which is more expensive (Krishnamurthy et al., 2010).

5. Empirical domains and reported performance

Large-scale wind interpolation provides one of the clearest demonstrations of ensemble spatial interpolation under non-Gaussian and nonstationary behavior. In simulation, under a purely Gaussian stationary field, cokriging with Matérn cores and DeepKriging achieve similar MSPE, approximately $\mathcal P=\{x_j\}_{j=1}^{N_s}$07 versus $\mathcal P=\{x_j\}_{j=1}^{N_s}$08, but DeepKriging is much faster. Under heavy-tailed or non-Gaussian simulated fields, DeepKriging’s MSPE, approximately $\mathcal P=\{x_j\}_{j=1}^{N_s}$09, dramatically beats Matérn cokriging at approximately $\mathcal P=\{x_j\}_{j=1}^{N_s}$10. Under nonstationary settings, DeepKriging attains MSPEs on the order of $\mathcal P=\{x_j\}_{j=1}^{N_s}$11, versus $\mathcal P=\{x_j\}_{j=1}^{N_s}$12–$\mathcal P=\{x_j\}_{j=1}^{N_s}$13 for cokriging variants. On real wind data over the Middle East region at $\mathcal P=\{x_j\}_{j=1}^{N_s}$14 locations, the root mean squared prediction error for the $\mathcal P=\{x_j\}_{j=1}^{N_s}$15 component falls from about $\mathcal P=\{x_j\}_{j=1}^{N_s}$16 for cokriging LMC to about $\mathcal P=\{x_j\}_{j=1}^{N_s}$17 for Bivariate DeepKriging, while similar gains occur for the $\mathcal P=\{x_j\}_{j=1}^{N_s}$18 component, from $\mathcal P=\{x_j\}_{j=1}^{N_s}$19 to $\mathcal P=\{x_j\}_{j=1}^{N_s}$20 (Nag et al., 2023).

Probabilistic precipitation correction gives a second large-scale example. Using monthly gauge totals and satellite precipitation over CONUS, the two best ensemble methods were stacking with a QR combiner and stacking with a QRNN combiner; their quantile skill score improvements over the linear-QR reference ranged from $\mathcal P=\{x_j\}_{j=1}^{N_s}$21 to $\mathcal P=\{x_j\}_{j=1}^{N_s}$22 across quantile levels from $\mathcal P=\{x_j\}_{j=1}^{N_s}$23 to $\mathcal P=\{x_j\}_{j=1}^{N_s}$24 (Papacharalampous et al., 2024). Frost prediction with limited local data provides a contrasting short-range application: the weighted-average ensemble of off-site ANNs consistently outperformed simple averaging and matched or exceeded IDW and ordinary kriging in RMSE, while attaining a detection rate up to $\mathcal P=\{x_j\}_{j=1}^{N_s}$25 (Zhou et al., 2022).

Operational weather post-processing highlights the interpolation of predictive distributions rather than scalar fields alone. In the EMOS case study based on ECMWF wind speed forecasts, forecast-based semi-local EMOS using raw forecasts as features achieved the best compromise at unobserved locations, with CRPS improvements of $\mathcal P=\{x_j\}_{j=1}^{N_s}$26–$\mathcal P=\{x_j\}_{j=1}^{N_s}$27 over regional EMOS, sharper 96.23% intervals, coverage close to nominal, and nearly uniform PIT histograms (Baran et al., 2024). In temperature post-processing with COSMO-DE-EPS, adaptive EMOS achieved out-of-sample CRPS $\mathcal P=\{x_j\}_{j=1}^{N_s}$28, compared with $\mathcal P=\{x_j\}_{j=1}^{N_s}$29 for EMOS and $\mathcal P=\{x_j\}_{j=1}^{N_s}$30 for GMA, while interval widths at 81% nominal shrank from about $\mathcal P=\{x_j\}_{j=1}^{N_s}$31 under EMOS to $\mathcal P=\{x_j\}_{j=1}^{N_s}$32 under adaptive EMOS, with coverage within 1–2% of nominal (Scheuerer et al., 2013). Spatially informed neural post-processing via convolutional autoencoders further showed that adding low-dimensional bottleneck features from full forecast fields reduced test-set CRPS from $\mathcal P=\{x_j\}_{j=1}^{N_s}$33 for the plain DRN to $\mathcal P=\{x_j\}_{j=1}^{N_s}$34 for DRN+ConvAE at bottleneck size $\mathcal P=\{x_j\}_{j=1}^{N_s}$35, while combined T2M+Z500 input at $\mathcal P=\{x_j\}_{j=1}^{N_s}$36 yielded CRPS approximately $\mathcal P=\{x_j\}_{j=1}^{N_s}$37; positive skill appeared at 96% of stations (Lerch et al., 2022).

Non-atmospheric examples confirm that the same principles are not confined to meteorology. In Spatialize’s non-gridded mining example, expert Global Ordinary Kriging with a nested spherical variogram and an automated PyKrige search both reached $\mathcal P=\{x_j\}_{j=1}^{N_s}$38, while ESI-Kriging produced a sharp, geologically credible map without global variogram modeling (Egaña et al., 23 Jul 2025). In epidemic tracking, EnOSI assimilated data every 10 time steps with $\mathcal P=\{x_j\}_{j=1}^{N_s}$39 ensemble members and successfully detected the emergence of a new Denver focus at $\mathcal P=\{x_j\}_{j=1}^{N_s}$40, even though the forecast ensemble had not yet discovered it; no quantitative error statistics were tabulated, and the evidence was presented through truth-versus-forecast-versus-analysis plots (Krishnamurthy et al., 2010).

6. Limitations, misconceptions, and methodological directions

A common misconception is to equate spatial interpolation exclusively with a single global kriging model. The methods summarized here use local IDW or local kriging inside stochastic forests, bootstrap DNNs on radial-basis spatial embeddings, weighted ensembles of station-specific neural submodels, stacked quantile regressors, and forecast-feature clustering for semi-local predictive distributions. This suggests that ensemble spatial interpolation is better understood as a class of aggregation strategies than as a single covariance specification (Egaña et al., 23 Jul 2025, Nag et al., 2023, Papacharalampous et al., 2024, Baran et al., 2024).

The principal methodological tensions are tuning, representation, and nonstationarity. In clustering-based EMOS, the choice of $\mathcal P=\{x_j\}_{j=1}^{N_s}$41 and feature dimension may require tuning per application; clustering purely on forecasts may misclassify sites with subtle climatological differences; and rolling window length trades off adaptivity versus statistical stability (Baran et al., 2024). In quantile-regression ensembles, quantile crossing must be prevented by enforcing $\mathcal P=\{x_j\}_{j=1}^{N_s}$42 for $\mathcal P=\{x_j\}_{j=1}^{N_s}$43 (Papacharalampous et al., 2024). In Spatialize, larger $\mathcal P=\{x_j\}_{j=1}^{N_s}$44 improves stability at linear cost, and hyperparameters such as $\mathcal P=\{x_j\}_{j=1}^{N_s}$45, $\mathcal P=\{x_j\}_{j=1}^{N_s}$46, nugget, range, and variogram model are selected by automated cross-validation (Egaña et al., 23 Jul 2025). In EnOSI, the time-invariant covariance is cheaper, but if the spatial correlation structure changes drastically over time, the recommendation is to switch to an EnKF or periodically re-estimate $\mathcal P=\{x_j\}_{j=1}^{N_s}$47 and $\mathcal P=\{x_j\}_{j=1}^{N_s}$48 (Krishnamurthy et al., 2010).

A further limitation concerns the information lost by pointwise interpolation of gridded forecast fields to stations. The ConvAE results show that classical interpolation from full fields to local station predictors discards spatial coherence beyond the local vicinity, while bottleneck features can preserve large-scale patterns such as ridges, troughs, and blocking regimes (Lerch et al., 2022). This suggests a broader design principle: ensemble spatial interpolation need not operate only on physical neighborhoods; it may also be constructed in feature space, latent space, or forecast-similarity space. Proposed extensions in the cited literature include dynamic or online clustering, additional covariates such as terrain height and roughness, multivariate post-processing, interpolation for dual-resolution ensembles, and dimensionality reduction through spatial embedding or autoencoders (Baran et al., 2024, Papacharalampous et al., 2024).