Filling the Missings: Spatiotemporal Data Imputation by Conditional Diffusion

1 Introduction

Spatiotemporal data integrates temporal dynamics with spatial relationships to represent complex real-world phenomena. This data type has proven valuable across diverse domains, from environmental science to urban systems [Yi et al., 2016; Li et al., 2017; Wu et al., 2019]. Despite its importance, the collection of spatiotemporal data faces significant challenges in real-world deployments. Sensor malfunction, system instability, and environmental interference often lead to missing or corrupted measurements. These data quality issues compromise the reliability of information extracted from spatiotemporal systems and impair the performance of downstream applications [Lian et al., 2023; Tan et al., 2018; He et al., 2025a].

* Corresponding author: Junhua Gu, Yude Bai.

Various methods have been proposed to handle missing data in spatiotemporal datasets. Traditional approaches rely on statistical and data analytic techniques, and make assumptions about the underlying probability distributions and use techniques such as linear interpolation [Challu et al., 2023], expectation-maximization (EM) algorithm [Cong and Li, 2024; Bao et al., 2025], or k-nearest neighbors (KNN) [Xu et al., 2024] to estimate missing data. Their effectiveness is constrained by their inherent mathematical foundations, which rely on simplified models of relationships between variables. When applied to real-world data, where relationships between spatial and temporal components often follow complex, nonlinear patterns, these methods fail to reconstruct the true data distribution accurately.

Deep learning methods have also been applied for spatiotemporal data imputation to better process the inherent complexity in spatiotemporal relationships. Recurrent neural networks (RNN) [Shi et al., 2024] model temporal relationships through hidden state propagation, while graph neural networks (GNN) [Agafonov, 2020; Jin et al., 2023] learn spatial dependencies through node-level representations. These methods, however, face fundamental limitations. RNN [Chen et al., 2022a] operate through an autoregressive mechanism, which causes early prediction errors to propagate through subsequent time steps. GNN [Cini et al., 2021] depend on features from neighboring nodes, allowing errors to disseminate across the graph structure. The recursive nature of their processing leads to error accumulation over time. Their reliance on historical data also constrains their ability to capture the full range of possible values during imputation.

Alternatively, data imputation can be treated as a generative task [Lin et al., 2024]. Recent work has shown that diffusion probabilistic models (DDPM) [Kingma, 2013] and generative adversarial networks (GAN) [Goodfellow et al., 2020]. However, current DDPM-based methods face two key challenges. First, they fail to model both temporal and spatial aspects simultaneously [Wu et al., 2019] - for example, Tashiro et al. [Tashiro et al., 2021] consider only temporal features, ignoring crucial spatial relationships. Second, they rely on problematic pre-imputation approaches - Liu et al. [Liu et al., 2023] use linear interpolation that creates artificial fluctuations, causing generated data to deviate from actual patterns. Our work addresses these limitations by developing a model that integrates both spatial and temporal information while avoiding error-prone pre-imputation steps.

In this paper, we introduce CoFILL (Conditional Diffusion Model based on Temporal-Frequency Spatiotemporal Imputation), a novel approach that addresses the limitations of current methods. Our model enhances imputation by utilizing temporal measurements from sensors and geographical connections between monitoring locations. The temporal data captures changing patterns over time, while spatial relationships reveal location interdependencies. The architecture features modules for temporal and frequency domain feature extraction to capture rapid changes and persistent patterns. A cross-attention mechanism combines these features into rich conditional information. This information guides the noise prediction network to transform random noise into meaningful data values through iterative denoising, which enables robust imputation by learning the underlying data distribution rather than direct estimation.

Our contributions are summarized as follows:

• We propose a spatiotemporal imputation framework, CoFILL, based on diffusion models. It applies the nonrecursive property of diffusion models to effectively reduce error accumulation during the imputation process.
• We develop a dual-stream feature processing architecture to model complete spatiotemporal dynamics. We extract temporal and frequency information through separate pathways, then integrate them by cross-attention to capture both rapid fluctuations and underlying patterns.
• We conduct extensive experiments on three real-world datasets to validate our approach. The results demonstrate that CoFILL outperforms current state-of-the-art methods in spatiotemporal data imputation tasks.

