Spatiotemporal Transformer Networks for Reconstructing Historical Landsat Time Series

1. Introduction

The Landsat program provides one of the most valuable long-term Earth observation records, delivering more than five decades of moderate-resolution imagery suitable for monitoring environmental change [1]. Its radiometric consistency, standardized surface reflectance products, and open data policy have established Landsat as a foundational dataset for land cover monitoring, ecosystem analysis, and change detection studies. The continuity of the archive enables retrospective analyses of landscape dynamics across multiple decades [2,3]. Nevertheless, despite its historical depth and spatial reliability, a fundamental constraint is the limited revisit frequency, which, when combined with the inconsistent availability of cloud-free acquisitions, results in sparse and irregular time series. These limitations pose persistent challenges for applications requiring temporally dense and continuous observations [4].

Reconstructing missing observations in optical satellite time series is essential for enabling reliable multi-temporal analysis. While interpolation or smoothing methods perform well for filling the gaps in high temporal resolution datasets such as MODIS [5], the highly sparse and irregular clear-sky sampling of Landsat, caused by cloud cover, shadows, atmospheric effects, and sensor revisit, makes Landsat time series reconstruction considerably more challenging [6,7]. Harmonic-based approaches, particularly those implemented within the Continuous Change Detection and Classification (CCDC) framework, have been widely used for Landsat time series modeling and change detection over the past two decades [4,8]. However, because these models rely on a limited set of harmonic components to describe seasonal behavior, they often struggle to represent complex, non-stationary, or highly variable temporal signals [9]. These challenges motivate the development of more flexible modeling frameworks capable of learning irregular and non-linear temporal dynamics directly from observations.

Deep learning (DL) methods have shown substantial potential for improving remote sensing applications, including time series analysis [10]. However, important challenges remain when applying these techniques to irregularly sampled data [11]. Most sequence modeling architectures, particularly recurrent neural networks (RNNs), are primarily designed for regularly spaced observations and are not inherently suited to handle uneven temporal intervals. Within the remote sensing domain, DL-based research on irregular time series remains limited, and most existing studies apply these architectures for classification directly on the sparse, irregularly spaced observations rather than for reconstructing dense, temporally continuous reflectance trajectories. For example, the transformer-based classifier of Zhang et al. [12] and the DL multiple-classifier framework of Dou et al. [13]. Where reconstruction has been attempted, methods such as the forward–backward LSTM of Zhou et al. [11] and the masked spatial–temporal attention network of Zhang et al. [14] restore missing values only within the existing set of acquired Landsat images, so the output remains constrained to the sensor’s coarse, irregular acquisition dates rather than being resampled onto a temporally continuous, regularly spaced grid suitable for fine-scale multi-temporal analysis. These limitations underscore the need for DL frameworks specifically designed for reconstructing highly sparse Landsat time series.

To address this gap, our previous work investigated a closed-form continuous-depth neural network integrated within an RNN architecture (CFC-mmRNN) for historical Landsat time series reconstruction [15]. The CFC design retains the advantages of continuous-time modeling inspired by ordinary differential equations (ODEs) but avoids the computational burden associated with numerical solvers, resulting in a more efficient implementation [16]. Although originally developed for general sequence modeling, CfC can also be applied to remote sensing time series reconstruction. When evaluated against the harmonic-based Continuous Change Detection (CCD) method under different observation densities, CFC-mmRNN produced consistently more accurate and temporally consistent reflectance estimates across spectral bands. It maintained superior performance and better preserved the overall temporal trajectory even under moderately sparse sampling conditions. Nevertheless, its robustness diminished when the number of available clear-sky observations fell below approximately four per year, indicating reduced reliability in extremely data-scarce scenarios. This limitation is critical for long-term Landsat archives, where irregular and infrequent acquisitions are common, particularly in regions with frequent cloud cover.

To address this challenge, subsequent research explored the use of a 1-dimensional (1D) Transformer-based network [17] for Landsat time series reconstruction [18]. Unlike CFC-mmRNN architecture, the Transformer relies on self-attention mechanisms, allowing it to capture long-range temporal dependencies and dynamically weigh the importance of individual time steps [12]. These capabilities make it particularly well-suited for modeling irregular and sparse time series data. In the study, the 1D Transformer was applied to reconstruct Landsat time series in a sample area of the Canadian Prairies from 1985 to 2023, and its performance was compared with that of the CFC-mmRNN method [18]. The model was evaluated across multiple spectral bands, seasons, and varying observation densities. Results showed that the 1D Transformer consistently outperformed CFC-mmRNN under low-density conditions, preserving the temporal structure of the reconstructed time series even when data were reduced to an average of only two clear observations per year. Both models performed comparably under high-density scenarios, highlighting the 1D Transformer’s robustness in sparse-data situations. These findings underscore the importance of selecting deep learning architectures that are well-suited to the specific nature of irregular and sparse time series remote sensing data.

