Spatio-temporal DeepKriging (STDK)

• Spatio-temporal DeepKriging is a deep learning framework that extends the spatial basis embedding logic of classical kriging to jointly model space and time.
• It utilizes a two-stage pipeline where a feed-forward network first interpolates historical data and then QLSTM or QConvLSTM modules forecast future quantiles.
• Empirical results show significant improvements in MSPE and predictive coverage over conventional methods, with applications ranging from PM2.5 forecasting to European precipitation analysis.

Spatio-temporal DeepKriging (STDK) is a deep neural approach to spatio-temporal interpolation and probabilistic forecasting that replaces explicit covariance modeling in classical kriging and Gaussian-process formulations with supervised learning on spatio-temporal basis embeddings. In its canonical form, STDK extends the basis-expansion logic of DeepKriging from purely spatial prediction to a space-time index $(\mathbf{s},t)$, and, in the 2023 formulation called “Space-Time.DeepKriging,” organizes prediction into a two-stage pipeline: basis-embedded interpolation or imputation over observed time points, followed by sequence forecasting with quantile LSTM or quantile convolutional LSTM modules (Chen et al., 2020, Nag et al., 2023).

1. Definition and historical positioning

The immediate precursor to STDK is DeepKriging, which introduced the core design pattern of representing spatial dependence through a basis-function embedding of location and then learning a nonlinear predictor with a deep neural network. In that formulation, the spatial random effect is motivated by a truncated Karhunen–Loève expansion, so that spatial structure enters the model as engineered features rather than through a fitted covariance matrix (Chen et al., 2020).

STDK generalizes that blueprint to spatio-temporal prediction. The 2023 “Spatio-temporal DeepKriging for Interpolation and Probabilistic Forecasting” paper defines a real-valued spatio-temporal random field $Y(\mathbf{s},t)$, noisy observations $Z(\mathbf{s},t)=Y(\mathbf{s},t)+\epsilon(\mathbf{s},t)$, and a two-stage procedure in which a basis-embedded feed-forward network first reconstructs or interpolates a complete historical signal at a target location, after which QLSTM or QConvLSTM models forecast future quantiles from that reconstructed history (Nag et al., 2023).

In a narrow architectural sense, STDK refers to this basis-expansion-plus-neural-network family. In a broader literature sense, it sits inside a larger movement from covariance-based kriging toward learned spatio-temporal field inference. That broader family also contains graph-based inductive kriging models. However, papers such as STA-GANN are best viewed as conceptually related but architecturally distinct: they solve the same high-level task of inferring spatio-temporal signals at unobserved sensors, yet do so through dynamic graph construction, message passing, phase alignment, and adversarial transfer rather than through canonical DeepKriging-style basis embeddings (Li et al., 22 Aug 2025).

2. Basis-embedded spatio-temporal formulation

The canonical STDK formulation begins with the decomposition

$Z(\mathbf{s},t)=Y(\mathbf{s},t)+\epsilon(\mathbf{s},t), \qquad Y(\mathbf{s},t)=\mu(\mathbf{s},t)+\gamma(\mathbf{s},t),$

where $\epsilon(\mathbf{s},t)$ is white noise and $\gamma(\mathbf{s},t)$ is the latent spatio-temporal residual process. The residual is motivated through a truncated basis expansion,

$\gamma(\mathbf{s},t)=\sum_{p=1}^{\infty} w_p \boldsymbol{\phi}_p(\mathbf{s},t) \approx \sum_{p=1}^{Q} w_p \boldsymbol{\phi}_p(\mathbf{s},t).$

Although tensor-product bases $\boldsymbol{\phi}_p(\mathbf{s},t)=\phi_i(\mathbf{s})\psi_j(t)$ are discussed conceptually, the implemented STDK embedding stacks spatial and temporal basis evaluations,

$\boldsymbol{\phi}(\mathbf{s},t)= \{\phi_1(\mathbf{s}),\ldots,\phi_G(\mathbf{s}),\psi_1(t),\ldots,\psi_H(t)\}^T,$

and then concatenates them with optional covariates into $\mathbf{X}_\phi(\mathbf{s},t)$ (Nag et al., 2023).

This stacked embedding is central to the method’s interpretation. It reduces dependence learning to a feature-learning problem over a structured basis layer, while avoiding explicit specification of a spatio-temporal covariance function. The paper explicitly states that this stacking does not imply separability, because the downstream DNN shares weights across all basis coordinates and can learn nonlinear interactions between space and time features (Nag et al., 2023).

