2.1 MODIS surface reflectance products MOD09GA (version 6.1) and MYD09A1 (version 6.1)

The MOD09GA product is a daily surface reflectance dataset with a 500 m spatial resolution obtained from the morning satellite EOS-terra (Barnes et al., 1998; Vermote et al., 2011). It contains a quality band that indicates whether the pixel is affected by snow, cloud, or cloud shadow. The MOD09GA dataset is the fundamental input data of SDC500, and cloud-contaminated pixels were removed in the pre-processing stage. The MYD09A1 dataset provides 8 d surface reflectance at a 500 m spatial resolution, acquired by the afternoon satellite EOS-aqua. Each pixel of MYD09A1 contains the best observation among the 8 d. Both the MOD09GA and MYD09A1 data were accessed from the Land Processes Distributed Active Archive Center: https://lpdaac.usgs.gov (last access: 15 August 2023).

2.2 MODIS land cover product MCD12Q1 (version 6.1)

MCD12Q1 is a MODIS land cover product that combines the observations from both EOS-terra and EOS-aqua (Sulla-Menashe and Friedl, 2018; Friedl et al., 2002). It adopts the International Geosphere–Biosphere Programme (IGBP) land cover classification system, which divides the global land surface into 16 ecosystems. In this paper, MCD12Q1 is adopted when building the a priori database for BRDF parameters for different land cover types.

2.3 MODIS BRDF parameter product MCD43A1 (version 6.1)

The MODIS BRDF parameter product MCD43A1 is currently the most widely used and acknowledge remote sensing BRDF dataset (Schaaf et al., 2002, 2012), providing three BRDF kernel model parameters, f_iso, f_vol, and f_geo, with a spatial resolution of 500 m, which can be used for BRDF correction. The MCD43A1 product is adopted for generating the parameters for BRDF correction, and it can be easily accessed and processed through the Google Earth Engine (GEE) platform.

2.4 MODIS nadir BRDF-adjusted reflectance product MCD43A4 (version 6.1)

The MODIS nadir BRDF-adjusted reflectance (NBAR) product MCD43A4 is generated to facilitate vegetation monitoring and phenological studies, in which the MOD09 surface reflectance data are corrected to nadir-viewing angles using images in a 16 d sliding widow (Schaaf et al., 2002; Vermote et al., 1997). However, as the inversion of the BRDF model needs the cumulation of multiple clear observations in a 16 d window, the resulting NBAR dataset suffers from missing data and is limited in its applications. The MODIS NBAR product (MCD43A4) is chosen for a comparison with our BRDF-corrected results because it is also derived from MODIS data and also normalized to a standard sun/view geometry.

3.1 Generation of the a priori database

The MODIS BRDF product MCD43A3 is adopted to generate the a priori database for BRDF and phenology as it is currently the most widely used and acknowledged remote sensing BRDF dataset. Directly using the MCD43A3 dataset in a per-pixel and per-date way is not practical due to the presence of noise and missing values. So, we derive robust statistics from the MCD43A3 dataset and utilize them for BRDF correction.

The following formula expresses the kernel-driven RossThick-LiSparse-Reciprocal (RTLSR) model:

where $K_{\mathrm{vol}}({\mathit{\theta }}_{\mathrm{s}},{\mathit{\theta }}_{\mathrm{v}},\mathit{\phi })$ and $K_{\mathrm{geo}}({\mathit{\theta }}_{\mathrm{s}},{\mathit{\theta }}_{\mathrm{v}},\mathit{\phi })$ are two kernel functions representing volumetric scattering and geometric–optical-shadowing effects, respectively; f_iso, f_vol, and f_geo are the three model parameters, which constitute the main content in the MCD43A3 dataset.

Before doing the statistics, we apply a simple transform:

where ω₁ and ω₂ are called BRDF shape factors. They jointly decide the shape of the angular variation of the BRDF model which will be applied to the BRDF correction in the next section.

Land surface classification is also considered in building the a priori database. The MODIS land cover type product MCD12Q1 is adopted for this purpose, and the IGBP land cover classification system within the MCD12Q1 product is adopted, which divided the global land surface into 16 ecosystems.