Despite these advances, the 1D Transformer relies solely on temporal information, ignoring the spatial context. Incorporating spatial information from neighboring pixels can improve reconstruction accuracy, especially in heterogeneous landscapes, cloud- or shadow-contaminated areas, and post-SLC failure Landsat 7 imagery with linear data gaps [19]. To address this limitation, the present study extends the Transformer-based approach to a three-dimensional spatiotemporal architecture, integrating spectral, spatial and temporal information to improve Landsat time series reconstruction. Specifically, the framework incorporates a convolutional encoding mechanism, with spatial-spectral features extracted via a ResNet-based convolutional encoder and the encoded sequence processed by the Transformer network. Self-attention in this architecture operates only over the temporal dimension, following spatial-spectral feature extraction by the convolutional encoder; the model therefore does not perform joint spatiotemporal self-attention as in three-dimensional vision Transformers. The term 3D Transformer is retained hereafter to denote the three-dimensional (spatial, spectral, and temporal) structure of the input processed by this hybrid CNN–Transformer architecture, rather than the internal attention mechanism itself.

This study develops and evaluates a 3D Transformer-based network designed to handle irregular and sparse time series for reconstructing historical Landsat data from 1985 to 2023 over a sample area in Canada’s Prairie region. It investigates whether incorporating spatiotemporal dependencies can yield more accurate and robust reconstructions than models relying solely on temporal information, particularly under challenging conditions of very sparse observations, complex landscapes, or in the presence of SLC-induced striping in Landsat 7 imagery. A flowchart summarizing the methodology is presented in Figure 1.

2. Study Area and Dataset

The study was conducted in southeastern Alberta, Canada, within the Canadian Prairies, part of the northern Great Plains. The region is dominated by grasslands and agricultural lands, with elevations ranging from approximately 953 to 1,083 m above mean sea level. It was selected for its pronounced seasonal agricultural dynamics driven by crop phenology, harvesting cycles, and climate variability. These temporal variations provide a suitable framework for evaluating time-series reconstruction methods, particularly for distinguishing native and cultivated areas based on their seasonal behaviour.

To define evaluation sites, the region was divided into 10 × 10 km grid cells, from which three cells were randomly selected as test areas (Figure 2). To reduce spatial autocorrelation, training samples were extracted at least 2 km away from these tests locations. The dominant land-cover types in the first two test areas are cropland and grassland, while the third area is mainly characterized by forest, shrubland, and grassland (Figure 3).

This study uses atmospherically corrected and orthorectified surface reflectance imagery from Landsat 5 (TM), Landsat 7 (ETM+), Landsat 8 (OLI), and Landsat 9 (OLI-2), covering the period 1985–2023 and accessed via Google Earth Engine. The temporal coverage and operational timelines of these Landsat missions are illustrated in Figure 4. For Landsat 7, imagery acquired after the 2003 Scan Line Corrector failure was included; however, only data up to August 2017 were used to avoid potential inconsistencies caused by orbital drift [20]. The analysis focuses on spectral bands common to all sensors in the visible, near-infrared (VNIR; 0.43–0.90 µm), and shortwave infrared (SWIR; 1.55–2.35 µm) regions at 30 m spatial resolution. Thermal bands were excluded.

To support model training, two supplementary satellite datasets were incorporated. The first dataset consists of the MODIS 8-day surface reflectance products (MOD09A1.061 from Terra and MYD09A1.061 from Aqua) at 500 m spatial resolution, covering the 2004–2023 period [21,22]. The second dataset is the Harmonized Landsat–Sentinel-2 (HLS) product (version 2), which provides radiometrically and spatially consistent surface reflectance observations from the Operational Land Imager (OLI) onboard Landsat 8 and the MultiSpectral Instrument (MSI) onboard Sentinel-2A at 30 m spatial resolution for the 2017–2023 period [23]. For both datasets, six spectral bands, including blue, green, red, near-infrared (NIR), shortwave infrared 1 (SWIR1), and shortwave infrared 2 (SWIR2), were used in training the deep learning models.

3. Methodology

3.1. Transformer for Time Series Modeling

Originally developed for natural language processing [24], Transformer-based architectures have been increasingly applied to time series modeling due to their ability to handle sequences of varying lengths and capture long-range dependencies [18]. Their core component, the self-attention mechanism, enables dynamic weighting of different time steps, allowing the model to identify complex temporal patterns. By incorporating multi-head attention and positional encoding, Transformers can model temporal relationships without relying on recurrent structures [12]. The transformer architecture consists of several key components, including a masking layer, input and positional embeddings, stacked Transformer encoder blocks, and an output layer. The overall architecture of the Transformer-based time series reconstruction framework is depicted in Figure 5.

As the flowchart shows, reflectance values at each timestep are first processed through a masking layer, which handles missing observations by ensuring that the model focuses solely on valid data. The input embedding layer encodes each timestep’s predictors into a fixed-length feature vector, to which sinusoidal positional encoding is added to capture temporal order. Unlike simple Day-of-Year representations, this encoding generates position-dependent sine and cosine values that reflect the sequential structure of the time series. The Transformer comprises three stacked Transformer encoder blocks, each consisting of multi-head self-attention and feedforward layers, with four attention heads and a model embedding dimension of 128. Self-attention allows the model to weigh the relevance of different timesteps dynamically, improving its ability to preserve complex temporal patterns [17]. The position-wise feed-forward network within each encoder block has a hidden dimension of 512 and uses the ReLU activation function. To mitigate overfitting, a dropout rate of 0.01 was applied throughout the network, randomly disabling neurons during training. Layer normalization and residual connections are also applied to stabilize learning and enhance convergence. The final output layer is a fully connected dense layer that reconstructs the time series by predicting six values at each timestep, corresponding to the reflectance of the Landsat spectral bands (Blue, Green, Red, NIR, SWIR1, SWIR2); a sigmoid activation function constrains the predictions to the normalized valid range.