2 RelatedWorks

2.1 Spatiotemporal Data Imputation

Methods for spatiotemporal data imputation mainly fall into three main categories: statistical approaches, machine learning techniques, and deep learning models [He et al., 2024]. Statistical methods, such as mean and median imputation [Amiri and Jensen, 2016], assume the data follows specific probability distributions, which is too restrictive for real-world data. Machine learning approaches provide more flexibility like random forest (RF) [Li et al., 2014] and matrix factorization [Wang et al., 2018]. Other machine learning methods such as Bayesian Gaussian process tensor completion (BGCP) [Chen et al., 2019b] and time-regularized tensor factorization (TRTF) [Baggag et al., 2019] show promise in processing complex data. But they still struggle with nonlinear relationships and scalability to large datasets. Deep learning methods, such as recurrent neural networks (RNN) [Vaswani et al., 2017] and their variants like GRUD [Shahidi et al., 2024], excel at learning temporal patterns but accumulate errors over prediction sequences. Convolutional neural networks, such as TimesNet [Wu et al., 2022], only analyze time series in Euclidean space using Fourier transforms. Graph neural networks [Marisca et al., 2022; Chen et al., 2022b] improve imputation accuracy by modeling spatial relationships, but their reliance on fixed graph structures reduces their adaptability to new scenarios. These existing approaches cannot fully capture the nonlinear dependencies and complex patterns inherent in large-scale spatiotemporal data.

2.2 DiffusionModels for Imputation

Compared with traditional methods and deep learning methods, the generative data imputation method by diffusion models generates high-quality data. The conditional score-based diffusion model (CSDI) [Tashiro et al., 2021] uses a subset of observed data as conditions to generate residual values. Meanwhile, PriSTI [Liu et al., 2023] introduces spatiotemporal relationships and attention mechanisms to achieve better imputation results. GP-VAE [Fortuin et al., 2020] models temporal patterns using Gaussian processes in latent space, while E2GAN [Li et al., 2025] offers a simplified approach to generate plausible values for time series data. However, current diffusion models struggle to effectively extract and use conditional information - the known patterns and relationships in the existing data that guide the generation of missing data. This limitation impacts imputation accuracy because of the generative values that fail to maintain the underlying temporal and spatial dependencies of the data.

2.3 Temporal and Frequency Domain Feature

Feature extraction plays a vital role in spatiotemporal data imputation. Temporal feature extraction through recurrent structures like GRU-D [Shahidi et al., 2024] captures immediate state changes but cannot handle variations at different time scales. Meanwhile, frequency domain approaches via Fourier transforms, such as TimesNet [Wu et al., 2022], identify periodic and cyclical patterns but miss important transient changes. GNN-based methods [Marisca et al., 2022] focus on modeling spatial relationships but overlook these temporal dynamics. We hope to address these issues by integrating temporal and frequency domain features through a crossattention mechanism. The temporal features capture immediate state transitions, while frequency features reveal longterm trends.

3 ProblemFormulation

Consider a spatiotemporal system with N observation nodes (such as air monitoring stations or traffic sensors) measured over L time steps. The measurements form a sequence X1:L = X1,X2, . . . ,XL ∈ RN×L, where Xl ∈ RN represents the values observed at time step l. The spatial relationships between nodes are represented by an adjacency matrix A ∈ RN×N , where Ai,j = 1 indicates a connection between nodes i and j, and Ai,j = 0 indicates no connection. In real-world deployments, some observations may be missing due to sensor failures or system errors. We use a binary mask M ∈ 0, 1N×L to indicate the presence or absence of measurements, where mi,l = 1 denotes a valid observation at node i and time l, while mi,l = 0 indicates missing data. The spatiotemporal imputation problem aims to reconstruct missing data in X given the observed measurements and the adjacency matrix A. The goal is to develop an imputation method that preserves both the temporal evolution of measurements and their spatial relationships across nodes.

4 Methodology

4.1 Preprocessing of CoFILL