The cloud computing platform of the Google Earth Engine (GEE) has enabled us to efficiently access and process the large global MCD43A3 and MCD12Q1 datasets. The statistics are derived by computing the median values of f_iso, ω₁, and ω₂ for each MODIS land band, every 4 Julian days, and each land cover type for spatial cells of 100 km × 100 km with a grid space of 50 km. The time scope is from 2000 to 2021, and the spatial extent includes all land surfaces, as well as shallow seas covered by the MCD43A3 dataset. The median statistics are chosen for their resistance to noise, and the coarse spatial resolution ensures an ample sample size of MCD43A3 pixels. However, not all land cover types have valid median statistics in each spatiotemporal cell. For example, a cell in the tropical desert may not have any snow, cropland, or forest samples. In these cases, we filled the missing statistics with the average of neighbouring valid statistics and iterated the filling process until all spatiotemporal cells contained statistics for all land cover types. This process enabled us to generate a global a priori database for BRDF and phenology. The GEE code can be accessed at https://code.earthengine.google.com/363b4d94090048f9e28103ad3efebfdf (last access: 15 October 2023).

3.2 Land-cover-based BRDF correction

To utilize the a priori information on BRDF and phenology, the first step is to select one sequence of time series of f_iso, ω₁, and ω₂ from the a priori database. For each input MODIS pixel, we use the coordinate of the 500 m pixel in the MODIS projection grid to locate the nearest time series of 100 km × 100 km statistics, which corresponds to each of the 17 land cover types. We use the 500 m resolution land cover product MCD12Q1 to acquire the cover type of the pixel. In consideration of ambiguous classifications or mixed pixels, we also select several alternative cover types. For example, if the pixel is indicated as evergreen needleleaf forest in the MCD12Q1 product then evergreen broadleaf forest, deciduous needleleaf forest, mixed forests, and closed shrubland are candidates for alternative cover. The a priori statistics of f_iso, ω₁, and ω₂ for these several possible cover types are used to simulate the candidate sequences of time series according to the actual observation date and sun/view geometry. We then compute the correlation coefficient between the actual MODIS observation and the several simulated candidate sequences. The candidate sequence of f_iso, ω₁, and ω₂, which corresponds to the maximum correlation coefficient, is selected for the following BRDF correction.

The process of BRDF correction is also known as BRDF normalization. As MODIS is a wide-angle sensor, the sun/view angle of a pixel varies from day to day, resulting in fluctuations in the observed time series of the surface reflectance. The objective of BRDF normalization is to simulate a time series that appears as if it was taken from standardized sun/view angles. In this approach, the standardized sun/view angles are defined as nadir viewing and solar illuminating at 10:30 am local time, which is denoted as $({\mathit{\theta }}_{\mathrm{s}\mathrm{0}},{\mathit{\theta }}_{\mathrm{v}\mathrm{0}},{\mathit{\phi }}_{\mathrm{0}})$. Then, the normalization is performed with the following formula:

in which the ω₁ and ω₁ are derived from the a priori information.

3.3 Outlier detection and gap filling

In this study, the gap refers to both the unobserved area of land surface in the daily MODIS acquisition and the area blocked by cloud or heavy aerosols. These gaps need to be filled with a prediction model to build the seamless dataset. Pure mathematical algorithms, such as linear or nonlinear interpolation, can be efficient prediction models when good-quality, valid samples are available near the gap. However, filling the prolonged continuous gaps and gaps with high levels of noise existing in nearby valid samples can be challenging. In these cases, the ecosystem curve-fitting method is the most robust gap-filling algorithm because it utilized the a priori information of phenology. Therefore, we have designed a strategy that combines the mathematical and the phenology curve-fitting algorithms, as illustrated in the second box of Fig. 1.

Step (a) first excludes the invalid observations according to the quality flag in MOD09GA. Then, step (b) provides preliminary fill values with phenology fitting to increase the stability and error resistance of the gap-filling process. As the phenology fitting is more stable and less flexible than the B-spline fitting, we apply it to provide preliminary fill values which will be utilized in the B-spline fitting step to prevent overfitting. And in step (c) of B-spline fitting, the observations with a fitting absolute error larger than 2.5 times the mean are marked as outliers, which usually means undetected cloud contamination or abnormal change of surface state such as ephemeral snow. These outlier values severely destabilize the fitting process and increase error. Therefore, in this study, we remove the outliers from the observation series together with invalid values indicated by the quality flag in the MOD09GA dataset.

The phenology fitting of step (b) can be expressed as the optimization of the following expression:

where F_phe(t) is the a priori phenology curve, which is actually the series of f_iso selected in Sect. 2.2; y₁(t) refers to the preliminary prediction values for the time series; a₀ and a₁ are parameters to be optimized by the least-squares principle, which minimizes the square error between predictions and valid observations in the time series.

