Global 30 m seamless data cube (2000–2022) of land surface reflectance generated from Landsat 5, 7, 8, and 9 and MODIS Terra constellations

The Landsat series constitutes an unparalleled repository of multi-decadal Earth observations, serving as a cornerstone in global environmental monitoring. However, the inconsistent coverage of Landsat data due to its long revisit intervals and frequent cloud cover poses significant challenges to land monitoring over large geographical extents. In this study, we developed a full-chain processing framework for the multi-sensor data fusion of Landsat 5, 7, 8, and 9 and MODIS Terra surface reflectance products. Based on this framework a global 30 m resolution daily seamless data cube (SDC) of land surface reflectance was generated, spanning from 2000 to 2022. A thorough evaluation of the SDC was undertaken using a leave-one-out approach and a cross-comparison with NASA's Harmonized Landsat and Sentinel-2 (HLS) products. The leave-one-out validation at 425 global test sites assessed the agreement between the SDC with actual Landsat surface reflectance values (not used as input), revealing an overall mean absolute error (MAE) of 0.014 (the valid range of surface reflectance values is 0–1). The cross-comparison with HLS products at 22 Military Grid Reference System (MGRS) tiles revealed an overall mean absolute deviation (MAD) of 0.017 with L30 (Landsat 8-based 30 m HLS product) and a MAD of 0.021 with S30 (Sentinel-2-based 30 m HLS product). Moreover, experimental results underscore the advantages of employing the SDC for global land cover classification, achieving a sizable improvement in overall accuracy (2.4 %–11.3 %) over that obtained using Landsat composite and interpolated datasets. A web-based interface has been developed for researchers to freely access the SDC dataset, which is available at https://doi.org/10.12436/SDC30.26.20240506 (Chen et al., 2024).

1 Introduction

Earth observation (EO) data acquired by satellite sensors are fundamental to global land monitoring (Markham and Helder, 2012; Song et al., 2018; Wulder et al., 2022), providing critical information sources with unparalleled spatial and temporal coverage at a low cost. Over the past decades, satellite remote sensing has emerged as a prominent technology in Earth system science (Gong et al., 2023; Yang et al., 2013), contributing to the monitoring of land surface dynamics (Gong et al., 2013, 2019; Huang et al., 2017; H. Liu et al., 2020, 2021; Song et al., 2018), land surface phenology (Bolton et al., 2020; Piao et al., 2019), forest (Hansen et al., 2008, 2013), water (Ji et al., 2018; Pekel et al., 2016; Pickens et al., 2022; Sagan et al., 2020), and urbanization (Gong et al., 2012, 2020; X. Liu et al., 2020).

The Landsat series stands as the most enduring source of Earth observations, with a historical archive extending back to 1972 (Wulder et al., 2022). This longevity, combined with its relatively high spatial resolution, rigorous radiometric calibration, and free-access policy, has made Landsat a cornerstone for monitoring global terrestrial environments (Markham and Helder, 2012; Wulder et al., 2022). Nevertheless, the utility of Landsat data inevitably encounters certain limitations. A notable constraint is its relatively low temporal frequency, revisiting each area on Earth every 16 d (8 d when there are two Landsat satellites in orbit with an 8 d offset) (Zhu et al., 2015b). This issue is further compounded by the presence of cloud and cloud shadow, which can introduce significant temporal gaps in the acquisition of clear-sky observations, especially in cloudy regions (Zhu et al., 2016). Moreover, Landsat time series observations typically exhibit irregularities in both observation frequencies and acquisition dates due to the presence of cloud contamination and the geographically heterogeneous Landsat overpass coverage (Li and Roy, 2017). These irregularities present significant challenges when utilizing Landsat for large-scale monitoring of land cover and land use change (Potapov et al., 2020; Zhang et al., 2024). Therefore, the availability of Landsat datasets characterized by consistency in both temporal and spatial dimensions is crucial for facilitating various global environmental studies (Khan et al., 2024; Li et al., 2023; Pickens et al., 2022; Potapov et al., 2021a, b, 2022a, b; Song et al., 2021; Turubanova et al., 2023).

One conventional approach employed to mitigate data gaps in optical remote sensing is image compositing, which selects the highest-quality observations within a pre-defined time interval based on specific criteria to create seamless clear images at large scales (Jin et al., 2023; Qiu et al., 2023; White et al., 2014). Historically, image compositing has been mostly applied to coarse-resolution data with high temporal frequency (Qiu et al., 2023), such as that obtained by the Advanced Very High Resolution Radiometer (AVHRR) and Moderate Resolution Imaging Spectroradiometer (MODIS) sensors (Chuvieco et al., 2005; Cihlar et al., 1994; Holben, 1986; Huete et al., 2002; Wolfe et al., 1998). The use of image compositing for medium-resolution data (e.g., Landsat) was comparatively uncommon before the advent of Landsat free-access policy in 2008 (Qiu et al., 2023). In recent years, many image compositing algorithms have been developed for Landsat data (Frantz et al., 2017; Griffiths et al., 2019; Jin et al., 2023; Nelson and Steinwand, 2015; Qiu et al., 2023; Roy et al., 2010; White et al., 2014). Nevertheless, Landsat image compositing is not without its limitations. Firstly, due to the lack of frequent Landsat observations (especially in cloudy areas), it may take several months or even years to provide a composite Landsat image, which can cause problems if there are land cover or phenological changes (Zhu et al., 2015b). Furthermore, the compositing process may introduce distortions to the temporal dynamics of Landsat time series (Qiu et al., 2023), thereby hampering subsequent applications that depend on precise temporal information.

Landsat interpolation methods also provide the capability to generate seamless synthetic Landsat images (Brooks et al., 2012; Malambo and Heatwole, 2016; Yan and Roy, 2018, 2020; Zhu et al., 2015b). Linear interpolation is commonly employed to address missing values in Landsat time series (Defourny et al., 2019; Inglada et al., 2017; Tran et al., 2022), though it may not be highly effective for applications like land cover classification (Che et al., 2024). To improve performance, more sophisticated interpolation methods have been developed (Brooks et al., 2012; Malambo and Heatwole, 2016; Yan and Roy, 2018; Zhou et al., 2022; Zhu et al., 2015b). Nevertheless, a significant limitation of these methods is their dependence on numerous clear-sky Landsat observations for accurate time series estimation (Chen et al., 2021; Zhu et al., 2015b). This requirement poses a considerable obstacle to their large-scale applications, particularly in cloudy areas. Moreover, the performance of interpolation-based methods relies on the careful tuning and selection of model parameters, thereby encountering the challenge of balancing between over-fitting and under-fitting (Wu et al., 2022; Zhou et al., 2022). Large-scale remote sensing applications prefer processing algorithms capable of automatic adaptation to diverse input data conditions, eliminating the need for manual parameter tuning (Chen et al., 2023).

