Overview of the SDC processing
Figure 2 illustrates the overview of the SDC processing chain, comprising five key processing steps.
1. Gridding and reprojection
The UTM-based Military Grid Reference System (MGRS) was chosen as the projection system for the SDC product, which 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. Since the original Landsat coordinate system exhibits a half-pixel (15 meters) offset relative to the Sentinel-2 grid, we expanded and shifted the original MGRS grid by 15 meters in each direction to align with the Landsat coordinate system. The SDC product is gridded into this modified MGRS system, with a tile size of 109.83×109.83 km (3661×3661 Landsat pixels). Although it has been found that the overlaps among neighbouring 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 into the modified MGRS grid, using nearest neighbour resampling for Landsat and bilinear resampling for MODIS. The 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.
2. Landsat cloud masking
Cloud and cloud shadow masks are essential for removing contaminated Landsat pixels in SDC generation. We used the Fmask (Zhu and Woodcock, 2012) detection results as primary cloud and cloud shadow indicators. There are still few clouds and heavy aerosols in Landsat images 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 meters (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 predefined threshold, the entire image patch is then flagged as cloudy. The threshold was set at 2 in this study, as this value effectively eliminates most residual clouds without being overly aggressive.
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/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 are 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 days, with an 8-day offset from each other. Notably, there exists a discernible across-track overlap between each ETM+ observation and an adjacent OLI observation acquired one day apart (Roy et al., 2016a). Our cross-calibration approach uses matched (acquired one day 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 one day apart during 2014-2021 are reprojected into the modified MGRS grid using nearest neighbour 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. Table 3 shows an example of calibration coefficients obtained at three MGRS tiles in this study.
4. MODIS harmonization to Landsat bandpass
Harmonizing MODIS to Landsat bandpass reduces inconsistencies between Landsat and MODIS observations, which has been proved effective to improve 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 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 Figure 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.
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., 2016, 2010), 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 use these gap-filled images for 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 Figure 4, the input data for uROBOT consist of a matched time series of Landsat-MODIS image patches and acquired at , a MODIS image patch acquired at the prediction phase , and a partially contaminated Landsat image acquired at (case 1 in Figure 4). The output is the predicted reflectance values for the unobserved segments of the Landsat image at Should be entirely contaminated/unobserved, the scenario simplifies to a conventional spatiotemporal fusion problem (case 2 in Figure 4).
Webmaster: Shuang Chen
Email: schen17 at connect.hku.hk
Last Updated: Feb 11, 2024