The B-spline of step (c) can be expressed as the optimization of the following expression:

where y₂(t) is the prediction by spline; $B_{\mathrm{1}}\left(t\right),\mathrm{\cdots },B_{\mathrm{6}}\left(t\right)$ are six B-spline base functions; and $a_{\mathrm{0}},\mathrm{\cdots },a_{\mathrm{6}}$ are seven parameters to be optimized by the least-squares principle. In this step, the preliminary prediction values are included in the samples together with the valid observations, and the outliers are detected and removed in an iterative way.

In the last step, the combined spline and phenology curve fitting of step (d) is expressed as follows:

In this step, the prediction values by spline from step (c) are included in the samples, and there is no outlier detection iteration. The six B-spline base functions, as well as a demonstration phenology curve, are illustrated in Fig. 2. The demonstration phenology curve is in a red band (MODIS band 1) corresponding to a cropland pixel in the North China Plain.

Figure 2Illustration of the B-spline base functions (SP1,…, SP6) and the surface reflectance phenology curve (PHE) of MODIS band 1.

Download

3.4 Sliding-window smoothing

The gap-filling process produces a time series that includes both actual MODIS observations and predicted filling values. However, this series is not completely smooth due to noise in the observations and inconsistency between the two groups. To ensure that the dataset is smooth, we applied the SG filter algorithm (Chen et al., 2004) to the gap-filled series. This algorithm applies a filter window of 29 d in width, which is centred on each output date and slides through the time series. The algorithm also gives larger weight to the actual MODIS observations than the predicted filling values to ensure that the resulting time series remains as close to the original observations as possible.

3.5 Consideration of snow

The time series processing described above is based on the assumption that surface reflectance changes gradually over time, which is not always the case for natural land surfaces. Abrupt changes, such as snowfall, can result in an abrupt change in surface reflectance, and the melting of snow can also cause rapid changes that cannot be simulated with B-spline and phenology curves. Furthermore, the BRDF feature of snow is completely different from that of snow-free surfaces. To address these issues, a special treatment for snow was developed.

The year is divided into two parts, the snow season and the snow-free season, based on the snow indicator in the quality flag of the MOD09GA product. To ensure a robust segmentation, a median filter with a length of 12 d is applied to the snow indicator series. If there are no valid observations within the range of the median filter, the resulting snow status is copied from the nearest valid neighbour.

The BRDF correction, gap-filling, and smoothing processes are then separately applied to the snow season and snow-free season of the time series. During the processing of the snow season, only the observations marked as snow are used to build a smooth snow reflectance curve and vice versa. Finally, the smooth snow reflectance curve and smooth snow-free reflectance curve are combined to create the final time series of surface reflectance. As a result, the final time series may be discontinuous in the switching date from snow season to snow-free season.

4.1 Visualizing each step in the SDC500 processing

We selected temporal curves of four representative sites to evaluate the effect of each step in the proposed framework, respectively. Figure 5 presents band-2 (NIR) reflectance time series in four representative sites in 2010. The original MODIS sensor observations, the outlier detection and gap-filled curve, the first smoothing results, and the final SDC500 results are shown together to facilitate visual validation.

Site S9 (Fig. 5a) is located in the Amazon Plain in Brazil, representing evergreen forest land cover, as well as a case of scarce clear observation. We can see that, in Site S9, the first gap-filling result by phenology curve (green dash-dot line) is too flat as the a priori phenology curve has a in coarse resolution and lacks details. However, the final SDC500 (red line) effectively reflects the seasonal dynamic while avoiding over-fitting, giving a reasonable result in the tropical zone.

Site S10 (Fig. 5b) is located in the southern part of Africa, representing grassland in the wet season and sparsely vegetated surfaces in the dry season. In Site 10, clear observations are dense, and the reconstruction results are close to the unbiased average status of the surface.

Site S11 (Fig. 5c) is located in the temperate zone of the North China Plain, where winter wheat grows from October to June and where corn grows from July to September. The discrete series (orange cross) of Site S11 show double peaks in summer around DOY 130 and 220 and three valleys around DOY 50, 190, and 300. We can see that the first gap filling (green dash-dot line) cannot approximate the double peak feature of the crop phenology, but the second gap filling (grey line) and the final smoothing (red line) gradually approach the discrete series (cross). It illustrates that SDC500 has stable estimation capabilities near peak and valley values.