3.2. 3D (Spatiotemporal) Transformer for Time Series Modeling

To capture both spatial and temporal dependencies in Landsat image time series, the Transformer framework is extended with a convolutional encoding mechanism. Each input consists of a sequence of image patches (image cubes) with dimensions T×6×8×8, where T is the number of time steps, 6 is the number of spectral bands, and 8 corresponds to the window size. The 8×8 patch size was determined empirically through trial and error: smaller patches reduced reconstruction accuracy, while larger patches produced no further improvement.

In the 3D Transformer, spatial-spectral features are extracted from each image in the time series independently using a convolutional encoder applied via a Time Distributed layer. The encoder begins with a 3×3 convolutional layer with 32 filters, followed by three residual blocks based on the basic ResNet architecture [25], with increasing channel depth from 32 to 64 filters, interleaved with two max-pooling operations. Each residual block includes convolutional layers, batch normalization, ReLU activations, dropout for regularization, and skip connections that preserve low-level information while enabling deeper feature extraction. Spatial-spectral features are aggregated through global average pooling and projected to a 128-dimensional latent space to match the Transformer embedding size. The overall 3D Transformer framework is summarized in Figure 6, and the convolutional encoder and residual block architectures are illustrated in Figure 7 and Figure 8, respectively.

In parallel to the encoded patch sequence, the central pixel of each patch is extracted. Unlike the 1D Transformer, the masking condition in the 3D Transformer considers both the validity of the central pixel and the proportion of valid pixels within each patch. Specifically, a time step is considered valid and retained for further processing if either the central pixel is valid or at least 30% of the pixels in the corresponding patch contain valid observations. This threshold was selected empirically through trial and error during training. This formulation allows the model to leverage spatial context even when the central pixel is missing, while still filtering out severely corrupted patches. Consequently, although the 3D Transformer utilizes more data through spatial aggregation, its masking mechanism remains quality-aware by excluding time steps with insufficient valid information.

Both the encoded patch sequence and the central pixel time series are projected to the Transformer input dimension via fully connected embedding layers. To capture temporal order, a positional encoding is added to the patch embedding. The embeddings are then fused through a learned weighted sum:

$$\mathrm{Transformer}\mathrm{Input}=W_1\cdot \mathrm{EmbeddedPixel}+W_2\cdot \mathrm{EmbeddedPatch}$$

where $W_1$ and $W_2$ are trainable scalar weights, EmbeddedPixel is the embedded central pixel time series, and EmbeddedPatch represents the embedded spatial-spectral features. This design emphasizes the contribution of the central pixel while retaining the broader spatial context, and the learnable weights can be monitored during training to evaluate the relative importance of each component.

Importantly, after the embedding and weighted sum steps, a batch normalization layer is applied to the fused embeddings prior to the stack of Transformer encoder blocks, which are identical for both the 1D and 3D models. Within each encoder block, layer normalization is applied after the self-attention and feedforward sublayers, following standard Transformer design. This allows for a consistent comparison between models, isolating the contribution of the spatial encoder. This spatiotemporal design enables effective reconstruction of Landsat time series while capturing both local spatial features and global temporal dependencies.

3.3. Landsat Time Series Reconstruction: Training and Evaluation

The transformer-based frameworks for Landsat time series reconstruction were implemented through a structured multi-stage pipeline. First, cloud- and shadow-affected observations were identified and excluded from the Landsat archive to ensure data reliability (section 3.2.1.). Subsequently, target observations (reference time series) were generated using two distinct strategies. In the first strategy, combined Landsat and MODIS observations were modeled using Gaussian Process Regression (GPR) [26] to generate dense and continuous target time series at 500 m spatial resolution. In the second scenario, holdout time steps with valid observations from the HLS time series were used as target data. In the first scenario, the primary objective was to emphasize the temporal dynamics of the time series and train the models to accurately reconstruct the signal’s detailed temporal trajectory and structural pattern. In contrast, the second scenario helped to train the model for the higher 30 m spatial resolution, but reduced temporal detail with maximum accuracy, given that the HLS dataset exhibits high spatial and radiometric consistency with the original Landsat time series. These two target sets were used as supervisory targets for training the transformer models. After optimization, the trained networks were applied to estimate missing Landsat observations over the independent test areas within the study region. Models were systematically evaluated across spectral bands, seasonal conditions, land covers, observation densities, and temporal segments. The subsequent subsections describe each component of this workflow in detail.

3.3.1. Contaminated Observation Filtering

A key step in the preprocessing of Landsat imagery in this study was the removal of clouds and cloud shadows, as reliable reconstruction of Landsat time series depends heavily on careful data screening prior to model training. This was accomplished by applying the commonly used Fmask object-based algorithm [27] to mask affected pixels. While Fmask generally produces reliable cloud and shadow masks, it is not without limitations. For instance, it can sometimes confuse clouds with other transient features such as dense aerosols, smoke plumes, or flood-affected areas, potentially leading to misclassification. Nevertheless, given that Landsat time series inherently contain some level of noise and outliers, robust reconstruction methods must be capable of accommodating these uncertainties while preserving the underlying temporal patterns.

3.3.2. Training Samples Preparation