The spatial basis functions are compactly supported Wendland radial bases,

$Y(\mathbf{s},t)$0

with

$Y(\mathbf{s},t)$1

where $Y(\mathbf{s},t)$2 are spatial anchor points and $Y(\mathbf{s},t)$3 is set to $Y(\mathbf{s},t)$4 times the distance between adjacent anchors. Temporal dependence is encoded with Gaussian radial basis functions,

$Y(\mathbf{s},t)$5

where $Y(\mathbf{s},t)$6 are temporal anchor points and $Y(\mathbf{s},t)$7 equals the spacing between adjacent temporal anchors (Nag et al., 2023).

This basis-embedded view makes STDK closely related to low-rank geostatistical models such as FRK and to the original DeepKriging construction, but the predictor class is no longer linear in the basis evaluations. A plausible implication is that STDK should be read less as a neural approximation to a specific covariance model than as a nonlinear regression architecture whose inductive bias is supplied by multi-resolution spatio-temporal basis functions.

3. Two-stage architecture and probabilistic prediction

Stage 1 of STDK is a feed-forward deep neural network operating on $Y(\mathbf{s},t)$8. For an $Y(\mathbf{s},t)$9-layer architecture, hidden layers use ReLU, and the output is

$Z(\mathbf{s},t)=Y(\mathbf{s},t)+\epsilon(\mathbf{s},t)$0

with parameters estimated by minimizing the quantile check loss

$Z(\mathbf{s},t)=Y(\mathbf{s},t)+\epsilon(\mathbf{s},t)$1

For deterministic interpolation, the median quantile $Z(\mathbf{s},t)=Y(\mathbf{s},t)+\epsilon(\mathbf{s},t)$2 is used; for probabilistic interpolation, additional lower and upper quantiles are estimated (Nag et al., 2023).

Stage 2 treats the Stage-1 interpolated history at a target location as a completed time series,

$Z(\mathbf{s},t)=Y(\mathbf{s},t)+\epsilon(\mathbf{s},t)$3

and feeds that sequence into either QLSTM or QConvLSTM. The QLSTM uses the standard LSTM gating equations, while QConvLSTM augments forecasting with local $Z(\mathbf{s},t)=Y(\mathbf{s},t)+\epsilon(\mathbf{s},t)$4 spatial neighborhoods of Stage-1 predictions, followed by $Z(\mathbf{s},t)=Y(\mathbf{s},t)+\epsilon(\mathbf{s},t)$5 convolutions with 64 feature maps before the recurrent stack (Nag et al., 2023).

The uncertainty mechanism is quantile-based rather than likelihood-based. STDK does not produce a full predictive density in closed form; instead it estimates conditional quantiles and constructs prediction intervals from them. To avoid quantile crossing, the 2023 paper uses a specialized output transformation centered on the median:

$Z(\mathbf{s},t)=Y(\mathbf{s},t)+\epsilon(\mathbf{s},t)$6

This ensures ordered lower and upper quantiles around the central forecast (Nag et al., 2023).

The two-stage organization is an important defining feature. Interpolation and forecasting are not jointly optimized as a single end-to-end dynamical model; rather, Stage 2 is conditional on Stage 1. The 2023 paper explicitly notes that forecast uncertainty therefore does not propagate interpolation uncertainty through a full joint probabilistic model (Nag et al., 2023).

4. Implementations, applications, and empirical demonstrations