Site S12 (Fig. 5d) is located in the sub-Arctic coastal land of Greenland and is featured for the snow–thaw circle and long polar nights in winter. In Site S12, the reconstructed result (red line) is composed of two parts: the snow-free part (140 < DOY < 305) with low reflectance and the snow-covered part with high reflectance. In the DOY range of 280–300, there are some large fluctuations in the raw observations (orange cross); a small part of the observations is marked as snow covered, and the rest is marked as snow free according to the raw quality flag. Our algorithm judges this period to be more snow free according to the median filter so the final reconstructed result (red line) is closer to the lower envelope of raw observations. In the data range before DOY 50 and after DOY 300, the site enters the polar night, and there is no valid raw observation (orange cross). The prediction in the polar night may not be as accurate as the normal date, but the reconstructed result (red line) still gives a reasonable prediction.

Figure 5Comparison of BRDF-adjusted raw MODIS observations, gap filled values with the phenology curve, gap-filled values with the B-spline + raw observations, and the final smoothed series in typical sites: Site S9 evergreen rain forest (a), Site S10 grassland (b), Site S11 cropland (c), and Site S12 arctic (d).

Download

4.2 SDC500 in tropical and subtropical areas

In the tropical and subtropical monsoon climate zones (Figs. 6 and 7), heavy cloud cover significantly affects the satellite observations. Particularly during the monsoon season from May to August, there is continuous and dense cloud coverage, leading to prolonged gaps in valid observations. This poses a considerable challenge for seamless observation reconstruction. Simple linear interpolation performs poorly due to the extensive missing data. But SDC500 can effectively fill these prolonged gaps and reconstructs observation time series that closely reflect the actual conditions.

From a spatial patterns perspective, the reconstructed results of SDC500 demonstrate the preservation of effective observation values in the central cloud-free area in site S1 on DOY 226 compared to the original observations (Fig. 6b). Additionally, the reconstructed data in cloud-covered areas exhibit good spatial consistency with the raw MODIS data in cloud-free areas. This spatial consistency is particularly notable considering that the SDC500 is constructed on a pixel-wise basis. It indicates the reliability of the reconstruction algorithm.

From the temporal series perspective, despite significant cloud interference and numerous data gaps before and after this period, SDC500 successfully reconstructs the missing data and preserves the phenological characteristics of forests, as illustrated in Fig. 6a. In site S2 (Fig. 7a), which is mainly covered by evergreen broadleaf forest, reflectance remains relatively stable throughout the year, and SDC500 effectively eliminates noises and obtains stable phenological curves.

Figure 6Performance of SDC in site S1. (a) The NDVI curves for the central pixel; the orange points indicate the valid MODIS observations, and the red line indicates the SDC500 results. (b) The spatial pattern of reconstructed and raw image blocks of 200×200 pixels centred around the site (red: band 1, green: band 2, and blue: band 3).

Download

Figure 7Performance of SDC in site S2. (a) The NDVI curves for the central pixel; the orange points indicate the valid MODIS observations, and the red line indicates the SDC500 results. (b) The spatial pattern of image blocks of 200×200 pixels centred around the site (red: band 1, green: band 2, and blue: band 3).

Download

4.3 SDC500 in acrid areas

The SDC500 demonstrates excellent noise reduction performance in the acrid area (Figs. 8 and 9). Figure 8b illustrates that the raw MODIS observation on DOY 091 is affected by its large view angle, resulting in no data in the main part of the block, and the remaining upper-left corner is brighter than other raw images. However, these are effectively corrected in the SDC500 reconstructed results, which show consistency with observations on other days. The improved observation consistency on different dates mitigates the inherent data inconsistency arising from sensor component fluctuations, observation/view geometry, and other factors. Notably, in Fig. 9b, the reflectance of water bodies in the images remains consistent throughout the year, without any noticeable darkening or lightening, aligning closely with the physical reality. From the temporal perspective, the NDVI series in S3 (Fig. 8b) peaks at about DOY 250, while the NDVI series in S4 (Fig. 9b) peaks at about DOY 50, which is consistent with the fact that the two sites are in different hemispheres. Additionally, the influence of outlier values has been satisfactorily removed, such as around DOY 200 and 250 in Fig. 8a.

Figure 8Performance of SDC in site S3. (a) The NDVI curves for the central pixel; the orange points indicate the valid MODIS observations, and the red line indicates the SDC500 results. (b) The spatial pattern of image blocks of 200×200 pixels centred around the site (red: band 1, green: band 2, and blue: band 3).