High-quality training samples are essential for deep learning models. In this study, 5,000 spatiotemporal image cubes were randomly extracted from the training region (blue area in Figure 2). The majority (~85%) of these were grassland and cropland, which are the dominant land covers in the study area. Sampling was stratified using the 2020 North American Land Cover 30-meter dataset available on Google Earth Engine [28].

The extracted Landsat samples were temporally sparse and irregular. To train the transformer to reconstruct dense and regularly spaced time series, corresponding dense reference data were required. To address this, in the first training strategy, Landsat and MODIS observations spanning 2004–2023 (20 years) were combined and modeled using GPR, a robust method used for Earth observation time series [29]. MODIS data were included as the daily coverage increases the likelihood of clear-sky observations. To maintain consistency, Landsat data were resampled to 500 m to match the resolution of the MODIS surface reflectance product.

GPR is a non-parametric Bayesian method that defines a distribution over possible functions and adjusts its complexity according to the data using Bayesian inference [30]. This flexibility enables GPR to capture a wide range of patterns and dependencies, making it particularly suitable for interpolating irregular time series while preserving trends and uncertainties [31]. By operating in an effectively infinite-dimensional function space, GPR can adapt to varying levels of smoothness, providing robust reconstruction of missing observations. This approach was used to produce high-quality, continuous, and evenly spaced target series suitable for training.

Because the quality of the reference time series directly influences the training procedure and overall model performance, generating reliable reference data is essential. Although GPR can be robust to moderate noise and outliers by selecting the right covariance functions, remaining anomalies in the Landsat and MODIS observations can still affect the reconstructed dense and continuous time series. Therefore, the reference data were generated through a two-step procedure to ensure high quality training data.

In the first step, NDVI was computed for each time series and modeled using GPR. Observations that deviated significantly from the modeled NDVI trajectory were identified as outliers and removed. Specifically, outlier detection was based on the GPR predictive uncertainty: observations falling outside the confidence interval (predicted NDVI ± 3σ) were excluded. This NDVI-based filtering process was performed iteratively with three iterations to progressively refine the valid observations. In the second step, the remaining observations were treated as more reliable inputs, and each spectral band was modeled separately using GPR. Band-wise modeling was also conducted iteratively, using two iterations to further filter outliers and to minimize the influence of anomalous values.

Although MODIS data were incorporated to increase temporal density for GPR modeling, Landsat observations were prioritized because they represent the primary focus of this study. Accordingly, different weights were assigned during the GPR fitting process: Landsat observations were given a weight of 2, whereas MODIS observations were assigned a weight of 1. Seasonal variability was also accounted for by additional weighting. Accurate modeling during spring, summer, and fall is more important than in winter in the study area, as winter months are typically snow-covered and associated with higher noise and uncertainty. Furthermore, among the non-winter seasons, spring and fall are more challenging to model than summer due to transitional vegetation dynamics; consequently, more attention should be given to the spring and fall seasons. Based on these considerations, seasonal weights of 1, 2, and 3 were assigned to winter, summer, and spring–fall months, respectively.

A fundamental component of GPR is the kernel (covariance function), which defines the similarity between observations and governs the smoothness, periodic behavior, and overall structure of the modeled function [32]. The adopted kernel configuration in this study consisted of a periodic kernel, a coarse-scale RBF (Radial Basis Function) kernel, a fine-scale RBF kernel, and a noise kernel. The periodic kernel was employed to model annual seasonality, using a period of 365 days and a length scale of 120 to control the frequency and smoothness of seasonal variations. The coarse-scale RBF kernel was incorporated to capture long-term temporal trends, with a length scale of $5\times 365$, while the fine-scale RBF kernel accounted for intra-annual variability. The length scale of the fine-scale RBF was set to 30 for most bands, except for the SWIR bands and NDVI, where values of 20 and 50 were used, respectively. Finally, the noise kernel was included to account for random observational variability. The selection of kernel components and their characteristic scales was determined experimentally to best reflect the temporal behavior of the time series in the study area. Figure 9 illustrates an example of a grassland time series modeled using the GPR framework, showing Landsat and MODIS observations along with the modeled dense time series for the NIR and SWIR2 bands.

In the second training strategy, HLS time series from 2017 to 2023 (7 years) were directly employed as target data, enabling the model to learn from the Harmonized Landsat–Sentinel-2 dataset, which exhibits strong spatial and radiometric consistency with native Landsat surface reflectance observations. The use of HLS data instead of the original Landsat dataset as the training target was primarily motivated by its higher temporal density, which provides more frequent observations and facilitates more effective model training.

As a data augmentation strategy and to improve generalization across varying time-series lengths, image cube samples for both training strategies were segmented into 5-year intervals using a moving window with a fixed length of five years and a one-year stride. The model is therefore trained on these 5-year segments and produces reconstructed time series of the same duration. Under this configuration, each 20-year time series in the first training strategy was divided into 16 overlapping samples, while each 7-year HLS time series was divided into 3 overlapping samples. This resulted in 5000 × 16 (80,000) samples for the first training set and 5000 × 3 (15,000) samples for the second training set. Table 1 summarizes the key characteristics and configurations of each training strategy, along with the corresponding number of training samples.