The spatiotemporal fusion technique provides a promising solution which aims at incorporating more frequent coarse-resolution observations to enhance the temporal frequency of Landsat and generate synthetic Landsat-like dense time series images (Chen et al., 2023, 2021; Gao et al., 2006; Liu et al., 2022; Zhu et al., 2010, 2016). For example, the Terra/Aqua Moderate Resolution Imaging Spectroradiometer (MODIS) provides frequent coarse-resolution observations at a 250, 500, and 1000 m spatial resolution with a near-daily revisit frequency (Schaaf et al., 2002). The MODIS land bands have comparable center wavelengths to the Landsat Enhanced Thematic Mapper Plus (ETM+) sensor, making the MODIS data the ideal input for the spatiotemporal fusion with Landsat (Gao et al., 2006). Many different types of Landsat–MODIS spatiotemporal fusion methods have been developed (Chen et al., 2023; Gao et al., 2006, 2022; Goyena et al., 2023; Guo et al., 2020; Hilker et al., 2009a; Liu et al., 2019, 2022; Mizuochi et al., 2017; Shi et al., 2022; Wang et al., 2017, 2020; Zhang et al., 2013; Zhu et al., 2010, 2016; Zurita-Milla et al., 2008) and applied to land cover and land surface phenology monitoring (Abowarda et al., 2021; Battude et al., 2016; Chen et al., 2018; Gervais et al., 2017; Hilker et al., 2009b; Y. Li et al., 2017; Senf et al., 2015; Singh, 2011; Tian et al., 2013; Walker et al., 2012; Watts et al., 2011). The utility of multi-sensor data fusion in facilitating land cover and land use analyses has also been validated empirically in previous studies (Carrasco et al., 2019; Chen et al., 2017; Yin et al., 2019).

To the best of our knowledge, there is currently no global 30 m seamless dataset of land surface reflectance generated by fusing Landsat and MODIS products available to the community. Although there have been numerous studies dedicated to developing algorithms for missing data reconstruction and multi-sensor data fusion (Shen et al., 2015; Zhu et al., 2018), unified and generalized frameworks for effective and efficient Landsat–MODIS data fusion on a global scale have not yet been explored extensively. To address this need, in this study, we (i) developed a full-chain processing framework for the multi-sensor data fusion of Landsat 5, 7, 8, 9 and MODIS Terra surface reflectance products; (ii) generated a global 30 m daily seamless data cube (SDC) of land surface reflectance, covering the period from 2000 to 2022; (iii) evaluated the reconstruction accuracy of the proposed framework quantitatively using a leave-one-out strategy at 425 global test sites; (iv) evaluated the quality of the SDC quantitatively by cross-comparing it with the Harmonized Landsat and Sentinel-2 (HLS) products at 22 Military Grid Reference System (MGRS) tiles; (v) evaluated the performance of using the SDC for global-scale land cover classification against Landsat composite and interpolated datasets; and (vi) provided a web-based interface for researchers to freely access the SDC dataset.

2 Materials

2.1 Landsat Collection 2 level-2 surface reflectance products

We collected a comprehensive dataset comprising 6 564 546 Landsat Collection 2 level-2 surface reflectance (L2SR) images from the US Geological Survey (USGS) Earth Resources Observation and Science (EROS) Center, including data acquired by the Landsat 8 and 9 Operational Land Imager (OLI), Landsat 7 Enhanced Thematic Mapper Plus (ETM+), and Landsat 5 Thematic Mapper (TM). This dataset covers most of global land surface except Antarctica, spanning from 2000 to 2022.

The L2SR products are generated through a sequence of processing steps applied to Landsat raw data. These steps include reprojection, radiometric calibration, geometric correction, atmospheric correction, and cloud masking. Compared to previous Landsat Collection 1 products, Collection 2 products have markedly improved the Landsat absolute geolocation accuracy using Landsat 8 geolocational imaging performance harmonized with the European Space (ESA) Agency Global Reference Image (GRI) data (Crawford et al., 2023). The Landsat 5 and 7 TM and ETM+ data were atmospherically corrected using the Landsat Ecosystem Disturbance Adaptive Processing System (LEDAPS) (Masek et al., 2006), and the Landsat 8 and 9 OLI data were corrected using the Land Surface Reflectance Code (LaSRC) (Vermote et al., 2016). The Fmask algorithm (Zhu and Woodcock, 2012) was applied to detect cloud and cloud shadow in Landsat images. The L2SR products are spatially referenced using Worldwide Reference System-2 (WRS-2) path rows and provided in Universal Transverse Mercator (UTM) projection. Figure 1 illustrates the spatial and temporal distribution of all L2SR images used in this study. As listed in Table 1, we used blue, green, red, near-infrared (NIR), and two shortwave infrared (SWIR1 and SWIR2) bands for the generation of the SDC. Although Landsat sensors have a relatively narrow field of view (15°), the bidirectional reflectance distribution function (BRDF) normalization was found to be effective in making multi-temporal Landsat observations more consistent (Claverie et al., 2015; Roy et al., 2016b). We applied the C-factor technique and global constant BRDF coefficients provided by Roy et al. (2016a) to obtain the Landsat nadir BRDF-adjusted reflectance (NBAR) dataset for SDC generation.

Figure 1 Spatial and temporal distribution of L2SR images used in this study.

Table 1 Attributes of Landsat 5 TM, 7 ETM+, and 8–9 OLI and MODIS Terra products (Markham and Helder, 2012; Masek et al., 2020; Morisette et al., 2002).

2.2 MODIS nadir BRDF-adjusted reflectance (NBAR) products

We reprocessed the MODIS MOD09GA version 6.1 surface reflectance product to obtain a daily seamless (without missing values) 500 m MODIS NBAR dataset as detailed in Liang et al. (2024). The official MODIS MCD43A4 NBAR product was not employed due to persisting concerns associated with it, including the prevalence of missing data and residual influence of cloud and aerosols (Liang et al., 2024).

All MODIS Terra Surface Reflectance MOD09GA version 6.1 images for the period 2000–2022 were acquired from NASA Earthdata. The MOD09GA product is tiled in the MODIS sinusoidal system with a spatial resolution of about 500 m. A set of processing algorithms was applied to these MOD09GA images to derive a daily 500 m resolution seamless MODIS NBAR data cube (Liang et al., 2024), including three main stages: (i) land-cover-based BRDF correction with the kernel-driven RossThick LiSparse reciprocal (RTLSR) model using parameters derived from the MCD43A1 BRDF model parameter dataset and land cover maps from the MCD12Q1 land cover product, (ii) outlier removal and gap filling using the ecosystem curve-fitting method, and (iii) sliding-window temporal smoothing using the Savitzky–Golay filter. The generated MODIS NBAR seamless dataset provides seven spectral bands that are commonly used for terrestrial applications, six of which that have compatible bandwidths with Landsat sensors are employed for SDC generation as listed in Table 1.

2.3 The Harmonized Landsat and Sentinel-2 V2.0 surface reflectance product