Download

Figure 9Performance of SDC in site S4. (a) The NDVI curves for the central pixel; the orange points indicate the valid MODIS observations, and the red line indicates the SDC500 results. (b) The spatial pattern of image blocks of 200×200 pixels centred around the site (red: band 1, green: band 2, and blue: band 3).

Download

4.4 SDC500 in temperate-climate areas

The primary land cover type of site S5 (Fig. 10) is cropland, which thrives in summer and is harvested in autumn each year. Hence, the typical NDVI phenology of cropland can be used as a reference to assess the performance of SDC500. In this evaluation, we assume that the NDVI phenologies in neighbouring years are similar for most of the pixels in the image block. Therefore, comparing the NDVI patterns over 2 years allows us to gauge the credibility of the reconstructed image. The comparison of SDC500 and raw NDVI images for the same day in 2020 and 2021 is presented in Figs. 11 and 12. The reconstructed NDVI series in Fig. 11 exhibits a transition from low to high and then back to low in the DOY range 110–310 in 2020, which aligns with the temporal curve pattern observed in 2021. Consequently, the NDVI values for the same day in different years are comparable. For example, DOY 110 in 2021 (second line in Fig. 13) is relatively cloud free and can be considered to be reliable in terms of reference data. The reconstructed NDVI image for DOY 110 in 2020 (third line in Fig. 12) agrees well with that of DOY 110 in 2021, indicating the accuracy of the reconstructed result. Similarly, DOY 260 in 2021 (second line in Fig. 13) can serve as reference data to validate the reconstructed result of DOY 260 in 2020 (third line in Fig. 12). Furthermore, certain parts of the region on DOY 160 and 310 in 2020 (second line in Fig. 12) are not affected by clouds and can be used as reference data. The reconstructed results on DOY 160 and 310 in 2021 (3rd line in Fig. 13) effectively preserve the high-quality observations in these regions, exhibiting good spatial continuity with the reconstructed results in cloud-contaminated regions. As the SDC500 is produced on a per-pixel basis, the consistency between clear and cloud-contaminated pixels demonstrates the robustness of the algorithm. Additionally, the results reveal the successful restoration of low NDVI values caused by clouds while preserving the low NDVI values of bare-land roads. This indicates the successful reconstruction of spatiotemporally continuous missing pixels by the proposed algorithm.

Figure 10Performance of SDC in site S5. (a) The NDVI curves for the central pixel; the orange points indicate the valid MODIS observations, and the red line indicates the SDC500 results. (b) The spatial pattern of image blocks of 200×200 pixels centred around the site (red: band 1, green: band 2, and blue: Band 3).

Download

Figure 11NDVI curves of the central pixel in a cropland area (site S5) in the period of 2020 and 2021.

Download

Figure 12Assessment of spatial pattern in a cropland area (site S5). True-colour composite (red: band 1, green: band 4, and blue: band 3) and raw MODIS NDVI and SDC500 NDVI imageries on five dates in 2020.

Download

Figure 13Assessment of spatial pattern in a cropland area (site S5). True-colour composite (red: band 1, green: band 4, and blue: band 3) and raw MODIS NDVI and SDC500 NDVI imageries on 5 dates in 2021.

Download

Figure 14 presents a demonstration of the raw and reconstructed reflectance and NDVI time series in S6, which is in a mid-latitude evergreen broadleaf forest area. Like the other sites, the SDC500 dataset in S6 successfully restored the image series of evergreen broadleaf forest: the annual variation of NDVI is within the range 0.3–0.5; the influence of outliers, such as on DOY 100, 252, and 272 in Fig. 14a, has been satisfactorily eliminated.

Figure 14Performance of SDC in site S6. (a) The NDVI curves for the central pixel; the orange points indicate the valid MODIS observations, and the red line indicates the SDC500 results. (b) The spatial pattern of image blocks of 200×200 pixels centred around the site (red: band 1, green: band 2, and blue: band 3).

Download

4.5 SDC500 in snow-dominant areas

The SDC500 demonstrates excellent capability in dealing with snow and ice cover. Figure 15 illustrates its accurate representation of the melting and snowfall processes while retaining the intricate features of needleleaf forests hidden beneath snow and ice. Even in images where both snow and vegetation are present, the SDC maintains spatial continuity effectively. In the Arctic region depicted in Fig. 16, the SDC500 successfully captures the short melting process in summer (DOY 180–270), showing its ability to document the variations in snow and ice cover over time.