The 2023 STDK paper reports results on both simulated and real datasets. On the KAUST competition datasets, STDK achieved RMSPE $Z(\mathbf{s},t)=Y(\mathbf{s},t)+\epsilon(\mathbf{s},t)$7 on sub-competition 2a and $Z(\mathbf{s},t)=Y(\mathbf{s},t)+\epsilon(\mathbf{s},t)$8 on sub-competition 2b, corresponding to a $Z(\mathbf{s},t)=Y(\mathbf{s},t)+\epsilon(\mathbf{s},t)$9 improvement and a $Z(\mathbf{s},t)=Y(\mathbf{s},t)+\epsilon(\mathbf{s},t), \qquad Y(\mathbf{s},t)=\mu(\mathbf{s},t)+\gamma(\mathbf{s},t),$0 improvement over the second-best team, respectively. In the nonstationary simulation, “Space-Time.DeepKriging” obtained interpolation MSPE $Z(\mathbf{s},t)=Y(\mathbf{s},t)+\epsilon(\mathbf{s},t), \qquad Y(\mathbf{s},t)=\mu(\mathbf{s},t)+\gamma(\mathbf{s},t),$1 versus $Z(\mathbf{s},t)=Y(\mathbf{s},t)+\epsilon(\mathbf{s},t), \qquad Y(\mathbf{s},t)=\mu(\mathbf{s},t)+\gamma(\mathbf{s},t),$2 for GpGp, while QConvLSTM achieved forecasting MSPE $Z(\mathbf{s},t)=Y(\mathbf{s},t)+\epsilon(\mathbf{s},t), \qquad Y(\mathbf{s},t)=\mu(\mathbf{s},t)+\gamma(\mathbf{s},t),$3, MPIW $Z(\mathbf{s},t)=Y(\mathbf{s},t)+\epsilon(\mathbf{s},t), \qquad Y(\mathbf{s},t)=\mu(\mathbf{s},t)+\gamma(\mathbf{s},t),$4, and coverage $Z(\mathbf{s},t)=Y(\mathbf{s},t)+\epsilon(\mathbf{s},t), \qquad Y(\mathbf{s},t)=\mu(\mathbf{s},t)+\gamma(\mathbf{s},t),$5 (Nag et al., 2023).

The large PM$Z(\mathbf{s},t)=Y(\mathbf{s},t)+\epsilon(\mathbf{s},t), \qquad Y(\mathbf{s},t)=\mu(\mathbf{s},t)+\gamma(\mathbf{s},t),$6 application in that paper used monthly U.S. EPA data from January 1999 to December 2022 at more than $Z(\mathbf{s},t)=Y(\mathbf{s},t)+\epsilon(\mathbf{s},t), \qquad Y(\mathbf{s},t)=\mu(\mathbf{s},t)+\gamma(\mathbf{s},t),$7 space-time locations. Reported interpolation performance was MSPE $Z(\mathbf{s},t)=Y(\mathbf{s},t)+\epsilon(\mathbf{s},t), \qquad Y(\mathbf{s},t)=\mu(\mathbf{s},t)+\gamma(\mathbf{s},t),$8 for STDK versus $Z(\mathbf{s},t)=Y(\mathbf{s},t)+\epsilon(\mathbf{s},t), \qquad Y(\mathbf{s},t)=\mu(\mathbf{s},t)+\gamma(\mathbf{s},t),$9 for GpGp. For forecasting, QConvLSTM achieved MSPE $\epsilon(\mathbf{s},t)$0, MPIW $\epsilon(\mathbf{s},t)$1, and coverage $\epsilon(\mathbf{s},t)$2, while QLSTM had MSPE $\epsilon(\mathbf{s},t)$3 and MPIW $\epsilon(\mathbf{s},t)$4. The authors also report runtimes of about $\epsilon(\mathbf{s},t)$5 minutes on a Tesla P100 GPU for the competition datasets, $\epsilon(\mathbf{s},t)$6 minutes on GPU for PM$\epsilon(\mathbf{s},t)$7 interpolation, and $\epsilon(\mathbf{s},t)$8 minutes on a 16-core CPU (Nag et al., 2023).

A later supplementary paper implemented STDK in PyTorch for daily European precipitation and emphasized interpolation plus probabilistic forecasting on irregular station data. In that implementation, interpolation used four spatial Wendland resolutions with nominal counts $\epsilon(\mathbf{s},t)$9, $\gamma(\mathbf{s},t)$0, $\gamma(\mathbf{s},t)$1, and $\gamma(\mathbf{s},t)$2, pruned to $\gamma(\mathbf{s},t)$3 spatial basis functions after removing water-body bases, together with $\gamma(\mathbf{s},t)$4 temporal Gaussian basis functions. The resulting input dimension was approximately $\gamma(\mathbf{s},t)$5, and the network output three values per sample: predictive median and the lower and upper bounds of a $\gamma(\mathbf{s},t)$6 prediction interval (Nag, 16 Sep 2025).