NASA's Harmonized Landsat and Sentinel-2 (HLS) V2.0 products were used in cross-comparison with the generated SDC product for quantitative assessment. HLS products combine observations from the Landsat Operational Land Imager (OLI; since 2013) and the Sentinel-2 Multi-Spectral Instrument (MSI; since 2015), providing global surface reflectance data at a 30 m spatial resolution with a theoretical revisit interval of 2–3 d at the Equator and even more frequent revisits in areas of higher latitudes (Claverie et al., 2018). The creation of HLS products involves four major processes (Claverie et al., 2018): (i) atmospheric correction and cloud masking, (ii) geometric resampling and geographic registration, (iii) BRDF normalization, and (iv) bandpass adjustment. HLS products are gridded into the UTM Military Grid Reference System (MGRS) used by the Sentinel-2 products. The HLS S30 product (Sentinel-2-based 30 m product) is derived from 10 and 20 m Sentinel-2 bands using overlapping-area-weighted averaging, and the L30 product (Landsat-based 30 m product) is reprojected to the same Sentinel-2 grid using cubic convolution interpolation (Claverie et al., 2018). Both HLS L30 and S30 products are atmospherically corrected using the Land Surface Reflectance Code (LaSRC) (Vermote et al., 2016). The cloud mask used in HLS products is a combination of the mask derived from the Fmask algorithm and the mask derived from the LaSRC algorithm. The HLS L30 and S30 products are delivered as NBAR, using the C-factor technique and global constant BRDF coefficients provided by Roy et al. (2016a). A bandpass adjustment is applied to the S30 product using a global constant set of coefficients (Claverie et al., 2018).

3 Methods

Figure 2 illustrates the overview of the SDC processing chain, comprising five key processing steps.

Figure 2 Overview of the SDC processing chain for each task unit.

3.1 Gridding and reprojection

The UTM-based Military Grid Reference System (MGRS) was chosen as the projection system for the SDC product and was also adopted by ESA's Sentinel-2 products and NASA's HLS products (Claverie et al., 2018). It is noteworthy that our adopted grid slightly deviates from the Sentinel-2 grid. The original Landsat coordinate system aligns the UTM coordinate origin with a pixel center, while the Sentinel-2 grid aligns it with a pixel corner (Claverie et al., 2018). We expanded the Sentinel-2 tiles by 15 m in each direction to align them with the original Landsat coordinate system, minimizing the need for resampling Landsat data. The SDC product is gridded onto this modified MGRS system, with a tile size of 109.83 km×109.83 km (3661×3661 Landsat pixels). Although it has been found that overlaps among neighboring MGRS tiles may result in resource wastage to some extent (Bauer-Marschallinger and Falkner, 2023), implementing this UTM-based projection system serves to minimize the need for resampling Landsat data, thereby reducing the introduction of additional errors (Dwyer et al., 2018).

The metadata for all Landsat and MODIS images has been pre-indexed into a database. For each SDC generation task unit, the metadata regarding all source data falling within specified spatial and temporal ranges can be efficiently retrieved from the database. Subsequently, all involved Landsat and MODIS source data are reprojected and gridded onto our UTM-based grid using nearest-neighbor resampling for Landsat and bilinear resampling for MODIS. MODIS data are resampled to 30 m spatial resolution to streamline subsequent processing steps, and the computational costs for this upscaling operation are negligible. Table 2 lists the input data products across distinct time periods utilized in SDC generation.

Table 2 Input product specifications across distinct periods.

3.2 Landsat cloud masking

Cloud and cloud shadow masks are essential for removing contaminated Landsat pixels in SDC generation. We used Fmask (Zhu and Woodcock, 2012) detection results as primary cloud and cloud shadow indicators. There are still a few clouds and heavy aerosols in Landsat images that remain undetected by the current Fmask method, which may have noteworthy implications for subsequent data processing procedures (Chen et al., 2021). To mitigate this issue, an enhanced cloud filtering approach is employed to reduce residual clouds and cloud shadows in Landsat imagery, which comprises three major steps:

• 1.

  Firstly, the cloud and cloud shadow masks generated by the Fmask algorithm are expanded by a margin of 150 m (Claverie et al., 2018). This dilation process is designed to exclude potentially contaminated pixels adjacent to the initially detected cloud and shadow areas.

• 2.

  Secondly, a brightness-threshold filter combined with a spatial filter is applied to remove remaining highly reflective pixels, which is achieved by cross-comparing with MODIS NBAR data. This filter operates on a patch-wise basis, with each patch measuring 20 × 20 Landsat pixels. For a given image patch, we commence by computing the ratio of Landsat reflectance (summed over the six spectral bands) to MODIS reflectance for each pixel. Following this, if the median of all these Landsat–MODIS reflectance ratios within this image patch surpasses a pre-defined threshold, the entire image patch is then flagged as cloudy. The threshold was set to 2 in this study as this value effectively eliminates most residual clouds without being overly aggressive.

  $$\begin{array}{}\text{(1)}& \mathrm{\text{median}}\left\{{\displaystyle \frac{{\sum }_i{\mathit{\rho }}_i^{\text{Landsat}}(x,y)}{{\sum }_i{\mathit{\rho }}_i^{\text{MODIS}}(x,y)}},\mathrm{\dots },\begin{array}{l}\text{for all pixels}\\ \text{within the}\\ \text{image patch}\end{array}\right\}>\mathrm{2},\end{array}$$

  where (x,y) indicates the pixel location and ${\mathit{\rho }}_i^{\text{Landsat}}(x,y)$ and ${\mathit{\rho }}_i^{\text{MODIS}}(x,y)$ are the surface reflectance of band i for the corresponding Landsat pixel and MODIS pixel, respectively.

• 3.

  Further, a time series outlier detection technique based on the Hampel filter combined with a spatial filter is employed to detect temporal outliers using the vegetation index (VI) of Landsat time series (Claverie et al., 2018).

  $$\begin{array}{}\text{(2)}& \text{VI}\left(t\right)={\displaystyle \frac{{\mathit{\rho }}_{\mathrm{NIR}}^{\text{Landsat}}\left(t\right)}{{\mathit{\rho }}_{\mathrm{Red}}^{\text{Landsat}}\left(t\right)}},\end{array}$$

  where t indicates the time phase and ${\mathit{\rho }}_{\mathrm{NIR}}^{\text{Landsat}}\left(t\right)$ and ${\mathit{\rho }}_{\mathrm{Red}}^{\text{Landsat}}\left(t\right)$ are the NIR and red band surface reflectance of Landsat at time phase t.

For each sample VI(t) of the Landsat time series, the Hampel filter computes the median of the VIs in a temporal window (with the center sample excluded).

$$\begin{array}{}\text{(3)}& \begin{array}{rl}& {\text{VI}}_{\mathrm{median}}=\\ & \mathrm{\text{median}}\mathit{\left\{}\text{VI}\right(t-\mathrm{\Delta }t),\mathrm{\dots },\text{VI}(t-\mathrm{1}),\text{VI}(t+\mathrm{1}),\mathrm{\dots },\text{VI}(t+\mathrm{\Delta }t\left)\mathit{\right\}},\end{array}\end{array}$$

where Δt represents the temporal window size.

Then, it estimates the scale of natural variation (SNV) of each sample by deriving the median of the absolute deviations of the VIs in the temporal window from the median.

$$\begin{array}{}\text{(4)}& \begin{array}{rl}& \text{SNV}\\ & =\mathrm{\text{median}}\mathit{\left\{}\right|\text{VI}(t-\mathrm{\Delta }t)-{\text{VI}}_{\mathrm{median}}|,\mathrm{\dots },|\text{VI}(t-\mathrm{1})-{\text{VI}}_{\mathrm{median}}|,\\ & \left|\text{VI}\right(t+\mathrm{1})-{\text{VI}}_{\mathrm{median}}|,\mathrm{\dots },\left|\text{VI}\right(t+\mathrm{\Delta }t)-{\text{VI}}_{\mathrm{median}}|\mathit{\}}\end{array}\end{array}$$

