Uncertainty quantification methods
We organize UQ methods along two axes by whether uncertainty is intrinsic to the model (produced during training) or applied post hoc, and what type of uncertainty each method captures. Table 9.1 provides a summary and we discuss each family below. The methods below represent different points on the tradeoff between statistical rigor and computational feasibility, and no existing method occupies the ideal corner of that space.
| Method | Captures | Overhead | EO failure modes | Calibrated? |
|---|---|---|---|---|
| Gaussian Processes | Both | \(O(n^3)\) | Kernel misspecification or non-stationarity | ✓\(^*\) |
| Heteroskedastic NN | Aleatoric | \(1\times\) | Blind to epistemic; weak variance signal | × |
| Distributional regression | Aleatoric | \(1\times\) | Distribution family misspecification | × |
| Neural Processes | Both | \(1\times\) | VI approximation error | \(\sim\) |
| Deep Ensembles | Epistemic | \(N\times\) | Correlated members from data | × |
| MC Dropout | Epistemic | \(K\times\) | Systematic underestimation | × |
| Cheap ensembles\(^\dagger\) | Epistemic | \(1.2\)–\(2\times\) | Same as ensembles, reduced | × |
| Quantile Regression | Aleatoric | \(1\)–\(Q\times\) | No calibration guarantee; crossing | × |
| Conformal Prediction | Marginal | \(1\times\) | Only marginal coverage; exchangeability | ✓\(^\ddagger\) |
| CQR | Both\(^\S\) | \(1\times\) | Marginal not conditional coverage | ✓\(^\ddagger\) |
Intrinsic methods
Gaussian Processes
Gaussian processes (GPs) provide the theoretically principled baseline for spatial prediction with uncertainty [1]. Under a well-specified kernel, the GP posterior yields exact, calibrated predictive distributions including at unobserved locations, which is the inference setting for map generation. The classical limitation is computational as exact GP inference scales as \(O(n^3)\) in the number of observations (pixels usually), which is prohibitive for the millions-to-billions of data points typical of global EO datasets.
Sparse approximations based on inducing points [2], [3] reduce the cost to \(O(nm^2)\) where \(m \ll n\), making GPs feasible for moderately large datasets. However, the choice of kernel function is significant and the placement of inducing points is itself also an optimization problem that interacts with kernel specification. A fixed inducing set may adequately represent covariance structure in data-dense regions while underrepresenting data-sparse ones, precisely where uncertainty estimates matter most. For non-stationary EO data, this means the user must jointly navigate kernel selection, inducing point placement, and potential domain partitioning, which often require domain expertise.
EO data are non-stationary (tropical forests and savannas have fundamentally different covariance structures), and misspecifying the kernel undermines the calibration that is the GP’s principal advantage. The geostatistical literature, where GPs appear under the name kriging, has a long history of addressing non-stationarity and change-of-support problems [4], [5] and these insights remain relevant even as the ML community increasingly favors neural approaches.
Heteroskedastic and Distributional Regression
The simplest neural approach to uncertainty is to add a second output head predicting the variance alongside the mean, trained via a Gaussian negative log-likelihood loss [6]. This captures aleatoric uncertainty cheaply (because the overhead is a single additional output neuron) but is blind to epistemic uncertainty. If the model is systematically wrong in some region (e.g., due to training data gaps), the variance head will not reflect this. The training signal for the variance head is also weak since the model must learn the noise structure from the same data it uses to learn the function, and the loss landscape for the variance parameter is poorly behaved when the model is misspecified.
More expressive approaches predict the full parameters of a flexible distribution rather than just a Gaussian mean and variance. For skewed variables such as biomass, where a Gaussian assumption is often inappropriate even after log-transformation, mixture density networks [7] or distributions with explicit skewness and kurtosis parameters (e.g., sinh-arcsinh normal) [8] offer a richer representation of predictive uncertainty.
Variational Inference and Neural Processes
Neural processes (NPs) [9], [10] can be understood as amortized variational approximations to Gaussian processes. They learn to map a set of context observations to a predictive distribution, much as kriging maps observations to a conditional Gaussian, but with the kernel learned implicitly rather than specified parametrically. This connection to the geostatistical framework that underpins decades of spatial prediction is underappreciated, in part because the NP literature has developed largely within the meta-learning community using different vocabulary and training conventions. From the perspective of an EO researcher already familiar with kriging, the core proposition is not a new uncertainty framework but a scalable approximation to an existing one, with \(O(n)\) inference cost and the capacity to learn non-stationary covariance structure directly from data.
The standard caveats of variational inference apply. The approximate posterior can underestimate uncertainty if the variational family is too restrictive, and the ELBO objective can trade off data fit against KL divergence in ways that are not transparent. In practice, design choices specific to the NP framework can produce well-calibrated prediction intervals. As an example, [11], [12] demonstrated that attentive neural processes produce residual standard deviations near the ideal value for biomass mapping across ecologically distinct regions, including under active deforestation. However, these remain among the few validated applications of NPs in large-scale EO mapping, and broader adoption has been limited by the absence of mature tooling and the unfamiliarity of the meta-learning paradigm.
Post-training methods
Ensembles
Deep ensembles [13] train \(N\) models with different random seeds or bootstrap samples and take the spread as an uncertainty estimate. Their dominance in practice is largely explained by their simplicity. Ensembles require no changes to the base model architecture, loss function, or training pipeline. One simply trains the same model multiple times and measures disagreement. Random forests offer an even lower barrier, since every tree is already an ensemble member and quantile regression forests [14] formalize this by retaining the full distribution of training samples at leaf nodes. For those already using tree-based methods for mapping, uncertainty estimates are essentially a byproduct.
As an example, [15] estimated uncertainty in global canopy height predictions using an ensemble of heteroskedastic CNNs, where each member predicts both a mean and a variance via Gaussian negative log-likelihood. Aleatoric uncertainty is estimated by averaging the individual variance predictions, and epistemic uncertainty by the variance of the individual means with a decomposition that follows the law of total variance. The approach illustrates both the appeal and the cost, as the heteroskedastic heads required no post-hoc calibration step, but inference at global scale meant running ten forward passes over the full domain.
However, ensemble spread measures model disagreement, which is related to but distinct from prediction interval width. For bootstrap based ensembles, including random forests, samples drawn from spatially autocorrelated data overlap substantially, and [16] showed that the resulting confidence intervals are only valid under conditions that such data systematically violate. Deep ensembles face a more severe form of the same problem class as members typically share the entire training set and differ only in random initialization, so the sole source of disagreement is optimization stochasticity. When the training data contain spatially structured biases as EO data invariably do, all members learn similar biases in the same regions and agree with each other where they are collectively wrong. In both cases, the mechanism parallels the well documented effect by which spatial autocorrelation inflates apparent model accuracy [17], [18]: the same data structure that makes test-set \(R^2\) overly optimistic makes ensemble spread overly narrow.
Several approaches reduce the computational cost of ensembles without addressing the spread-versus-interval issue. Monte Carlo dropout [19] approximates an ensemble by applying dropout at test time, but systematically underestimates uncertainty because it explores a narrow region of weight space. Snapshot ensembles [20] collect checkpoints along a cyclic learning rate schedule. Stochastic Weight Averaging (SWA/SWAG; [21]) averages weights or samples from a Gaussian approximation to the weight posterior. BatchEnsemble [22] shares most parameters across members, reducing cost to approximately \(1.2\times\) a single model. FiLM-Ensemble [23] modulates network activations via feature-wise linear modulation to produce diverse ensemble members from a single base network, achieving calibration comparable to explicit ensembles at a fraction of the memory cost. These techniques change the cost–calibration tradeoff but do not resolve the correlation-induced miscalibration when applied to spatially structured data.
A complementary direction is multi-model ensembling at the foundation model level. As the ecosystem of geospatial foundation models grows, disagreement between independently pretrained representations of the same location could be a proxy for epistemic uncertainty, analogous to classical ensembles but operating in embedding space rather than prediction space. This is appealing because the diversity between foundation models (different training data, different sensors, different pretraining objectives) is far greater than the diversity between ensemble members that differ only in random seed, which is the limitation that can make classical deep ensembles miscalibrated for spatial data. Recent work has shown that fusing embeddings from multiple geospatial foundation models yields complementary information that improves downstream task performance [24], confirming that these models capture meaningfully different representations. However, whether embedding-space disagreement translates into well-calibrated predictive uncertainty remains an open question, and the computational cost of running multiple foundation models at inference may rival or exceed that of traditional ensembles.
Quantile Regression
Rather than deriving prediction intervals from a point-estimate model, quantile regression directly targets conditional quantiles by optimizing a pinball (check) loss [25]. One can train separate models for the \(\tau\)-th and \((1-\tau)\)-th quantiles, or use a single model with multiple quantile heads, to obtain prediction intervals at any nominal level. Quantile Random Forests [14] retain the full distribution of training samples at leaf nodes to estimate conditional quantiles non-parametrically.
The appeal is directness as the model is optimized for the quantity of interest (interval bounds) rather than deriving intervals from variance estimates. The limitation is that there is no guarantee the resulting intervals are calibrated. A model trained for the 5th and 95th percentiles may achieve empirical coverage far from 90%. Quantile crossing (the upper bound falling below the lower bound) can also occur, though isotonic regression or architectural constraints can mitigate this.
Conformal Prediction
Conformal prediction [26] provides distribution-free coverage guarantees. Given exchangeable calibration data, a conformal procedure produces prediction sets that contain the true value with at least the nominal probability. The accessible tutorial by [27] has driven recent adoption in the EO community [28], [29].
The appeal is theoretical cleanliness as coverage is guaranteed without distributional assumptions, and practical cheapness since conformal methods operate as a post-hoc wrapper around any base predictor. The limitation for spatial mapping is that the guarantee is marginal. The coverage holds on average over the full test distribution, not conditionally for specific regions or subpopulations. For a map user who needs to know how uncertain a given pixel is in the Amazon, a guarantee that averages over the entire globe is of limited value.
More fundamentally, standard conformal prediction assumes exchangeability of calibration and test points, which spatial autocorrelation violates. Spatially correlated data are not exchangeable: nearby points carry redundant information, and the effective calibration sample size is smaller than the nominal one. This can lead to coverage that is correct marginally but arbitrarily poor conditionally [30]. To address this limitation, [31] proposed localized conformal prediction, which emphasizes a local region around each test point when computing prediction sets, retaining the finite-sample marginal coverage guarantee while offering additional local coverage guarantees under suitable assumptions. Adapting such methods to the data volumes and non-stationarities of global EO maps remains an open problem, but they represent a promising direction beyond purely marginal guarantees.
Conformalized quantile regression (CQR; [32]) offers a middle ground. CQR uses conformal calibration to adjust model-produced quantile intervals, inheriting the adaptivity of the base quantile model while gaining the marginal coverage guarantee. [33] applied conformal prediction with a regression model to forest attribute mapping, using split conformal inference to construct prediction intervals and evaluating empirical coverage on withheld test data. While CQR does not solve the conditional coverage problem for spatial data, it is a practical post-hoc calibration approach that is currently available.
The methods above are framed primarily for continuous-valued predictions, but many EO maps are categorical like land cover, crop type, or change detection, and UQ for classification raises distinct issues. The natural uncertainty output is a vector of class probabilities, typically from a softmax layer, but these are usually miscalibrated and neural networks tend to produce overconfident probability vectors where the maximum class probability is systematically higher than the true likelihood of being correct [34]. Temperature scaling is the standard fix and is effective for its simplicity in classification where it is appropriate. Focal loss is a training time alternative that produces better-calibrated models than cross entropy by implicitly maximizing prediction entropy [35]. Platt scaling and isotonic regression provide further options when miscalibration is non-uniform across the probability range [36].