To reconstruct time series longer than five years, a moving-window approach is applied, in which overlapping 5-year segments are processed sequentially, and their outputs are aggregated to form the full series. Except at the boundaries, only the three central years of each reconstructed segment are retained, as edge regions tend to be less reliable than the central portion. The choice of a 5-year segment length was determined empirically through trial and error, with this configuration yielding the best performance in terms of reconstruction accuracy and stability.

For model training, the central pixel in the image cube samples served as input to the 1D Transformer, while the full cube was used for the 3D Transformer. In the first strategy, Landsat observations used for training were upsampled to 500 m resolution, whereas in the second strategy, the original resolution was maintained. To account for varying observation densities, a portion of the input observations was randomly dropped during each epoch in both training strategies, with dropout ratios randomly selected between 10% and 90%. To reduce computational load, all the time series were resampled at 7-day intervals. Finally, the two training sets were combined, with 85% of the data used for training and 15% for validation.

3.3.3. Training

The models were trained using a weighted Huber loss [33] tailored for multiband time series data, in which each timestep is assigned a seasonally dependent weight. Winter months (December–February) were assigned a lower weight (0.5), while all other months were assigned a weight of 1.0, thereby reducing the influence of winter observations during training. Invalid values in the target time series, denoted by −2, were masked and excluded from the loss computation. This masking strategy is applied specifically to the second training dataset, which contains missing values, allowing the model to leverage only the available observations from the HLS dataset as reference data. The Huber loss improves robustness to outliers by treating small errors quadratically and large errors linearly, reducing the impact of extreme deviations. The loss is computed over all valid timesteps and spectral bands, normalized at the sample level, and subsequently averaged across the batch.

The described architectures (Section 3.1 and Section 3.2) resulted in 629,894 and 807,912 trainable parameters for the 1D and 3D Transformer models, respectively. Both models were trained with the Adaptive Moment Estimation (Adam) optimizer [34], using the same training settings to enable a fair comparison between the two architectures. Training was conducted with a batch size of 5 over 100 epochs, totaling 1,615,000 iterations. All training and inference processes were carried out in Python 3.11 using TensorFlow 2.15 on an NVIDIA A30 GPU (NVIDIA Corporation, Santa Clara, CA, USA).

3.3.4. Evaluation

The performance of the models was assessed through two distinct evaluations. In the first evaluation, 1,500 samples randomly selected from the test areas (red regions in Figure 2) were used as independent test data. The predictive accuracy of each model was then compared against the GPR-generated reference time series for these samples. In the second evaluation, a subset of Landsat images from test areas was reserved to enable an independent comparison of the Transformer-based models. Specifically, 15% of images from each test area were set aside, including only those with at least 60% clear observations per season to ensure reliable assessment. This threshold minimized the influence of clouds, shadows, haze, or cloud-adjacency effects, which can introduce artifacts and bias evaluation results. Based on this criterion, 141, 100, and 78 images were retained for the first, second, and third test areas, respectively. Variations in the number of eligible images across test areas are partly due to differences in cloud cover patterns but primarily reflect the distinct satellite acquisition paths, which influence revisit times, scene overlaps, and coverage frequency. Root mean square error (RMSE) was used in both evaluations as the primary metric for model accuracy. Because GPR is also used to generate one of the two training targets (Section 3.2.2), the first evaluation shows how closely each model reproduces the GPR-modeled reference rather than measured reflectance. Although the test samples come from independent areas not used in training, the reference itself remains a GPR estimate rather than an actual observation. The second evaluation avoids this limitation: the held-out Landsat images were not used, directly or indirectly, in GPR modeling or model training, and therefore provide a genuinely independent benchmark against actual observed reflectance.

For both evaluations, results were analyzed over the complete time series as well as for individual seasons: spring (March–May), summer (June–August), fall (September–November), and winter (December–February). To investigate the impact of missing observations due to clouds or shadows, 10% to 90% of the time steps were randomly removed, generating ten scenarios ranging from the original observation density to progressively sparser levels. The average annual number of clear-sky observations was also calculated to determine the minimum needed for reliable reconstruction. In the second evaluation, model performance was further analyzed across different landcover types and satellite sources. For visual inspection, predicted images, error maps showing deviations from the reference, and error histograms were generated to highlight model accuracies.

4. Results and Discussion

4.1. Time Series Reconstruction Performance (GPR-Based Evaluation)

This section presents results from the first evaluation strategy, comparing the 1D and 3D Transformer reconstructions of 1,500 resampled 500 m Landsat time series from the test areas against the GPR-generated reference (Table 2 and Figure 10). As noted in Section 3.2.4, these results indicate how closely each model reproduces the GPR-modeled reference rather than measured reflectance; Section 4.2 presents the independent evaluation against actual Landsat observations. Overall, the 3D Transformer consistently outperformed the 1D Transformer across all six spectral bands, achieving lower RMSE values for the test samples. Notably, the magnitude of improvement is relatively uniform across bands, with an average reduction in RMSE of approximately 10%. To further evaluate model performance under varying seasonal conditions, reconstruction accuracy across different seasons is summarized and illustrated in Figure 11. The results indicate that the 3D Transformer consistently achieves higher accuracy than the 1D model across all seasons; however, the magnitude of the improvement varies seasonally. Overall, the largest performance gains are observed in spring, followed by fall and winter, whereas improvements during summer are minimal, suggesting that both models perform comparably for this season. This behavior indicates that when observation density is high and temporal variability is relatively low, the two approaches yield similar results as expected. In contrast, under conditions of sparse observations and increased temporal complexity, the 3D Transformer demonstrates a clear advantage.