If the center sample, VI(t), differs from the median, VI_median, by more than 5 SNV, it is flagged as an outlier. No filter is applied if there are fewer than 3 samples within a 60 d temporal window. To eliminate isolated outlier pixels that generate a speckle effect, the sample pixel is flagged as an outlier only if the majority of its surrounding pixels are also flagged as outliers, cloud, or cloud shadow.

3.3 Landsat cross-sensor calibration

To ensure the temporal continuity of the generated SDC dataset, a Landsat cross-sensor calibration approach was employed to reduce the data inconsistencies between the input Landsat OLI and TM and ETM+ products. The linear regression models have been widely used to reduce cross-sensor reflectance difference (Chastain et al., 2019; Claverie, 2023; Claverie et al., 2018; Roy et al., 2016a; Shang and Zhu, 2019). It has been found that a single set of linear transformation coefficients is not proper for global-scale applications (Olthof and Fraser, 2024; Shang and Zhu, 2019). Therefore, our approach aims at building multiple transformation models for each MGRS tile and each spectral band separately.

The Landsat 7 ETM+ and Landsat 8 OLI share a sun-synchronous orbit at an altitude of approximately 710 km. Both satellites revisit a specific area on Earth every 16 d, with an 8 d offset from each other. Notably, there exists a discernible across-track overlap between each ETM+ observation and an adjacent OLI observation acquired 1 d apart (Roy et al., 2016a). Our cross-calibration approach uses matched (acquired 1 d apart) ETM+ and OLI observations in these overlapped regions to build linear regression models for each MGRS tile and each spectral band separately, adjusting ETM+ (and TM) observations to match the OLI spectral responses (Roy et al., 2016a).

For each MGRS tile, ETM+ and OLI observations in the overlap regions acquired only 1 d apart during the years 2014–2021 are reprojected onto the modified MGRS grid using nearest-neighbor resampling. Pixels flagged as cloud, cloud shadow, and snow/ice were discarded. Then, the remaining candidate pixels are used to build linear regression models to adjust reflectance differences for each spectral band. If the number of available pixels is insufficient to build the regression model, candidate pixels from adjacent MGRS tiles would be included. This step helps improve stability and mitigate the issue of spatial inconsistency between neighboring MGRS tiles. Table 3 shows an example of calibration coefficients obtained at three MGRS tiles in this study:

$$\begin{array}{}\text{(5)}& {\mathit{\rho }}_i^{\text{OLI}}=a\times {\mathit{\rho }}_i^{\text{ETM+}}+b\text{for each band}i.\end{array}$$

Table 3 Calibration coefficients obtained at three MGRS tiles in this study.

3.4 MODIS harmonization to Landsat bandpass

Harmonizing MODIS to the Landsat bandpass reduces inconsistencies between Landsat and MODIS observations, which has been proven effective in improving the reconstruction accuracy of subsequent gap-filling and spatiotemporal fusion processes (Chen et al., 2023; Gevaert and García-Haro, 2015; Shi et al., 2022). Here, a cross-resolution data harmonization approach is employed for harmonizing MODIS to the Landsat OLI bandpass. This method involves utilizing matched Landsat and MODIS observations to establish multiple linear transformation models for each spectral band and each local image patch.

In contrast to previous methods that construct distinct transformation models for each land cover type (Cao et al., 2020; Shen et al., 2013; Yang et al., 2020), our approach adopts a patch-wise harmonization strategy with an overlapping mechanism to tackle spatial heterogeneity. This strategy avoids the necessity for high-accuracy land cover maps while concurrently ensuring computational efficiency. As illustrated in Fig. 3, Landsat images are paired with MODIS images acquired on corresponding dates, with the exclusion of contaminated Landsat pixels, such as clouds and cloud shadows. Subsequently, for each image patch and each spectral band, a linear transformation model is constructed utilizing candidate pixels within that image patch. The MODIS reflectance is then adjusted using the obtained transformation coefficients to generate more harmonized MODIS data. In areas of patch overlap, the transformation coefficients are averaged for the final adjustments. This approach employs multiple regional transformation models to better account for material-dependent spectral characteristics that vary across regions and uses an overlapping mechanism to enhance spatial consistency between neighboring image patches.

Figure 3 An illustration of patch-wise harmonization with overlapping.

3.5 Unified gap filling and spatiotemporal fusion

Existing spatiotemporal fusion algorithms generally require cloud-free seamless Landsat images as input (Gao et al., 2006; Shi et al., 2022; Zhu et al., 2010, 2016), which may harm their data efficiency and performances, especially in cloudy areas where there are few cloud-free Landsat images available. Therefore, previous studies (Chen et al., 2018; Liu et al., 2021) applied gap-filling algorithms to partly contaminated Landsat images first and then used these gap-filled images for the subsequent fusion process. Different from these approaches, we propose the Unified, ROBust, OpTimization-based spatiotemporal reconstruction model (uROBOT), which can tackle both gap-filling and spatiotemporal fusion problems in a unified manner.

As shown in Fig. 4, the input data for uROBOT consist of a matched time series of Landsat–MODIS image patches, D_F and D_C, acquired at $\mathit{\{}{T}_{\mathrm{1}},\mathrm{\dots },T_n\mathit{\}}$; a MODIS image patch, C_p, acquired at the prediction phase, T_p; and a partially contaminated Landsat image, F_p, acquired at T_p (case 1 in Fig. 4). The output is the predicted reflectance value for the unobserved segments of the Landsat image, $F_{\mathrm{p}}^{-}$, at T_p. Should F_p be entirely contaminated/unobserved, the scenario simplifies to a conventional spatiotemporal fusion problem (case 2 in Fig. 4).

Figure 4 Input–output settings of uROBOT. The superscripts + and − indicate the corresponding observed and unobserved (or contaminated) segments of Landsat image at the prediction phase. The uROBOT can tackle both gap-filling (case 1) and spatiotemporal fusion (case 2) problems in a unified manner.

3.5.1 Preliminaries of the uROBOT model

Similar to the previous spatiotemporal fusion model (Chen et al., 2023), the basic assumption of uROBOT is that the MODIS image, C_p, can be accurately approximated by a linear combination of other similar MODIS images in the input time series data given by

$$\begin{array}{}\text{(6)}& C_{\mathrm{p}}\approx D_c\mathit{\alpha },\end{array}$$

where C_p is the MODIS image patch acquired at the prediction phase, T_p; D_C is a matrix that stacks the input MODIS time series image patches acquired at $\mathit{\{}{T}_{\mathrm{1}},\mathrm{\dots },T_n\mathit{\}}$; and α is a sparse vector that selects out MODIS image patches similar to C_p and combines them to approximate C_p.

Then, it is assumed that the representation coefficient α can be transferred to corresponding Landsat images, and we obtain the following estimation for F_p:

$$\begin{array}{}\text{(7)}& {\widehat{F}}_{\mathrm{p}}=D_F\mathit{\alpha },\end{array}$$

where ${\widehat{F}}_{\mathrm{p}}$ is the estimated Landsat image patch at T_p and D_F is a matrix that stacks the input Landsat time series image patches acquired at $\mathit{\{}{T}_{\mathrm{1}},\mathrm{\dots },T_n\mathit{\}}$.

The representation coefficients α can be obtained by solving an optimization problem with two extra regularization terms. Detailed explanations regarding these two regularization terms will be provided subsequently.

$$\begin{array}{}\text{(8)}& \begin{array}{rl}\underset{\mathit{\alpha }}{min}& |C_{\mathrm{p}}-D_C\mathit{\alpha }{|}_{\mathrm{2}}^{\mathrm{2}}+\mathit{\lambda }|\mathit{\alpha }{|}_{\mathrm{1}}+\mathit{\beta }|F_{\mathrm{interp}}-D_F\mathit{\alpha }{|}_{\mathrm{2}}^{\mathrm{2}}\\ & +\mathit{\mu }|F_{\mathrm{p}}^{+}-D_F^{+}\mathit{\alpha }{|}_{\mathrm{2}}^{\mathrm{2}},\end{array}\end{array}$$

where λ, β, and μ are scalars; α is a vector that flattens the representation coefficients; C_p is a vector that flattens the MODIS image patch at T_p; D_C is a matrix that stacks the MODIS time series image patches acquired at $\mathit{\{}{T}_{\mathrm{1}},\mathrm{\dots },T_n\mathit{\}}$; D_F is a matrix that stacks the Landsat time series image patches acquired at $\mathit{\{}{T}_{\mathrm{1}},\mathrm{\dots },T_n\mathit{\}}$; F_interp is a vector that flattens an interpolated Landsat image patch at T_p (further elucidation on this matter is subsequently presented); $F_{\mathrm{p}}^{+}$ is a vector that flattens the observed part of F_p; and $D_F^{+}$ is a matrix that stacks the corresponding observed parts of D_F.

Beyond the previous spatiotemporal fusion model (Chen et al., 2023), the uROBOT model accepts the observed part $F_{\mathrm{p}}^{+}$ (if there is one as in case 1 in Fig. 4) as extra input. Hence, the regularization term $\mathit{\mu }\left(\right|F_{\mathrm{p}}^{+}-D_F^{+}\mathit{\alpha }\left|\right)$ allows uROBOT to exploits the observed segments $F_{\mathrm{p}}^{+}$ to better reconstruct the target image, F_p. This feature enables the uROBOT model to handle both gap-filling and spatiotemporal fusion problems in a unified manner. Additionally, uROBOT utilizes the interpolated Landsat image, F_interp, in the temporal continuity penalty term, which further improves the performance of uROBOT. F_interp is an interpolated patch of Landsat imagery at time T_p, generated by combining Landsat observations acquired at nearby dates, with weights determined by the corresponding MODIS data. Therefore, the regularization term $\mathit{\beta }\left(\right|F_{\mathrm{interp}}-D_F\mathit{\alpha }{|}_{\mathrm{2}}^{\mathrm{2}})$ serves to enhance the temporal continuity of the final predicted Landsat image, ${\widehat{F}}_{\mathrm{p}}=D_F\mathit{\alpha }$.

All the constraint terms in Eq. (8) contribute to addressing gradual and step changes. To handle extreme conditions such as ephemeral land cover changes, the uROBOT model also distributes the approximation residuals into the prediction (Chen et al., 2023), and the final prediction is formulated as follows:

$$\begin{array}{}\text{(9)}& {\widehat{F}}_{\mathrm{p}}=D_F\mathit{\alpha }+(C_{\mathrm{p}}-D_C\mathit{\alpha }).\end{array}$$

In regions with frequent cloud cover, the scarcity of cloud-free observations can pose a challenge. To address this, the temporal continuity constraint $\mathit{\beta }|F_{\mathrm{interp}}-D_F\mathit{\alpha }{|}_{\mathrm{2}}^{\mathrm{2}}$ and the residual distribution in Eq. (9) ensure that our estimations are consistent with C_p and are at least as accurate as the interpolated results, F_interp.

3.5.2 Implementation of uROBOT for SDC reconstruction

As shown in Fig. 4, the uROBOT model reconstructs seamless Landsat images for each prediction phase separately. At each prediction phase, T_p, the uROBOT model takes three main steps to reconstruct the corresponding Landsat image, F_p:

• 1.

  Firstly, the uROBOT model constructs an interpolated image, F_interp, using the weighted combination (Zhu et al., 2010) of clear-sky Landsat pixels acquired nearest to the prediction phase, T_p.

  $$\begin{array}{}\text{(10)}& \begin{array}{rl}{F}_{\mathrm{interp}}(x,y)& =w_{\mathrm{1}}\times F_{\mathrm{1}}(x,y)+w_{\mathrm{2}}\times F_{\mathrm{2}}(x,y)\\ & \text{for each pixel location}(x,y),\end{array}\end{array}$$

  where (x,y) indicates a given pixel location, F₁(x,y) is the cloud-free Landsat pixel acquired nearest to and before T_p, F₂(x,y) is the cloud-free Landsat pixel acquired nearest to and after T_p, and the weights w₁ and w₂ are obtained using corresponding MODIS pixels.

  $$\begin{array}{c}{\displaystyle }{w}_{\mathrm{1}}={\displaystyle \frac{\left(C_{\mathrm{2}}\right(x,y)-C_{\mathrm{p}}(x,y){)}^{\mathrm{2}}}{\left(C_{\mathrm{1}}\right(x,y)-C_{\mathrm{p}}(x,y){)}^{\mathrm{2}}+(C_{\mathrm{2}}(x,y)-C_{\mathrm{p}}(x,y){)}^{\mathrm{2}}}}\\ \text{(11)}& {\displaystyle }\text{and}{\displaystyle }{w}_{\mathrm{2}}={\displaystyle \frac{\left(C_{\mathrm{1}}\right(x,y)-C_{\mathrm{p}}(x,y){)}^{\mathrm{2}}}{\left(C_{\mathrm{1}}\right(x,y)-C_{\mathrm{p}}(x,y){)}^{\mathrm{2}}+(C_{\mathrm{2}}(x,y)-C_{\mathrm{p}}(x,y){)}^{\mathrm{2}}}},\end{array}$$

  where C_p(x,y) is the corresponding MODIS pixel acquired at T_p, C₁(x,y) is the corresponding cloud-free MODIS pixel acquired on the same date as F₁(x,y), and C₂(x,y) is the corresponding cloud-free MODIS pixel acquired on the same date as F₂(x,y).

• 2.

  Secondly, the uROBOT model utilizes all input time series data and the spatial information of $F_{\mathrm{p}}^{+}$ to do similar image matching and approximation (Chen et al., 2023), by solving the optimization problem in Eq. (8) to obtain the representation coefficient α.