CoFILL performs a two-step preprocessing on the input series to ensure high-quality inputs for training. We first mask some of the observed values in the original data based on a random masking strategy. This not only simulates missing patterns but also enhances the model's generalization to diverse scenarios [Tashiro et al., 2021]. Two sets of initial imputed data are then generated from the masked data. The first set uses forward interpolation, as shown in Figure 1. It fills in missing data by using the known values from the previous time step, thus preserving the continuity and periodicity of the spatiotemporal data. This also reduces information loss due to over-smoothing. The second set fills the missing data with Gaussian noise. We introduce noise consistent with the original data distribution to enhance data diversity and introduce appropriate randomness for the model [Fortuin et al., 2020; He et al., 2025b]. These two imputed data sets are adopted as inputs for subsequent modules, which further supports CoFILL in learning spatiotemporal dependencies and missing data distributions.

4.2 Conditional InformationModule

We extract features from the imputed data sets of preprocessing to provide effective conditional information for the noise prediction network. These features, by identifying the spatiotemporal dependencies from both the observed and preimputed data, guide the network to predict noise more accurately. The process works as follows. Initially, we map the imputed data X1 to a latent space Hin = Conv(X1) ∈ Rd×N×L, where Conv(·) refers to a 1 × 1 convolution operation, and d denotes the channel size. Next, we extract features from both the temporal and frequency domains, and the combined feature is obtained through cross-attention.

In the temporal domain, data reflects instantaneous changes, which is crucial for short-term prediction. The temporal feature extraction module separately obtains both temporal and spatial dependencies. We map X1 to the latent space Hin and pass it through a Temporal Convolutional Network (TCN) to learn temporal dependencies. The TCN employs padded 1-dimensional gated causal convolutions, ensuring that the output and input maintain the same shape. The convolution kernel Γτ ∈ RK×Cin×Cout generates outputs P,Q ∈ RCout×N×L of the same shape. The definition and specific formulas of the TCN structure can be, respectively, expressed as follows:

where ⊙ represents the Hadamard product, and σ is the sigmoid activation function, introducing nonlinearity and acting as a gating mechanism. We denote the output of the TCN as ¯Hin ∈ Rd×N×L, the convolution kernel size is k, the dilation size is 1, which indicates short-term feature extraction, and the convolution is expanded along the time axis. We use a Graph Convolutional Network (GCN) to learn spatial dependencies. The definition and specific formulas of the GCN structure can be, respectively, expressed as follows:

where Wi ∈ RCg in×Cg out is a learnable parameter, σ is the sigmoid activation function, and φ(·) represents the aggregation function that determines how neighboring node features are aggregated into the target node. The graph-normalized adjacency matrix is Agcn = D-1/2 (A + I)D-1/2 ∈ RN×N , where I is the identity matrix, and D is the diagonal matrix with Dii = ∑j ∂(A + I)ij . Agcn is the static adjacency matrix, o is the convolution order of the graph, and Conv2D indicates that convolution is used for channel fusion. We denote the output of both TCN and GCN as ˜Hin ∈ Rd×N×L. By processing dynamic characteristics in the temporal dimension with the TCN, we extract the temporal changes in the data. At the same time, we adopt the GCN to learn the relationships between nodes in the spatial dimension. This combination effectively identifies the complex changes in spatiotemporal data over time.

In the frequency domain, data reveals periodicity and trends, which are essential for long-term forecasting and recognizing periodic patterns. Many time series models fail to account for the frequency impact, which leads to an inability to extract the intrinsic features of the series. To this end, we apply the Discrete Cosine Transform (DCT) to extract frequency domain features ˆHin. Unlike Fourier transforms, DCT avoids the Gibbs phenomenon. It also provides superior energy compaction compared to Wavelet transforms [Jiang et al., 2023]. Its specific formula is as follows:

where Hin ∈ Rd×N×L, ˆHin ∈ Rd×N×L, T represents the time step. As m→0, the corresponding DCT component captures low-frequency signals, reflecting long-term stable trends. Conversely, as m→T, it captures high-frequency components, indicating periodic structures. Together, the different frequency components in DCT represent both the longterm trends and periodic patterns within the sequence.