Figure 12 presents the reconstructed red-band time series for a representative grassland pixel over the period 2007–2012, as produced by the 1D and 3D Transformer models, along with the corresponding GPR-based reference time series. Although both approaches exhibit comparable performance over much of the time span, the 3D Transformer provides noticeably more accurate reconstructions during intervals characterized by sparse observations, particularly in the spring of 2008 and 2010 (highlighted by the orange circles in Figure 12). This improvement can be attributed not only to the enhanced modeling capacity of the 3D architecture, but also to its ability to leverage spatiotemporal context by incorporating information from neighboring pixels across adjacent months.

As described earlier, model performance under varying observation densities was evaluated by randomly removing between 10% and 90% of the available observations and reconstructing the resulting time series using each method. Given the similar behavior observed across spectral bands and the number of plots required to show all band and season combinations, only the NIR-band results are presented in Figure 13. Consistent with earlier findings, the smallest performance gains are observed during summer, when the time series are denser and exhibit lower temporal complexity. While increasing the dropout rate leads to higher RMSE values for both models, this degradation is more pronounced for the 1D Transformer. Moreover, the relative advantage of the 3D Transformer becomes increasingly evident as observation density decreases, indicating that the 3D architecture is more robust to sparse temporal sampling compared with the 1D approach.

4.2. Image Reconstruction Performance (Hold-Out Landsat Evaluation)

This section presents the results of Landsat image reconstruction based on the Landsat holdout images. As described in the methodology, the 3D Transformer incorporates both the encoded patch sequence and the central pixel time series. The relative contribution of these two components is controlled by trainable scalar weights learned during model optimization. After training, the optimized weights were 55.8% for the central pixel component and 44.2% for the patch sequence component. These values indicate that although temporal information contributes slightly more to the reconstruction process, spatial context also plays a substantial role, confirming the importance of integrating spatiotemporal dependencies for improved time series reconstruction.

Table 3 and Figure 14 summarize the reconstruction accuracy using the hold-out test images from the test areas, quantified using RMSE, for both models. The results indicate that the 3D Transformer consistently outperformed the 1D model across all three test areas. The highest reconstruction accuracy was observed for the SWIR bands, followed by the NIR, Red, Green, and Blue bands, and this performance ranking remains consistent across all study areas. The improvement rates, as reported in Table 3, vary across spectral bands and study areas, ranging from 3% to 15%.

A more detailed evaluation was conducted to compare model performance across the three test areas under various conditions, including different seasons, land-cover types, time periods, and observation densities. Although minor variations are observed among the test areas, the overall performance patterns of both Transformer models remain consistent, in agreement with the findings reported in previous analysis [18]. Consequently, given the large number of experimental scenarios and corresponding figures, the remaining analyses are summarized using performance metrics averaged across the three test areas. Exceptions are made for the single-image reconstruction examples (Figure 18 and Figure 20) and the spatial error analysis (Figure 19), which are presented for the first and third test areas.

Figure 15 presents the average reconstruction accuracy of the test images across different seasons and spectral bands. Consistent with the time series reconstruction results discussed in Section 4.1, the 3D Transformer achieves higher reconstruction accuracy than the 1D Transformer in all seasons. Summer is only slightly improved if at all. In agreement with our previous study [18], differences between models are most evident during transitional seasons, when observation density is lower and temporal dynamics are more complex due to rapid reflectance changes in vegetated landcover types relative to summer and winter conditions. This pattern is further confirmed by the average improvements across all spectral bands: the 3D Transformer improves image reconstruction accuracy by 29%, 4%, 20%, and 6% in spring, summer, fall, and winter, respectively.

Figure 16 summarizes the average test image reconstruction accuracy, measured by RMSE, across different land-cover types. Since the study area is generally covered by snow during the winter months, including winter test images introduces additional noise and variance, which often leads to substantially higher errors. Therefore, only test images from non-winter months were used to analyze and compare the accuracy. The results indicate that both models achieve the highest reconstruction accuracy for forested areas, followed by shrubland, grassland, and cropland. This pattern is expected, as approximately 75% of the forested pixels correspond to needleleaf (evergreen) species, with the remaining portion dominated by broadleaf deciduous vegetation, resulting in relatively low temporal variability within this class. In contrast, croplands exhibit the largest reconstruction errors, reflecting their higher temporal variability and more complex phenological dynamics compared with the other land-cover types. This complexity stems from the diversity of crop types and rotations across fields, each with distinct planting and harvest schedules that produce abrupt, management-driven reflectance transitions rather than the gradual seasonal changes characteristic of natural vegetation.

In a subsequent analysis, we evaluated the capability of each model to reconstruct Landsat time series for SLC-on and SLC-off periods. The SLC-on period was from January 1997 to December 2009, using Landsat 5 observations. The SLC-off period was from January 2004 to December 2016, using Landsat 7 imagery acquired after the scan line corrector (SLC) failure on May 31, 2003 (see Figure 4). It is also worth noting that the average number of clear observations per year was approximately 12.5 for the SLC-on period and 9.6 for the SLC-off period. This experiment was primarily designed to evaluate the robustness of each model under SLC-off conditions and to assess its ability to mitigate data gaps caused by the Landsat 7 sensor failure.