• 3.

  Then, the last step is to reconstruct the target $F_{\mathrm{p}}^{-}$ using the obtained coefficients α and all input data, as

  $$\begin{array}{}\text{(12)}& {\widehat{F}}_{\mathrm{p}}^{-}=D_F^{-}\mathit{\alpha }+(C_{\mathrm{p}}^{-}-D_C^{-}\mathit{\alpha }),\end{array}$$

  where ${\widehat{F}}_{\mathrm{p}}^{-}$ is an estimation of the unobserved/contaminated segments of F_p and $D_F^{-}$, $C_{\mathrm{p}}^{-}$, and $D_C^{-}$ are the corresponding masked parts of D_F, C_p, and D_C.

3.6 Quantitative assessment using the leave-one-out approach

The leave-one-out assessment approach has been widely used to evaluate the accuracy of reconstructed surface reflectance values in previous gap-filling and spatiotemporal fusion studies (Chen et al., 2011; Gao et al., 2006; Zhu et al., 2016, 2010). This approach initially excludes a certain Landsat image from the input data, subsequently using the remaining input data to reconstruct the excluded image and thereafter evaluating discrepancies between the originals and reconstructions using standard metrics, such as the correlation coefficient (CC), root mean square error (RMSE), mean absolute error (MAE), rMAE (MAE normalized by true surface reflectance values), and Roberts edge (Edge) spatial features (Zhu et al., 2022). We calculated the normalized difference in the Edge spatial features between a reconstructed image and actual Landsat image. The normalized metric values range from −1 to 1, indicating the under- or over-estimate of spatial details. The average normalized metric value of pixels with an Edge value higher than the 90th percentile in the actual Landsat image was used to represent the spatial accuracy of the reconstructed image (Zhu et al., 2022).

A total of 425 test sites that are distributed across the globe are selected randomly, each of which covers a 6 km×6 km area. These test sites were grouped by their dominant land cover types using the FROM_GLC land cover map (Gong et al., 2013; C. Li et al., 2017). Figure 5 shows the spatial distribution and corresponding dominant land cover types of the 425 global test sites. For each Landsat image at each test site, we exclude it from the input data and apply SDC reconstruction using the remaining data, separately. Then, the reconstruction accuracy is evaluated by comparing the reconstructed images with the original Landsat images. The accuracy assessment is conducted in four years (2001, 2004, 2012, and 2021) representative of different input data conditions.

Figure 5 Spatial distribution and corresponding dominant land cover types of the 425 global test sites.

3.7 Quantitative assessment by cross-comparing with HLS products

Another assessment approach is to cross-compare the gap-filled SDC dataset with actual observations from other sensors. NASA's HLS products provide dense 30 m observation data by harmonizing Landsat OLI (since 2013) and Sentinel-2 MSI (since 2015) products, making themselves good reference data to evaluate the SDC product for the period 2016–2022. We selected 22 MGRS tiles representative of different land cover types as test sites (Chen et al., 2023) and evaluated the agreement between SDC and HLS products in the year 2021. Figure 6 shows the spatial distribution of the 22 MGRS tiles in the cross-comparison. Each MGRS tile covers a 109.8 km×109.8 km area. Since the L30 product is derived from Landsat OLI data, we employed the leave-one-out validation strategy (the same as in the last section) for the cross-comparison with L30 data. The least-squares regression method was used to reduce spectral bandpass differences between SDC and HLS data for each spectral band and each test site.

Figure 6 Spatial distribution of the 22 MGRS tiles involved in the cross-comparison with HLS.

4 Results

4.1 Global 30 m daily seamless data cube (SDC) of land surface reflectance

Based on the developed framework, a global 30 m 23-year (2000–2022) daily surface reflectance SDC dataset was generated by combining multi-sensor observations from Landsat TM, ETM+, and OLI and MODIS Terra products. The generated SDC dataset is tiled into the UTM-based grid as described in Sect. 3.1. This gridding system includes 18 466 tiles (each of which includes 3661×3661 Landsat pixels), covering most of global land surface except Antarctica as shown in Fig. 7.

Figure 7 Distribution of the 18 466 tiles used for SDC generation.

Figures 8–11 depict four examples of SDC time series in comparison with the HLS data (not used in the generation of SDC), with a specific focus on land cover changes. Figure 8 illustrates a crop harvest event that took place in April 2021 in Egypt. Remarkably, the SDC time series demonstrates a more adept representation of the various crop harvest stages compared to the L30 time series. This enhanced performance is attributed to the incorporation of more frequent temporal information from MODIS. Figure 9 presents the second case in the Poyang Lake wetland, a region prone to frequent cloud cover. Although there are only limited cloud-free Landsat observations, the temporal phases of land–water transition are effectively captured in the SDC time series. The third case typified a tundra region in Canada. As shown in Fig. 10, the SDC time series accurately reflects the snow season and vegetation growth in the tundra ecosystem. The fourth case in Fig. 11 illustrates a flood event that occurred in India. This region experienced flooding in October 2013, a period when Sentinel-2 data were unavailable. During December, the water gradually receded as depicted in the figure. The results demonstrate the proficiency of SDC data in monitoring rapid land cover changes. Further, the SDC dataset exhibits robust consistency in both spatial and temporal dimensions, with an extended temporal coverage dating back to 2000.

Figure 8 False-color composites of SDC and HLS data and time series of red SR (red) and NIR SR (blue) of the central pixel located at 30.4851° N, 31.9383° E (Egypt; tile 36RUU) from 15 February to 15 July 2021.

Figure 9 False-color composites of SDC and HLS data and time series of red SR (red) and NIR SR (blue) of the central pixel located at 29.0965° N, 116.1107° W (China; tile 50RMT) from 15 February to 15 July 2021.

Figure 10 False-color composites of SDC and HLS data and time series of red SR (red) and NIR SR (blue) of the central pixel located at 56.2778° N, 110.9034° W (Canada; tile 12VWH) from 15 February to 15 June 2021.

Figure 11 False-color composites of SDC and HLS data and time series of red SR (red) and NIR SR (blue) of the central pixel located at 23.2923° N, 68.7947° E (India; tile 42QVL) from 15 November 2013 to 15 March 2014.

4.2 Accuracy assessment results using the leave-one-out approach

Table 4 displays the results of quantitative assessments for the 425 global test sites across four typical years, utilizing different input data settings. Notably, the SDC attains its peak reconstruction accuracy in 2021, when incorporating Landsat ETM+, OLI, and MODIS data. The reconstruction accuracy of SDC is comparatively diminished in 2012, when Landsat TM and OLI observations are not available. Table 5 illustrates the SDC reconstruction accuracy across test sites characterized by different land cover types. The reconstruction accuracy of SDC in waterbodies and tundra areas is observed to manifest relatively higher error levels. Figure 12 displays scatterplots depicting the predicted SDC and actual Landsat surface reflectance values in the leave-one-out assessment. The majority of data points closely align with the 1:1 line, indicating a robust consistency between the predicted and actual reflectance values. Overall, these results substantiate that reconstructed SDC surface reflectance values achieve a high level of accuracy on a global scale.

Table 4 Accuracy of SDC reconstruction at the 425 global test sites (^∗: scan-line corrector failure).

Table 5 Accuracy of SDC reconstruction at the 425 global test sites for different land cover types.

Figure 12 Scatterplots of the actual and predicted values for the test sites of different land cover types.