That supplementary study reported interpolation over approximately $\gamma(\mathbf{s},t)$7 million space-time locations at 10-day intervals for Europe, after a 10-day moving average and standardization of precipitation observations from 2014-01-01 to 2024-12-31. Reported interpolation metrics were MSPE $\gamma(\mathbf{s},t)$8, PICP $\gamma(\mathbf{s},t)$9, and MPIW $\gamma(\mathbf{s},t)=\sum_{p=1}^{\infty} w_p \boldsymbol{\phi}_p(\mathbf{s},t) \approx \sum_{p=1}^{Q} w_p \boldsymbol{\phi}_p(\mathbf{s},t).$0. For the forecasting comparison on a $\gamma(\mathbf{s},t)=\sum_{p=1}^{\infty} w_p \boldsymbol{\phi}_p(\mathbf{s},t) \approx \sum_{p=1}^{Q} w_p \boldsymbol{\phi}_p(\mathbf{s},t).$1 subregion, ConvLSTM obtained MSPE $\gamma(\mathbf{s},t)=\sum_{p=1}^{\infty} w_p \boldsymbol{\phi}_p(\mathbf{s},t) \approx \sum_{p=1}^{Q} w_p \boldsymbol{\phi}_p(\mathbf{s},t).$2, PICP $\gamma(\mathbf{s},t)=\sum_{p=1}^{\infty} w_p \boldsymbol{\phi}_p(\mathbf{s},t) \approx \sum_{p=1}^{Q} w_p \boldsymbol{\phi}_p(\mathbf{s},t).$3, and MPIW $\gamma(\mathbf{s},t)=\sum_{p=1}^{\infty} w_p \boldsymbol{\phi}_p(\mathbf{s},t) \approx \sum_{p=1}^{Q} w_p \boldsymbol{\phi}_p(\mathbf{s},t).$4, whereas STDK obtained MSPE $\gamma(\mathbf{s},t)=\sum_{p=1}^{\infty} w_p \boldsymbol{\phi}_p(\mathbf{s},t) \approx \sum_{p=1}^{Q} w_p \boldsymbol{\phi}_p(\mathbf{s},t).$5, PICP $\gamma(\mathbf{s},t)=\sum_{p=1}^{\infty} w_p \boldsymbol{\phi}_p(\mathbf{s},t) \approx \sum_{p=1}^{Q} w_p \boldsymbol{\phi}_p(\mathbf{s},t).$6, and MPIW $\gamma(\mathbf{s},t)=\sum_{p=1}^{\infty} w_p \boldsymbol{\phi}_p(\mathbf{s},t) \approx \sum_{p=1}^{Q} w_p \boldsymbol{\phi}_p(\mathbf{s},t).$7 (Nag, 16 Sep 2025).

The supplementary paper is also notable for its implementation detail and its ambiguities. It states that the interpolation model has “10 hidden layers,” while explicitly listing five hidden layers of width 100 and four hidden layers of width 50, which sums to nine. It also under-specifies the forecasting STDK architecture and defers several equations, including the exact quantile activation, to the main STDK paper. This suggests that the PyTorch supplement is best read as an implementation-oriented application note rather than a complete standalone methodological specification (Nag, 16 Sep 2025).

5. Relation to adjacent neural kriging frameworks

Within the larger neural kriging literature, STDK occupies the basis-expansion branch. Graph-based methods such as SATCN, INCREASE, and STA-GANN address the same family of tasks—inferring spatio-temporal signals at unobserved locations—but they do so through masked graph aggregation, relation-aware recurrent units, dynamic graph learning, or adversarial transfer rather than through fixed spatio-temporal basis embeddings (Wu et al., 2021, Zheng et al., 2023, Li et al., 22 Aug 2025).

SATCN represents kriging as masked graph message passing plus temporal convolution. Its first SAN layer prevents unsampled nodes from sending messages, after which TCN blocks model temporal dynamics over the imputed node features. INCREASE constructs target-specific inductive representations from multiple relations—spatial proximity, functional similarity, and transition probability—then uses relation-aware GRUs and multi-relation attention. STA-GANN extends the graph-based line further with Dynamic Data-Driven Metadata Graph Modeling, the Decoupled Phase Module for timestamp shift, and adversarial transfer learning for unknown-sensor generalization (Wu et al., 2021, Zheng et al., 2023, Li et al., 22 Aug 2025).