As shown in Figure 2 B), we apply cross-attention to combine the temporal-domain feature ˜Hin and the frequencydomain feature ˆHin into a mixed feature representation. Temporal and frequency-domain features offer different perspectives. Temporal features reveal short-term behavior, while frequency-domain features identify long-term trends. The fusion of both perspectives enables the system to process multiscale data within a unified framework. The temporal feature Q retains time sequence, local changes, and spatial features. It encodes immediate information along the time axis. The frequency-domain feature KV focuses on global periodicity and trends, which reduces noise interference while preserving key signal characteristics. The combination of both captures both local temporal information and global periodic trends. The time-frequency domain attention weights and the final conditioning information Ccon are defined as:

4.3 Noisy EstimationModule

The noise prediction network ε in this module is employed during the reverse denoising process to progressively sample noise, as illustrated in Figure 2 C). It uses conditional information Ccon to guide the extraction of spatiotemporal dependencies in the noise data Xin. Additionally, it leverages an attention mechanism to mitigate the randomness introduced by Gaussian noise injected into the data distribution. The model accepts three inputs, including noise information Xin = (X1 ∥ ˜Xt ), conditional information Ccon, and the adjacencymatrix A. It first generates temporal features Xtem through a temporal attention module ΥT (·). The data is then passed into a spatial attention module ΥS(·), equipped with a GNN, to produce spatial features Xspa. This process is formalized as follows:

Attntem(·) and Attnspa(·) represent temporal and spatial attention, respectively. The mechanisms are the dot-product multi-head self-attention from the transformer architecture. We also put conditional information into the attention calculations to reduce the noise accumulated during the diffusion steps. Specifically, for Attntem(·), the dot-product attention Attntem(QT , KT , VT ) = softmax((QTKT T )/√d) · VT is rewritten as Attntem(AT , VT ) = AT · VT , where AT = softmax((QT KT T )/√d) denotes the attention weights. The attention weight AT is calculated by the conditional information Ccon. We set the inputs QT , KT , and VT as follows:

where W Q T ,W K T ,W V T ∈ Rd×d are learnable projection parameters. Similarly, the spatial attention Attnspa(AS, VS) computes the attention weight in the same manner:

Algorithm 1 Training of CoFILL

Input: Incomplete observed dataX, the adjacencymatrixA, the number of iterations N, the number of diffusion steps T, noise levels sequence ¯αt.

Output: Optimized noise prediction model ϵθ

end for

The noise prediction network consists of multiple stacked attention layers that progressively refine spatial representations. Each layer processes features through a gated activation unit, producing spatial features Xspa that split into two paths. The primary path forms residual connections feeding into subsequent layers, while a secondary path creates skip connections that aggregate features across network depths. The network combines these accumulated features and processes them through two 1-dimensional convolution layers, which contains only the values of the interpolation targets.

4.4 Training and Prediction

We apply a random masking strategy, during the training, to the input observations X to derive imputation targets ˜X, while the remaining observations are used to extract preimputation features X1 as conditional information. Similar to CSDI, we implement masking strategies including point strategy, block strategy, and hybrid strategy. As shown in Algorithm 1, the noise prediction model of CoFILL is trained by iterating over the observed data. Specifically, once obtaining the training imputation targets ˜X and conditional information X1, the training objective for spatiotemporal imputation is:

In the imputation, we deploy the trained noise prediction model ϵθ and the available observation masks M to perform imputation. The imputation targets consist of all missing data in the spatiotemporal data. As shown in Algorithm 2, the trained model will fill in the missing data in the spatiotemporal data by performing reverse diffusion.

5 Experiments

5.1 Datasets

We use three spatiotemporal datasets for imputation algorithms. The AQI-36 dataset [Zheng et al., 2014] serves as a good benchmark due to its inherent missing data patterns and complex spatiotemporal dependencies. The hourly PM2.5

Algorithm 2 Imputation of CoFILL

Input: A sample of incomplete observed data X, the adjacency matrix A, the number of diffusion steps T , the optimized noise prediction model ϵθ.

Output: Missing data of the imputation target Xout

measurements across Beijing's 36 monitoring stations over a full year handle both short-term fluctuations and long-term season trends. The traffic datasets, METR-LA and PEMSBAY [Li et al., 2017], present different scales to test our model's capabilities. METR-LA's 207 sensors and PEMSBAY's 325 sensors represent large-scale networks with dense spatial correlations. Their five-minute sampling frequency verifies model ability to face rapid temporal variations, while their different geographic contexts - Los Angeles and San Francisco Bay Area - validate results across urban traffic.