4.3 Accuracy assessment results by cross-comparing with HLS products

Table 6 presents the quantitative cross-comparison results between SDC and HLS L30 and S30 products at the 22 MGRS tiles. The results indicate a higher level of agreement between SDC and HLS L30 products compared to that between SDC and HLS S30 products. This discrepancy may be attributed to the spectral differences between Landsat OLI and Sentinel-2 MSI instruments. Notably, since the metric values listed in Table 6 are calculated for test sites at a different spatial scale (109.8 km×109.8 km here and 6 km×6 km in Sect. 4.2), the listed root mean square deviation (RMSD) and CC values are not directly comparable to those in Tables 4 and 5. As described in Sect. 3.1, the gridding system of SDC has a 15 m offset compared to the Sentinel-2 gridding system used by HLS products. This 15 m offset may also cause systematic deviations in the cross-comparison.

Table 6 Cross-comparison results of between SDC and HLS L30 and S30 products at the 22 MGRS tiles.

Figure 13 presents scatterplots depicting the reconstructed SDC data in comparison to HLS data. The blue band in the plots indicates higher deviations, known to be sensitive to atmospheric conditions. Most data points lie near the 1:1 line, indicating a high degree of agreement between SDC and HLS products.

Figure 13 Scatterplots of SDC and HLS reflectance values in the cross-comparison.

4.4 The effectiveness of Landsat cross-sensor calibration

The evaluation of our Landsat cross-sensor calibration method involved a comparison of surface reflectance differences before and after the bandpass adjustment. The surface reflectance differences are quantified by the RMSDs between ETM+ and OLI observations (temporally interpolated to align with the acquisition dates of ETM+). This comparison is conducted for six spectral bands across six MGRS tiles (10SEH, 17RKQ, 32UNA, 36MWS, 49SGV, and 55HCC). Figure 14 presents the surface reflectance differences before and after the cross-calibration process for six spectral bands at the six test sites. Additionally, Fig. 15 illustrates two examples of calibrated blue and NIR surface reflectance time series. The results demonstrate the effectiveness of our method in reducing the surface reflectance differences between ETM+ and OLI observations, aligning them more closely using obtained local-scale linear transformation models.

Figure 14 The mean RMSDs between ETM+ and OLI observations before and after the cross-calibration for six spectral bands at six MGRS tiles.

Figure 15 Time series of original and calibrated surface reflectance: (a) located at 23.2923° N, 68.7947° E; (b) located at 23.2923° N, 68.7947° E. The blue lines are the estimated curves of the OLI observations based on the Fourier approach (Dash et al., 2010; Shang and Zhu, 2019).

4.5 The effects of remaining Landsat sensor differences

The effects of remaining Landsat sensor differences on the SDC reconstruction accuracy are investigated in this section. We evaluated the reconstruction accuracy with different input data settings in three typical years (2001, 2004, and 2021) using the 425 global test sites from Sect. 3.6. The accuracy assessment results presented in Table 7 suggest that incorporating more data from different sensors helps to improve the overall reconstruction accuracy. However, a more in-depth analysis reveals that while incorporating ETM+^∗ (^∗ denoting a scan-line corrector failure) images in the input data improved the accuracy of TM image reconstructions in 2004, it simultaneously diminished the accuracy of OLI image reconstructions in 2021. In a similar vein, the inclusion of OLI somewhat reduced the accuracy of ETM+^∗ image reconstructions in 2021 as evidenced in Table 8 and Fig. 16. Despite this, these fluctuations in accuracy were relatively minor.

Table 7 Accuracy of SDC reconstruction (of all images) with different input data settings. Metrics were averaged over the six bands, and the best results are marked in bold. (^∗: scan-line corrector failure).

Table 8 Accuracy of SDC reconstruction (of images from specific sensors) with different input data settings. Metrics were averaged over the six bands, and the best results are marked in bold.

Figure 16 Comparison of the SDC reconstruction accuracy with different input data settings. The symbols + (increase) and − (decrease) indicate the change in accuracy with the additional input data.

4.6 The effectiveness of Landsat cloud masking

To mitigate cloud mask omission errors in Fmask results, SDC generation incorporated an enhanced cloud masking approach involving a brightness-threshold filter and a time series outlier filter. Figure 17 presents two cases of cloud masking results by Fmask and the enhanced cloud masking method. These cases reveal the presence of cloud or heavy aerosol pixels that remain undetected by the Fmask algorithm. Meanwhile, the enhanced cloud masking method adeptly identified and filtered out the majority of these previously undetected cloud pixels, demonstrating its effectiveness in ensuring the quality of the input Landsat data for SDC generation.

Figure 17 Comparison of cloud masks derived by Fmask and the enhanced cloud masking approach. The yellow shading in the cloud masks signifies cloud-free pixels. It can be seen that thin clouds undetected by Fmask were mostly screened out by the enhanced cloud masking.

4.7 Temporal continuity of reconstructed SDC time series

The SDC product provides an extensive historical dataset of multi-spectral and medium-resolution data covering the period from 2000 to 2022. As indicated in Table 2, this prolonged temporal coverage encapsulates multiple distinct phases, each characterized by varied input data configurations. The temporal continuity and consistency of SDC data emerge as pivotal factors for its efficacy in monitoring long-term land dynamics. Figure 18 presents two distinct scenarios of SDC time series extending from 2000 to 2022. The first scenario pertains to a forest situated in a mountainous region, characterized by strong phenological changes and rugged terrain features. The second scenario depicts a desert location exhibiting minimal temporal variations and near-Lambertian surfaces. It can be seen that there are no discernible discontinuities between the different phases in the SDC time series.

Figure 18 Two cases of SDC and Landsat time series: (a) located at 35.4549° N, 113.4381° E and (b) located at 39.0850° N, 81.2460° E.

5 Discussion

5.1 Global-scale land cover classification utilizing the SDC: a comparative analysis with Landsat composite and interpolated datasets

Image compositing is a conventional approach employed to mitigate data gaps in optical remote sensing, which selects the highest-quality observations within a defined time interval based on specific criteria to create seamless clear images (Qiu et al., 2023). Contrary to best-pixel composite methods, there are also some interpolation methods to generate synthetic images based on harmonic time series fits (Zhou et al., 2022; Zhu et al., 2015b). Figure 19 displays the comparison of SDC images with Landsat composite and interpolated images. The Landsat composite images are generated by the NLCD (Jin et al., 2023) and MAX-RNB (Qiu et al., 2023) algorithms, and the Landsat interpolated images are generated using the HAPO algorithm (Zhou et al., 2022). Notwithstanding the good visual quality apparent in these Landsat composite and interpolated images, it is evident that certain temporal change information was lost during the process of image compositing and interpolation.

Figure 19 False-color composites of Landsat, Landsat composite, and SDC images at 35.3104° N, 97.7974° W.

Moreover, we conducted a comparative analysis to evaluate the effectiveness of using Landsat composite and interpolated images versus using SDC data for land cover classification. To perform this assessment, we utilized the validation sample set provided in a previous study (C. Li et al., 2017), which consists a total of 32 946 data points distributed globally. To facilitate a clear and equitable comparison, we employed solely the six-band surface reflectance time series from Landsat composite interpolated images and SDC data as input features for the classification. Subsequently, we conducted a fivefold cross-validation procedure independently for each input data configuration using the same LightGBM classifier with default parameter settings.