These models are not canonical STDK in the narrow architectural sense. The relation is conceptual rather than nominal: they share the objective of replacing explicit covariance-based kriging with learned spatio-temporal inference, but their dependence modeling is graph-induced rather than basis-embedded. The STA-GANN analysis in particular is explicit that the method is “conceptually yes” but “architecturally no” with respect to STDK: it belongs to the same broader movement away from explicit covariance specification, yet it is better characterized as inductive graph neural kriging than as standard DeepKriging (Li et al., 22 Aug 2025).

At the same time, STDK itself has begun to develop internal extensions. A 2026 paper proposes a density-adaptive multi-quantile STDK variant with trainable cluster-adaptive spatial bases and cluster-aware conformal calibration. In that formulation, the standard STDK input $\gamma(\mathbf{s},t)=\sum_{p=1}^{\infty} w_p \boldsymbol{\phi}_p(\mathbf{s},t) \approx \sum_{p=1}^{Q} w_p \boldsymbol{\phi}_p(\mathbf{s},t).$8 is retained, but fixed regular-grid spatial bases are replaced with basis functions whose centers and scales are initialized from spatial sampling density, and prediction intervals are calibrated within spatial clusters rather than by a single global conformal correction (Kim et al., 4 Jun 2026).

This extension is important because it directly targets a weakness of standard basis-grid STDK under clustered sampling. The paper reports that under clustered observation patterns, standard STDK plus global conformal calibration can severely undercover, whereas the adaptive-basis variant improves both coverage and worst-site reliability. Under uniform sampling, by contrast, gains are small or absent. This suggests that the principal fault line inside the STDK family is no longer only covariance versus deep learning, but also fixed geometric bases versus sampling-adaptive bases (Kim et al., 4 Jun 2026).

6. Statistical precursors, limitations, and open directions

Two classical frameworks are especially relevant for situating STDK. FRK provides a low-rank spatial random effects model with automatically generated or user-specified spatial and spatio-temporal basis functions, BAU discretization, EM estimation, and exact Gaussian uncertainty propagation over arbitrary supports. The SPDE advection-diffusion approach provides a sparse-precision, physics-informed route to nonseparable and asymmetric spatio-temporal dependence, together with kriging means, conditional variances, and conditional simulations (Zammit-Mangion et al., 2017, Clarotto et al., 2022).

Compared with FRK, STDK inherits the basis-expansion intuition but discards the Gaussian latent-variable structure, support-change machinery, and analytic posterior uncertainty. Compared with the SPDE approach, STDK replaces explicit transport-diffusion dynamics and sparse precision matrices with a more flexible but less interpretable neural approximation. These contrasts are not merely philosophical: FRK and SPDE retain principled distinction between measurement error and latent-process uncertainty, while STDK’s uncertainty is quantile-based and operational rather than model-based (Zammit-Mangion et al., 2017, Clarotto et al., 2022).

The main limitations explicitly identified in the STDK literature are therefore structural. The 2023 paper states that Stage-2 forecasting is conditional on Stage-1 interpolation and does not propagate interpolation uncertainty through a joint model; forecast maps can be less spatially smooth because forecasting is performed location by location; basis choices are fixed rather than learned; and QConvLSTM is more computationally expensive than QLSTM (Nag et al., 2023). The PyTorch supplementary application adds practical limitations of under-specified implementation details, including an inconsistent hidden-layer count and a forecasting architecture that is summarized more as a pipeline than as a fully specified model (Nag, 16 Sep 2025).

A more recent limitation is geometric rather than probabilistic. Under highly non-uniform spatial sampling, fixed regular-grid bases can over-allocate representational capacity to sparse regions while under-resolving dense clusters. The density-adaptive conformal extension addresses exactly this issue by learning cluster-adaptive spatial bases and calibrating intervals locally, which indicates that standard STDK’s fixed spatial basis layer is not always efficient under clustered observation designs (Kim et al., 4 Jun 2026).

A plausible implication is that subsequent STDK developments will continue to hybridize three strands that already coexist in the literature: basis-embedded DeepKriging, graph-based inductive kriging, and classical low-rank or physics-informed spatio-temporal structure. The most technically natural directions suggested by the existing papers are adaptive basis construction, localized uncertainty calibration, explicit treatment of unseen-location generalization, and mechanisms that recover some of the uncertainty semantics and support-awareness that remain stronger in FRK and SPDE-based kriging than in current deterministic STDK formulations.