5.2 Baselines for Comparison

Traditional Statistical Baselines. The fundamental statistical approaches include Mean, which imputes values using historical node averages, and Daily Average (DA), which uses time-specific averages. We also include K-Nearest Neighbors (KNN) [Bai et al., 2023; Xu et al., 2024], which leverages geographical proximity for imputation, and Linear interpolation for time series reconstruction. For robust sequential data processing, we incorporate Kalman Filter (KF), which estimates missing data through state-space modeling.

Matrix and Tensor-Based Baselines. These intermediate approaches bridge classical statistics and modern machine learning by decomposing complex data structures. MICE [White et al., 2011] uses conditional specifications with chained equations for multivariate imputation. Vector Autoregression (VAR) predicts single-step values through temporal dependencies. Time-Regularized Matrix Factorization (TRMF) [Yu et al., 2016] specializes in temporal data patterns, and Bayesian Augmented Tensor Factorization (BATF) [Chen et al., 2019a] integrates spatiotemporal knowledge through tensor decomposition.

Deep Learning Baselines. For probabilistic modeling with uncertainty quantification, we employ V-RIN [Mulyadi et al., 2021] and GP-VAE [Fortuin et al., 2020]. To leverage adversarial learning, we use rGAIN [Yoon et al., 2018], which enhances GAIN through bidirectional processing. We further employ BRITS [Cao et al., 2018] for RNN-based temporal modeling and GRIN for graph neural network processing with bidirectional GRU. Our evaluation also encompasses recent diffusion-based methods: CSDI [Tashiro et al., 2021], which captures feature interactions through transformers, and PriSTI [Liu et al., 2023], which implements spatiotemporal attention mechanisms.

5.3 EvaluationMetrics

We investigate the performance by three complementary metrics - Mean Absolute Error (MAE), Mean Squared Error (MSE), and Continuous Ranked Probability Score (CRPS) [Berrisch and Ziel, 2023]. MAE and MSE provide measures of imputation accuracy by calculating the distance between predicted values and ground truth. CRPS evaluates the quality of probabilistic predictions by measuring how well the predicted probability distributions align with given values.

5.4 Experimental Design

To ensure reproducibility and meaningful evaluation, we establish a methodical data partitioning strategy that respects the temporal characteristics of our datasets. The AQI-36 dataset requires careful temporal consideration due to its seasonal nature. We distribute the data across seasons by selecting March, June, September, and December for testing, which captures seasonal variations throughout the year. The validation set draws from the final 10% of data in February, May, August, and November to maintain seasonal representation. The remaining months form the training set.

To evaluate model performance under different missing data scenarios, we design two types of data corruption schemes. For AQI-36, we implement a simulated failure (SF) pattern that replicates real-world sensor malfunction distributions. For PEMS-BAY and METR-LA traffic datasets, we create controlled missing data scenarios through mask matrices. These scenarios include random point missing (Point), where we mask 25% of observations uniformly at random, and structured block missing (Block), which combines 5% random masking with continuous missing segments. These segments span 1 to 4 hours per sensor and occur with 0.15% probability, simulating extended sensor outages.

Our implementation runs on an NVIDIA RTX 4090 GPU with 24GB VRAM. We optimize the model using Adam with cosine annealing learning rate decay. The learning rate starts at 10-3 and decays to 10-5 . Table 1 details the complete hyperparameter configuration for each dataset. To establish statistical significance, we conduct five independent runs with different random seeds for each experiment.

5.5 Comparative Performance Study

Our experimental evaluation demonstrates CoFILL's consistent superiority across multiple datasets and missing data scenarios. The quantitative results in Table 2 reveal CoFILL's strong performance in terms of MAE and MSE metrics. For probabilistic imputation assessment, Table 3 presents the CRPS comparisons. CoFILL achieves the best performance in 12 out of 15 experimental configurations spanning different metrics and missing patterns. The improvement becomes particularly apparent when compared with PriSTI, the stateof-the-art method. On the METR-LA dataset with Block, CoFILL reduces both MAE and MSE by 10.22%. On the PEMS-BAY point-missing pattern, CoFILL performs on par with PriSTI (only ≤1% difference with identical MSE).