Figure 15Performance of SDC in site S7. (a) The band-2 reflectance curves for the central pixel; the orange points indicate the valid MODIS observation, and the red line indicates the SDC500 results. (b) The spatial pattern of image blocks of 200×200 pixels centred around the site (red: band 1, green: band 2, and blue: band 3).

Download

Figure 16Performance of SDC in site S8. (a) The band-2 reflectance curves for the central pixel; the orange points indicate the valid MODIS observation, and the red line indicates the SDC500 results. (b) The spatial pattern of image blocks of 200×200 pixels centred around the site (red: band 1, green: band 2, and blue: band 3).

Download

5.1 Quantitative assessment of the BRDF correction method

In this study, we perform BRDF correction using land cover information from the MCD12Q1 dataset and kernel model parameters from the MCD43A1 dataset. Thus, the uncertainty in these two datasets may result in uncertainty in the BRDF-corrected result. As there is no established ground-truth dataset for land surface BRDF at the 500 m resolution, we evaluate the uncertainty of our proposed BRDF correction method through two approaches: (1) a visual comparison of the MODIS observations before and after the BRDF correction and (2) a quantitative comparison with the MODIS nadir BRDF-adjusted reflectance dataset MCD43A4.

Figure 17 presents a comparison of the time series of raw MODIS observations and BRDF correction results in four typical sites. In Fig. 17, the raw clear observations are directly extracted from the MOD09GA dataset, corresponding to the various sun/view geometries of the actual satellite overpass; the BRDF-corrected results are standardized to the nadir view angle and sun angle of 10:30 local time. As shown in Fig. 17, the BRDF-corrected results become more consistent with each other after normalization, with reduced noise from biased marginal values and, to some extent, alleviated abnormal fluctuations. We also noticed that there are still remanent fluctuations in the normalized reflectance series; this is partly due to the reduced spatial resolution of BRDF parameters in the a priori database and partly due to other sources of uncertainties such as geometric alignment error and atmosphere correction error.

In addition, we conducted quantitative validation between the proposed BRDF-corrected results and the MODIS nadir BRDF-adjusted reflectance dataset MCD43A4 in the four tiles (h11v09, h20v10, h27v05, and h16v02) corresponding to the typical sites S9, S10, S11, and S12. For each MODIS tile, a regular grid of 12×12 pixels was selected as sample pixels. The raw MODIS observations in 2010, along with their corresponding normalized values, were validated against the MODIS NBAR product MCD43A4. The results, shown in Fig. 18, demonstrate that the raw MODIS reflectance shows a remarkable difference compared to the MCD43A4 product, with an RMSE of 0.275 and a bias of 0.1721; in comparison, our BRDF correction method is in good agreement with the MCD43A4 product, with an RMSE of 0.056 and a bias of −0.0085. As SDC500 and MCD43A4 are derived using different methods but are consistent with each other, the consistency indicates that they both provide accurate results.

Figure 17Example of series of raw MODIS observations, BRDF-corrected observations, and the final smoothed series in four typical sites: site S9 evergreen rain forest (a), site S10 grassland (b), site S11 cropland (c), and site S12 arctic (d).

Download

Figure 18Quantitative evaluation of our proposed BRDF correction method: (a) raw MODIS reflectance vs. MODIS NBAR product and (b) SDC500 dataset vs. MODIS NBAR product.

Download

5.2 Quantitative assessment of the gap-filling method

An essential module in the data-processing algorithm is outlier detection and gap filling. In extreme situations of long-lasting cloud coverage, which commonly happens in tropical and monsoon areas, the gap-filling strategy basically determines the shape of the final output time series. In the past, it is common practice in time series reconstruction that the invalid values should be masked out and filled with linear interpolation of valid values before applying the smoothing filter (e.g. the SG filter) to the time series. In this study, we use a combination of phenology-based and spline-based fitting to fill the gaps. Besides, during the time series fitting, outlier values are detected and masked. As the outlier values are most likely incorrectly flagged cloud or cloud shadow pixels, the exclusion of them can further promote the stability of the reconstruction result.