Table 9 presents the results of overall classification accuracy for various input data configurations. It is discernible that there exist minor discrepancies in performance between the NLCD and MAX-RNB results. And the utilization of seasonal composite images led to a significantly higher level of classification accuracy compared to annual composite images. Remarkably, the highest classification accuracy was attained using SDC time series as input features, outperforming other input settings by a wide margin (2.4 %–11.3 %). Figure 20 shows two examples of land cover classification results when using SDC can correctly identify land cover types, while other data cannot. The SDC time series preserves the temporal information of the original Landsat data and remains stable when Landsat observations are sparse.

Table 9 Overall accuracy (OA) of land cover classification results using NLCD composite images, MAX-RNB composite images, HAPO interpolated images, and SDC time series as input features.

Figure 20 Examples of land cover classification results using different input data where (a) the ground truth is cropland and (b) the ground truth is grass.

5.2 The influence of Landsat 7 ETM+ scan-line corrector failure

The scan-line corrector (SLC) of the ETM+ sensor on board Landsat 7 failed permanently in May 2003, resulting in about 22 % of the pixels per scene not being scanned (causing wedge-shaped gaps) since that time (Chen et al., 2011). These wedge-shaped gaps are evident in the Landsat ETM+ images as displayed in Figs. 17–19. However, this issue poses no obstacle to SDC generation as the uROBOT model can exploit time series input data to fill these gaps and fuse them with MODIS in a unified manner. In the year 2012, when there are only Landsat ETM+ SLC-off images available, the uROBOT model can still utilize these ETM+ image patches as input and achieve a relatively high reconstruction accuracy as indicated in Table 4. Even in later years when Landsat OLI images become available, the incorporation of Landsat ETM+ images as input can still help improve the overall accuracy of the SDC generation as indicated in Table 7.

5.3 Limitations of SDC products

The SDC is not equivalent to actual daily 30 m Earth observation data. It is an estimation based on Landsat and MODIS time series observations. Reconstructing missing Landsat data is an under-determined problem, meaning there can be infinitely many possible solutions (Shen et al., 2015). Using 500 m MODIS images as “guidance” and applying the constraints presented in Eq. (8), we can narrow down the solution space and make more accurate estimations. However, achieving 100 % accuracy is not feasible since the information provided in the input data is usually incomplete. Additionally, the effective spatial resolution of MODIS observations changes significantly due to the variations in view angles (Pahlevan et al., 2017). Even after BRDF normalization and temporal smoothing, these effects cannot be perfectly mitigated. The effective temporal resolution of SDC depends on the quality of the input Landsat and MODIS data, which can vary in space and time.

The LEDAPS atmospheric correction algorithm for Landsat TM and ETM+ does not perform as well as the LaSRC for Landsat OLI (Vermote et al., 2018). Additionally, the Fmask algorithm for Landsat 5 and 7 is less effective compared to its performance for Landsat 8 and 9 due to the absence of the cirrus band (Zhu et al., 2015a). Our enhanced cloud masking approach identifies most of the previously undetected clouds. Nonetheless, there are still certain thin aerosols and cloud shadows in Landsat imagery that can introduce temporal noise and spatial artifacts into the SDC dataset. Implementing a more aggressive cloud masking strategy could minimize these impacts, yet it could also significantly reduce the number of available observations. Therefore, the development of a more accurate and robust algorithm for cloud and cloud shadow detection is essential for future improvements.

As shown in Table 8 and Fig. 16, remaining Landsat cross-sensor inconsistencies may slightly reduce the reconstruction accuracy. And the sensor differences between Landsat and MODIS sensors could also introduce errors. Future improvements may require a more effective method for cross-sensor calibration and data harmonization.

The influence of geographic registration errors was not considered in this study since Landsat Collection 2 products have significantly improved the absolute geolocation accuracy of Landsat data (Crawford et al., 2023). However, the co-registration accuracy between Landsat imagery and MODIS data could also influence the reconstruction accuracy of SDC dataset. Future improvements could introduce co-registration processes to address this concern.

The SDC dataset could also benefit from integrating more data sources. For example, the European Space Agency (ESA) launched two satellites of the Sentinel-2 mission (S-2A and S-2B) in 2015 and 2017, respectively. The Multi-Spectral Instrument (MSI) on board both satellites acquires multi-spectral data at a spatial resolution of 10 to 60 m depending on the wavelength with a 5 d revisit period at the Equator. The incorporation of Sentinel-2 could facilitate the generation of higher-quality SDC datasets with finer spatial resolutions.

6 Code and data availability

The SDC dataset is available at https://doi.org/10.12436/SDC30.26.20240506 (Chen et al., 2024) or on the project website (http://sdc.iearth.cloud/, last access: 16 November 2024), where a web-based interface is provided for all researchers to be able to freely access the SDC dataset. Notably, the SDC dataset is dynamically generated upon request to optimize data storage efficiency. Furthermore, leveraging the generated SDC dataset, we have produced global 23-year (2000–2022) annual land cover maps using the FROM_GLC classification system, which are readily accessible on the website. All code utilized in our analyses and the accompanying experimental data are available upon reasonable request.

7 Conclusion

In this study, a global-coverage 30 m resolution 23-year (2000–2022) daily seamless data cube (SDC) of land surface reflectance was developed based on the fusion of multi-sensor observations from the Landsat 5, 7, 8, and 9 and MODIS Terra constellations. SDC generation relies on a novel processing framework, which comprises a set of processing stages: gridding and reprojection, Landsat cloud masking, Landsat cross-sensor calibration, MODIS harmonization to Landsat bandpass, and a unified gap-filling and spatiotemporal fusion stage. The quality of the generated SDC dataset was evaluated using a leave-one-out approach and a cross-comparison with NASA's Harmonized Landsat and Sentinel-2 products. The leave-one-out validation at 425 global test sites assessed the agreement between the SDC with actual Landsat surface reflectance values (not used as input), revealing an overall mean absolute error (MAE) of 0.014 (the valid range of surface reflectance values is 0–1). The cross-comparison of SDC with HLS products at 22 Military Grid Reference System (MGRS) tiles revealed an overall mean absolute deviation (MAD) of 0.017 with L30 (Landsat 8-based 30 m HLS product) and a MAD of 0.021 with S30 (Sentinel-2-based 30 m HLS product).

The SDC has several advantages compared with other existing Landsat-based surface reflectance datasets: (i) it exhibits a higher observation frequency and enhanced capabilities for monitoring land cover changes, (ii) it is consistent in both spatial and temporal dimensions without missing values, (iii) it has a global coverage and a prolonged 23-year duration from 2000 to 2022, and (iv) validation results revealed a high level of accuracy in SDC reconstruction. Moreover, the experiment employing the SDC for global-scale land cover classification underscores its advantages compared with Landsat composite/interpolated datasets, achieving a sizable improvement in overall accuracy (2.4 %–11.3 %). The SDC will be a competitive analysis-ready surface reflectance dataset in many other studies related to global environmental monitoring.