As shown in Figure 17, the 3D Transformer consistently achieves higher reconstruction accuracy than the 1D model under both scenarios. Both approaches perform better during the SLC-on period, which is expected given the absence of striping artifacts and the greater number of clear observations compared with the SLC-off scenario. Notably, the relative performance gain of the 3D Transformer is larger during the SLC-off period, highlighting its superior capability under more challenging conditions such as the presence of scan-line artifacts. These findings suggest that the 3D Transformer is more effective than the 1D model in mitigating the adverse effects of striping associated with the SLC-off condition.

Figure 18 also presents two examples of Landsat image reconstruction using the 1D and 3D Transformer models under SLC-on and SLC-off conditions for the first test area. The SLC-on example corresponds to an image acquired by Landsat 5 on 3 July 2005, while the SLC-off example was acquired by Landsat 7 on the following day, 4 July 2005. As shown, both models achieve comparable and satisfactory performance for reconstructing the SLC-on image. In contrast, for the SLC-off image, the reconstruction produced by the 1D Transformer is noticeably affected by scan-line striping artifacts. Although residual striping is also present in the image reconstructed by the 3D Transformer, this model mitigates the SLC-off effects more effectively than the 1D approach.

Figure 18. Image reconstruction results obtained using the 1D and 3D Transformer models for two test images acquired by Landsat 5 (SLC-on) and Landsat 7 (SLC-off), illustrating the impact of scan-line striping on model performance.

Figure 18. Image reconstruction results obtained using the 1D and 3D Transformer models for two test images acquired by Landsat 5 (SLC-on) and Landsat 7 (SLC-off), illustrating the impact of scan-line striping on model performance.

Figure 19 presents the average error maps and corresponding histograms, showing the absolute differences between the reference and reconstructed values for 78 test images from path 39 and row 25 in the third test area generated by the 1D and 3D Transformer models. The results indicate that the 3D Transformer consistently had lower reconstruction errors across all six spectral bands, as evidenced by both the spatial error patterns and the associated error distributions.

A spatial contrast is observed between land-cover types, with forested areas, predominantly located in the southern portion of the test area, exhibiting lower reconstruction errors than grassland and shrubland regions, which are mainly concentrated in the middle and northern part. This pattern is consistent with the results shown in Figure 16 and is visually apparent in the error maps, where forested regions are dominated by lower error values. The error distributions further support this observation by revealing two distinct peaks, with the lower peak corresponding to forested areas. This discrepancy can be largely attributed to the reduced temporal variability and lower complexity of the time series associated with coniferous forests compared with grasslands and shrublands. Similar behavior has been reported in previous studies [15,18], where croplands exhibited higher reconstruction errors than grasslands due to stronger reflectance fluctuations and more complex temporal dynamics, which increase the difficulty of accurate reconstruction.

Figure 19. Average error maps and corresponding histograms of error values for reconstructed test images generated using the 1D Transformer (left) and the 3D Transformer (right).

Figure 19. Average error maps and corresponding histograms of error values for reconstructed test images generated using the 1D Transformer (left) and the 3D Transformer (right).

To visually assess the performance of each method further, Figure 20 presents examples of reconstructed test images across different seasons using both the 1D and 3D Transformer models for the NIR band predictions. Based on the reconstructed images alone, both models appear to achieve nearly identical reconstruction performance. However, the error maps reveal a slightly better performance of the 3D Transformer compared to the 1D model. This difference is particularly minimal during summer, where the error maps indicate highly comparable performance between the two models.

Figure 20. Example image reconstruction and error maps using 1D and 3D Transformer for the NIR band across seasons.

Figure 20. Example image reconstruction and error maps using 1D and 3D Transformer for the NIR band across seasons.

In the following, the average results of Landsat image reconstruction for the NIR band are presented across all three study areas under different observation density levels (Figure 21). As described earlier, between 10% and 90% of the available observations were randomly selected and removed, and the test images were reconstructed using the remaining observations to analyze the effect of density level on the performance of each model. Overall, the 3D Transformer outperformed the 1D Transformer in reconstructing test images across all seasons and most density levels. Consistent with previous findings, the largest performance improvements were observed in spring and fall, followed by winter and summer. Both methods exhibited similar performance during the summer. However, at high dropout rates (80% and 90%), the 1D model was slightly better. The results also indicate that the increase in RMSE with increasing dropout rates follows a similar trend for both methods, with a pronounced rise in error once more than 60% of the observations were removed, corresponding to fewer than five clear observations per year on average. Similar results were observed for the other bands.

Examining the dropout results for an example single pixel time series of cropland is shown in Figure 22. At high density levels, both methods show comparable performance in predicting test observations and reconstructing dense, regular time series. However, as the observation density decreases, the 3D model predicts test images more accurately, particularly for observations corresponding to transitional seasons (spring and fall).

For instance, at an 80% dropout rate, the 3D model clearly outperforms the 1D model in predicting test observations around spring 2010, spring 2012, fall 2013, and spring 2016, as highlighted by the orange circles. Furthermore, at a 90% dropout rate, which corresponding to fewer than two clear observations per year on average, the 3D model reconstructs the time series more accurately than the 1D model, especially during the 2010–2013 period, when the availability of clear observations is particularly limited. At this extreme sparsity level, despite the absence of clear winter observations in 2010, 2011, and 2012, the 3D model successfully preserves the overall temporal shape of the time series. This robustness can be attributed to its use of spatiotemporal contextual information from neighboring pixels, rather than relying solely on the spectral information of the central pixel.