To demonstrate the performance of our proposed gap-filling algorithm in comparison to the linear gap filling, we select an extreme case in the monsoon season in southern Asia. Figure 19a shows an area of 600×1000 pixels in MODIS tile h25v07, which is located in central India. The monsoon brings heavy cloud and precipitation in June and July every year and leaves almost no clear observations for optical remote sensing satellites in these two continuous months. Figure 19b presents the average cloud coverage of this area for each Julian day from 2000 to 2022, as well as the cloud coverage in the year 2019. We can see that there is a rare case of a clear window during Julian days 196 to 198 in 2019, which can serve as a reference truth value to validate the output of the time series reconstruction algorithm. From Fig. 19a it is found that the clear area is only in the centre part of the image; the surroundings are still contaminated with clouds. Thus, we manually outlined the clear area as a region of interest (ROI) to derive validation statistics.

Figure 19Study area of quantitative validation. (a) Original MODIS observation on 2019 DOY 197 in tile h25v07, line 1–600, column 801–1800. Red: band 1, green: band 2, and b: band 3. The region of interest (ROI) outlined by the red box indicates the clear area to derive validation statistics. (b) Average cloud coverage in the period of 2000–2022 and cloud coverage in 2019.

Download

Figure 20Reconstructed results on 2019 DOY 197 in four cases: (a) SDC500 with all valid input data (case 1), (b) linear interpolation with all valid input data (case 2), (c) SDC500 with reduced input data (case 3), and (d) linear interpolation with reduced input data (case 4).

Download

Table 2Quantitative validation of linear interpolation and SDC500 in different bands.

Download Print Version | Download XLSX

Figure 21Quantitative comparison between actual values and predicted values of linear interpolation (a) and SDC500 (b).

Download

We implemented the algorithm in four cases to reconstruct the image series in 2019. Case 1 represents the operational run of the SDC500 gap-filling and smoothing algorithm, with all the available MODIS observations as input. Case 2 is the same algorithm with an imaginary scenario of reduced input in which MODIS observations in Julian days 195 to 199 are excluded from the input dataset; the purpose of this scenario is to simulate a long gap which occurs in other years. Case 3 is a comparison case in which linear gap filling and SG filter smoothing are applied to the same input data as in case 1. Case 4 is a comparison case in which linear gap filling and SG filter smoothing are applied to the reduced input data, as in case 2. Figure 20 illustrates the reconstructed images of the four cases; the good results of Fig. 20a and b indicate that the SDC500's gap-filling results exhibit improved spatial continuity and greater robustness with reduced input. In terms of performance evaluation metrics, lower values of root-mean-square error (RMSE) and mean absolute error (MAE), as well as a higher correlation coefficient (CC), indicate better results. Table 2 and Fig. 21 provide the quantitative comparison between SDC500 and the result of linear interpolation. The SDC500 achieves am RMSE of 0.0496, a MAE of 0.0430, and a CC of 0.9574, demonstrating superior accuracy compared to the result of linear interpolation, which yields an RMSE of 0.0844, a MAE of 0.0631, and a CC of 0.8174.

5.3 Performance of SDC500 in capturing rapid disturbance

The SDC500 product provides high-resolution observations with accurate characterization of the long-term phenology of ecosystems while also retaining the useful localized feature of rapid disturbances. In this study, we chose to demonstrate the product's capability to detect such disturbances in two wildfire events. The first event is the largest recorded tundra fire that occurred near Anaktuvuk River (Jones et al., 2009), which was started by a lightning strike in July 2007 and spread across a large area of 1039 km² till September. The other event is a serious forest fire that occurred in Australia in July 2019, which was caused by drought and heat weave weather and raged for months. To enhance the extent of the burned area, we computed the normalized burnt ratio (NBR) index, with a lower NBR indicating a higher fraction of burned area. Spatially, the SDC500 was able to capture the continuous spreading pattern of the wildfires, as demonstrated in Figs. 22 and 23. Temporally, the NBR decreased sharply during the fires and returned to normal values after the Australia forest wildfire, as seen in Fig. 23b, while remaining low after the Arctic tundra wildfire, as shown in Fig. 22b. However, as the algorithm of SDC500 is based on temporal smoothing, the reconstructed time series of NBR demonstrate a gradual change, though there should be an abrupt change of NBR in the fire events. Data users need to be aware of the uncertainty caused by this smooth-out effect.

Figure 22The performance of SDC500 during the Arctic tundra wildfire in 2007 summer. (a) The spatial pattern of NBR. (b) The NBR series for a central pixel of burned area in 2006 (before burn) and 2007 (burn) and a non-fired pixel in 2007.

Download

Figure 23The performance of SDC500 during the Australia forest wildfire. (a) The spatial pattern of NBR. (b) A central pixel of burned area in 2018 (before burn) and 2019 (burn) and a non-fired pixel in 2019.