The comparative analysis reveals clear performance patterns across different methodological categories. Traditional statistical and machine learning approaches, such as Mean and KNN, show limited effectiveness due to their simplified assumptions about spatiotemporal relationships. The deep learning approaches, represented by GRIN and BRITS, face performance degradation from accumulated errors in their autoregressive processing. While recent diffusion-based models like CSDI and PriSTI demonstrate strong baseline performance, CoFILL still outperforms them with advanced conditional information construction and feature fusion methods.

5.6 Ablation Study

We evaluate the impact of key components of CoFILL by comparing the following settings.

• w/o Forward (Forward Interpolation): Pre-imputation information X1 is removed, and conditional feature extraction is performed directly on the missing data.
• w/o Temporal (Temporal Domain): The temporal module in the conditional feature extraction block is removed, leaving only the frequency domain module.
• w/o Frequency (Frequency Domain): The frequency module in the conditional feature extraction block is removed, leaving only the temporal module.
• w/o Cross (Cross Attention): The cross-attention mechanism between temporal and frequency domains is removed and replaced by direct addition.

Table 4 presents the impact of each setting on the performance. We only analyze AQI-36 and METR-LA due to the computational cost of re-running diffusion-based CoFILL, especially on the larger PEMS-BAY. This aligns with the experimental setup of the baseline PriSTI [Liu et al., 2023]. The removal of forward interpolation produces the most significant performance degradation across datasets. For AQI-36, this modification increases the MAE from 8.7 to 9.15, representing a substantial 0.45 unit deterioration. The temporal domain component emerges as the second most critical element, with its removal causing more severe performance decline than eliminating the frequency domain component from the conditional information module described in Section 4.2. The cross-attention mechanism shows the least impact on model performance, particularly evident in the METR-LA dataset where its removal results in minimal MAE changes.

5.7 Hyperparameter Sensitivity Analysis

This experiment analyzes two important parameters in the diffusion model: max noise level (βT ) and channel size of the hidden state (d) that represent fundamental trade-offs. The former controls the intensity of the noise added during the diffusion process. A high noise level (βT ) may lead to distorted generated data, while a low noise level (βT ) may make it difficult to learn the various relationships within the data. The latter determines the model's ability to represent and handle information at each time step. A larger d allows the model to learn more features, but it may also lead to higher overhead.

The experimental results are illustrated in Figure 3 (only for AQI-36 and METR-LA due to the same reasons discussed in Section 5.6). For AQI-36, the model maintains stable MAE and MSE values with maximum noise level βT between 0.2 and 0.4. The METR-LA dataset shows higher sensitivity to βT variations within this range. The channel size d produces opposing effects, including degrading performance on AQI36 (increasing d) but improving results on METR-LA. This behavior reflects the different complexities in air quality and traffic data. Furthermore, we set d to 16 for AQI-36 and 64 for METR-LA. We apply METR-LA's parameter to PEMSBAY due to their similar data characteristics, as detailed in Table 1. This enables CoFILL to outperform PriSTI across most scenarios, which obtains better performance even under PEMS-BAY's Point where the metrics remain close to PriSTI.

6 Conclusion

In this paper, we present CoFILL, a novel conditional diffusion framework for spatiotemporal data imputation. CoFILL first overcomes the error accumulation problem by implementing a non-recursive processing structure through our diffusion framework. Second, it resolves the challenge of capturing complex temporal dependencies by introducing a dualstream architecture that simultaneously processes temporal and frequency domain features. Extensive experiments on real-world datasets demonstrate CoFILL's superior performance. Our model achieves better imputation accuracy than SOTAs across various missing data scenarios. We will develop adaptive fusion and extend CoFILL to handle multiresolution spatiotemporal data in the future.

Acknowledgments

This work is supported by the National Natural Science Foundation of China (62302148), Hebei Natural Science Foundation (F2024202076), and Hebei Yanzhao HuangJintai Talents Program (Postdoctoral Platform) (B2024003002).