4.3. Comparison of 1D and 3D Transformer for Computational Performance

In this section, the processing time and computational burden of each method are evaluated. Using the input data described in Section 3.2.2., training the 1D model, with 629,894 parameters, required 237 seconds per epoch (approximately 6.6 hours per 100 epochs). In contrast, the 3D model, with 807,912 parameters, required 439 seconds per epoch (approximately 12.2 hours per 100 epochs). Thus, the 1D model is approximately 1.85 times faster than the 3D model during training.

In the prediction phase, however, the processing times of the two methods are identical. Both models reconstruct a single pixel time series with a length of 2,028 weeks (39 years) in 167 milliseconds. Consequently, reconstructing test images of size 310 × 310 × 2,028 × 6 (rows × columns × time steps × bands) requires 268 minutes (~4.47 hours) for either method, considering prediction time only.

Based solely on training time and computational cost, the additional burden of the 3D model may appear moderate, particularly given its overall superior reconstruction performance compared to the 1D model. However, computational cost is not the only practical consideration; data preparation constitutes an even more critical factor. Training the 1D Transformer requires only 2D time-series inputs (time steps × six bands), whereas the 3D Transformer requires 4D image chips (rows × columns × time steps × bands) for each sample. Preparing, downloading, and storing image chips is substantially more demanding than exporting single-pixel time series from platforms such as Google Earth Engine (GEE). For comparison, 5,000 input time-series samples stored in npz format (a standard NumPy storage format) for training and testing the 1D model occupy approximately 306 MB. In contrast, 5,000 image chips required for the 3D model occupy approximately 16,370 MB. Therefore, the training dataset for the 3D model is about 53.5 times larger than that of the 1D model. This substantial difference can lead to serious storage and data transfer limitations, particularly when training large-scale or global models that require hundreds of thousands of samples.

An important distinction of the 3D model in this study is that it operates as an image-sequence-to-point-sequence framework, where spatial context is used to reconstruct the time series of the central pixel. This design implies that, if extended to predict full image sequences (image patches) instead of individual pixel time series, the number of required model inferences could be significantly reduced. Specifically, for a chip size of 8 × 8, predictions could be generated for multiple pixels simultaneously, potentially reducing the total number of inferences by a factor proportional to the patch size. In contrast, the 1D model processes each pixel independently, requiring a separate inference for every location within the area of interest. However, generating full image-sequence predictions was not explored in this study, primarily due to challenges associated with constructing reliable target image sequences for training. In particular, residual noise sources such as clouds, shadows, haze, geolocation inconsistencies, and atmospheric effects may propagate into the target data and affect model performance. These challenges become more critical under low observation densities, where errors in both training and evaluation data can have a larger impact on reconstruction accuracy. Therefore, while the 3D framework shows clear potential for more efficient large-scale inference, further investigation is required to develop high-quality spatiotemporal training datasets and to assess their influence on model performance.

5. Conclusions

This study introduced a spatiotemporal Transformer framework for reconstructing historical Landsat time series and demonstrated that incorporating spatial context substantially improves reconstruction accuracy, robustness, and stability compared with a purely temporal 1D Transformer. Across all evaluations, including GPR-based time series reconstruction, hold-out Landsat image reconstruction, seasonal analyses, land-cover comparisons, and SLC-on/off scenarios, the 3D Transformer consistently outperformed the 1D model. The GPR-based results indicate how closely each model reproduces the GPR-modeled reference, while the hold-out image results confirm that this advantage extends to actual measured reflectance. The largest gains occurred during spring and fall, when phenological transitions produce rapid reflectance changes and observation density is typically low. Under these challenging conditions, the 3D model reduced the RMSE on average across the six spectral bands by 29% and 20% for image reconstruction during the spring and fall seasons, respectively. In addition, it preserved temporal structure even when fewer than two clear observations per year were available.

The 3D Transformer also demonstrated superior resilience to Landsat 7 SLC-off striping, producing smoother and more coherent reconstructions than the 1D model. Spatial error analyses confirmed that the 3D architecture reduces reconstruction errors across heterogeneous landscapes, particularly cropland areas where temporal variability was high. Although the 3D model requires substantially larger training datasets and longer training times, inference costs are identical for both models, and the 3D framework offers the potential for more efficient large-scale prediction if extended to image-sequence reconstruction.

Overall, the results demonstrate that combining convolutional spatial-spectral feature extraction with temporal self-attention provides a powerful and flexible approach for reconstructing long-term Landsat time series. By effectively capturing both local spatial structure and long-range temporal dependencies, the 3D Transformer enhances the reliability of reconstructed reflectance trajectories. This should enhance land-cover monitoring, change detection, and environmental analysis, particularly for periods with low-density observations. The model’s superior performance during transitional seasons, when phenological dynamics are most complex, suggests strong potential benefits for downstream applications that rely on spring and fall observations as key discriminatory periods. Future research will first examine the practical impact of these seasonal improvements on land-cover classification and change-detection workflows. A second direction will focus on developing high-quality spatiotemporal training datasets and extending the 3D framework toward full image-sequence prediction to enable more efficient large-scale or global implementations.