Download

5.4 Quality assessment (QA) and limitations of SDC500

This study proposes an advanced framework for reconstructing the spatiotemporally seamless MODIS surface reflectance dataset. The preliminary evaluations demonstrate the effectiveness of the SDC500 in addressing the challenge of missing key points, filling isolated and continuous gaps, as well as gaps at peaks and valleys. Moreover, the SDC500 successfully reduces noise from biased outlier values and eliminates abnormal fluctuations caused by observation conditions rather than true land surface dynamics.

Despite the satisfactory results achieved using the proposed methodology, uncertainties should be assessed as auxiliary information to support various applications of data users. Pixelwise quality assessment (QA) is the commonly adopted means to express the uncertainties in the dataset, which is also adopted here. Although the ideal way is to express the uncertainties quantitatively, it is not applicable in the current version of SDC500 for two reasons: (1) the data-processing chain contains many sub-algorithms, as well as feedbacks, which make the quantitative uncertainty assessment very difficult, and (2) there is still a lack of quantitative uncertainty indication for the input data of surface reflectance. Thus, the QA flag of SDC500 only gives a relative indicator to the quality of the reconstructed reflectance value. In the QA flag, which is 16 bits unsigned integer, the overall quality is indicated in the lowest 2 bits, i.e. bits 0–1, with 00 representing the best credibility and 11 representing the worst credibility. Figure 24 presents the average percentage for the overall quality bits in each MODIS tile over the year 2020. Pixels over ocean and pixels in polar night are excluded from the percentage statistics. It shows that the tiles with low-quality pixels are mainly distributed in the tropical and subtropical zones or in tiles over small islands, while the main reason for the decreased quality is heavy cloud coverage.

Figure 24The average percentage for the overall quality bits of QA flag in each MODIS tile over the year 2020.

Apart from the QA flag, some of the known issues have not been addressed in the current version of algorithm and may bring extra uncertainties to the final product. Firstly, topographic effect is very important and may bring large uncertainty in the remote sensing estimated surface parameters (Wen et al., 2018). However, the operational topology correction method is mostly empirical rather than quantitative, and topographic research in the resolution of 500 m is even more scarce. Thus, topology is totally ignored in our current version of the algorithm, and data users should take extra caution when using the SDC500 product in mountainous areas. Secondly, the mixed-pixel issue is also not addressed in our current version of the algorithm. Mixed pixels widely exist in 500 m resolution MODIS data. If the pixel is composed of the same or similar land covers, the influence is small. But if the pixel is located in a distinct boundary, such as between land and water or between rural and urban areas, the uncertainty of the reconstructed reflectance will increase, but this is not indicated in the QA flag.

Some other issues have only been partly addressed and need to be improved upon in future research. For example, the BRDF correction relies heavily on a priori parameters derived from the MCD43A1 dataset and the empirical RTLSR model. However, the RTLSR model, which was originally designed for vegetated land pixels, may not perform well in snow- or ice-covered surfaces and water bodies, leading to higher uncertainty in the SDC500 for these areas compared to vegetated and arid regions. The current granularity used to derive the statistics of BRDF parameters is 100 km × 100 km, which may be too coarse and result in some loss of accuracy. Our future research should focus on exploring enhancements to the BRDF model and parameterization. In terms of prolonged gaps caused by cloudy days, there are still extreme cases that cannot be adequately filled using the proposed algorithm. Visual inspection has identified that these extreme cases occur mostly in mountainous terrain within tropical rainforest climate zones. In these areas, the topographic cloud formed by the terrain effect persists throughout the year, resulting in almost no clear satellite images. Another noteworthy case is the transition between snow-covered and snow-free states in high-latitude areas. As described in Sect. 3.5, the observations in a pixel's time series are divided into two groups, snow covered and snow free, and are processed separately. When these two results are combined, the final series may exhibit discontinuities around the switching date. Furthermore, the current algorithm struggles to robustly reconstruct images near the date of surface status change, presenting a challenging research topic.

In general, the process of image reconstruction involves altering the original values of sensor observations and replacing them with approximated values. Therefore, the reconstructed dataset is primarily intended to support large-scale statistical research or discrete applications such as land cover classification rather than small-scale or in-depth research such as model tuning and validation. In the past, people have tended to be overcautious when sharing the reconstructed dataset, being fearful of misleading data users. However, the demand for seamless data cubes has prompted us to share the SDC500 dataset and to provide an opportunity for evaluation and improvement by the wider public.