A Novel Strategy to Reconstruct NDVI Time-Series with High Temporal Resolution from MODIS Multi-Temporal Composite Products

: Vegetation indices (VIs) data derived from satellite imageries play a vital role in land surface vegetation and dynamic monitoring. Due to the excessive noises (e.g., cloud cover, atmospheric contamination) in daily VI data, temporal compositing methods are commonly used to produce composite data to minimize the negative influence of noise over a given compositing time interval. However, VI time series with high temporal resolution were preferred by many applications such as vegetation phenology and land change detections. This study presents a novel strategy named DAVIR-MUTCOP (DAily Vegetation Index Reconstruction based on MUlti-Temporal COmposite Products) method for normalized difference vegetation index (NDVI) time-series reconstruction with high temporal resolution. The core of the DAVIR-MUTCOP method is a combination of the advantages of both original daily and temporally composite products, and selecting more daily observations with high quality through the temporal variation of temporally corrected composite data. The DAVIR-MUTCOP method was applied to reconstruct high-quality NDVI time-series using MODIS multi-temporal products in two study areas in the continental United States (CONUS), i.e., three field experimental sites near Mead, Nebraska from 2001 to 2012 and forty-six AmeriFlux sites evenly distributed across CONUS from 2006 to 2010. In these two study areas, the DAVIR-MUTCOP method was also compared to several commonly used methods, i.e., the Harmonic Analysis of Time-Series (HANTS) method using original daily observations, Savitzky–Golay (SG) filtering using daily observations with cloud mask products as auxiliary data, and SG filtering using temporally corrected composite data. The results showed that the DAVIR-MUTCOP method significantly improved the temporal resolution of the reconstructed NDVI time series. It performed the best in reconstructing NDVI time-series across time and space (coefficient of determination (R 2 = 0.93~0.94) between reconstructed NDVI and ground-observed LAI). DAVIR-MUTCOP method presented the highest robustness and accuracy with the change of the filtering parameter (R 2 = 0.99~1.00, bias = 0.001, root mean square error (RMSE)


Introduction
Satellite-derived VIs have been being widely used in monitoring vegetation conditions and dynamics on regional or global scales [1,2]. The satellite-derived normalized difference vegetation index (NDVI) calculated from the spectral reflectance of near infrared (NIR) and visible red bands is one of the most used Vis [2][3][4]. Accurate NDVI data with high temporal (e.g., daily) and spatial resolution was preferred and sometimes necessary in many applications like forest disturbances [5,6] and vegetation phenology detection [7,8].
A Moderate Resolution Imaging Spectroradiometer (MODIS), with a high temporal and moderate spatial resolution, has been a key instrument aboard the NASA Terra and Aqua satellites providing a long time series of global satellite observations since 2000 [9][10][11][12]. MODIS provide a near-daily global coverage multispectral imagery of the entire earth surface for 36 spectral bands including the visible red and near infrared (NIR) bands with moderate spatial resolution (up to 250 m) [13]. Compared to the Advanced Very High Resolution Radiometer (AVHRR), the MODIS instrument has a higher spatial resolution and can provide improved global dynamics and processes occurring on the land [13][14][15]. Compared to other satellites with higher spatial resolution, such as Landsat and Sentinel, MODIS has higher temporal resolution [16,17]. As a result, MODIS data plays a vital role in a wide range of regional or global monitoring efforts of the earth system, considering its high temporal resolution, moderate spatial resolution, and long time-series of continuous daily observations (20 years) [1,18]. However, the noises in these products that include cloud cover and other atmospheric contamination, off-nadir viewing effects and shadow effects are the primary issues for accurate monitoring of daily vegetation dynamics [19]. Even after screening the data using the quality flag of the MODIS MOD09 product, 40% of the daily NDVI data were still potentially contaminated by residual sub-pixel clouds [19][20][21].
Considering the excessive noises remaining in the daily MODIS NDVI product, temporal compositing methods over 8-or 16-day intervals are commonly used to produce composite products with minimal cloud cover and atmospheric contamination [5,13]. Currently, MODIS provides daily, 8-day and 16-day products, e.g., MOD09GQ, MOD09Q1 and MOD13Q1 products, respectively (National Aeronautics and Space Administration (NASA), https://modis.gsfc.nasa.gov/data/, accessed on 20 March 2021). Temporally-composited MODIS VI products (e.g., MOD09Q1 and MOD13Q1) have been widely used in the previous studies, e.g., crop yield or gross primary production (GPP) estimation [22,23], vegetation phenology detection [9,12], and land surface change detection [6,24]. The maximum value composite (MVC) method is one of the common compositing methods designed to select the highest VI value from series of daily observations during a given time interval to represent the VI value for that time period [25]. MODIS VI composites, e.g., 16-day MOD13Q1 NDVI products, are produced using a refined version of the MVC technique called constrained view MVC (CV-MVC) method that reduces angular and sun-target-sensor variations [8]. The CV-MVC method is designed to select observations closest to nadir view zenith angles and can produce more consistent and accurate datasets [8,[26][27][28][29].
However, these composite MODIS products with decreased temporal resolution inevitably introduce uncertainties in two ways. First, only one value is selected in a given time interval and subtle, shorter-term VI change information within this time interval will be lost [7,8,30,31]. Second, the nominal observation date of the selected value is deter-mined as a fixed day within the composite time interval. As for MODIS composite products, the nominal observation date is the first day of the time interval. In reality, the actual acquisition date of this selected value can vary from the first to the last day of the composite time interval. Guindin-Garcia et al. [32] found that the temporal intervals (the period between two consecutive observations) reached 15 days for 8-day composites and 30 days for 16-day composites. This introduces uncertainty when comparing the value to the ground observations on a certain day [14,30]. In addition, previous research has shown that the adoption of the nominal date of composites introduced temporal errors that potentially caused appreciable changes in the trajectory of NDVI time series [8,33,34]. As a result, it is expected to be inadequate to correctly describe phenological patterns [34,35] nor to estimate rapidly changing biophysical characteristics [32,34]. For example, MODIS composite data using the first date in the interval was found to result in an earlier estimated start of growth [34][35][36].
To eliminate the temporal uncertainty of composite data, many previous studies suggested that this influence can be mitigated by adopting several temporal pattern strategies [32,34,35,[37][38][39]. For example, Testa et al. [40] used the median date of a composite period as the observation date based on the consideration of that the actual acquisition date was close to the center of the compositing period in the most years. While Guindin-Garcia et al. [32] found that the actual acquisition dates within the compositing period change without any predictable pattern, Wang and Zhu [38] found that the actual acquisition date of the NDVI value was usually later than the mean date of a composite period in spring and earlier in fall. The researchers, e.g., Thayn and Price [35], Guindin-Garcia et al. [32], Wang and Zhu [38], suggested using the actual acquisition date. Temporal correction with the actual acquisition date can eliminate the temporal shift and some errors of composite data. However, the degradation of temporal resolution during the compositing processing inevitably missed some key information, especially when monitoring the rapid changing events. This is the reason why daily VI time series data instead of temporal composite data were still suggested to be used in many previous studies [5][6][7][8]. While, in reality, composite data were more widely used in previous studies than daily data considering the data storage and computing cost, as well as the excessive noises in daily data which might introduce higher uncertainty than the degradation of temporal resolution [2].
Many smoothing and denoising methods have been developed and applied to reconstruct the daily VI data [2,41]. Local smoothing methods, for example the Savitzky-Golay (SG) filtering method, can exhibit strong fidelity, but it is difficult to interpret the excessive noise time-series data and it present less smoothness and spatial continuity [42,43]. Some global methods, such as asymmetric Gaussian and Harmonic Analysis of Time Series (HANTS) methods, have robust smoothness and tend to be less sensitive to the noise, but they are unable to describe more detailed, subtle changes of time-series curves with strong fidelity [2,42,44]. The widely used, upper envelope method assumes that the most types of noise are negatively biased and is designed to address low values that represent some forms of noise, but cannot detect noises that are positively biased [21] nor distinguish the actual vegetation change trajectory with lower to moderate levels of noise (e.g., a thin cirrus cloud) [45]. Cloud contamination and other atmospheric effects generally decrease NDVI values, while solar and viewing geometry variation through time introduce higher NDVI values and even a phase shift in the NDVI series [46][47][48][49]. No single method always outperforms all others under all these different situations to obtain a daily NDVI time series, because each of the methods have trade-offs in their respective approaches (e.g., removal of noise or preservation of the details of the NDVI temporal dynamics) [3,4,20,50].
Some studies suggested a combination of data from different sensors (e.g., data from both the Terra and the Aqua satellites) [7,21] or from microwave or geostationary satellitebased sensors [51,52]. However, the combinations of these different observations can introduce uncertainties considering the different atmospheric conditions at their varying observation times, different sensor characteristics and product generation algorithms [52,53]. For example, some studies applied cloud mask data in an attempt to develop a cloud-free dataset (e.g., [14]). However, the 250 m cloud mask indicator is based on the visible channel data only and still includes uncertainty. Luo et al. [54] reported that the standard MOD35 cloud mask product is inadequate for masking cloudy pixels and it could only identify bright and thick clouds but misses a large amount of other cloud types. Wilson et al. [55] found that spatial variability in the processing path applied in the MOD35 algorithm affects the likelihood of a cloudy observation by up to 20% in some areas. Wang et al. [56] reported that at least 9.1% of clouds were missed by the MOD35 product. Sun et al. [57] found that the universal accuracy of the MOD35 cloud mask algorithm was 50% for vegetated regions. In addition, noise introduced by other factors such as solar and viewing geometry variations are still difficult to eliminate by only using cloud mask data. As a result, a challenge still exists to reconstruct the daily NDVI with high quality that considers all these factors.
To address the above challenges, this study presents a novel strategy named the DAVIR-MUTCOP (Daily Vegetation Index Reconstruction based on Multi-Temporal Composite Products) method to reconstruct high-quality NDVI time-series data by combining the MODIS daily product (MOD09GQ) with MODIS 8-day or 16-day composite product (MOD09Q1 or MOD13Q1) as an example without any other auxiliary data. The core of the proposed DAVIR-MUTCOP method is using the NDVI time-series of composite product corrected by actual acquisition date to cross-check the original daily NDVI observations and select the daily observations that are not significantly contaminated. Most of the selected daily NDVI observations are not included in the temporal composite products, considering that the composite products select only one value in each temporal composite window. The proposed method is based on an assumption that the composite products which are less influenced by noise can generally reflect the temporal trends of the daily NDVI observations [8,[26][27][28][29], and it is reliable to screen the remaining effective daily NDVI observations in the daily product because these effective daily observations are restricted by the variation tendency of the composite data. The development of the DAVIR-MUTCOP method is presented in detail in the methodology section and the associated MATLAB codes and example datasets are provided as attachments. A comparison between DAVIR-MUTCOP and commonly used, existing methods for NDVI time series reconstruction was conducted for two study areas in the United States to demonstrate the capability and advantages of the DAVIR-MUTCOP method.

Study Area
The methods included in this study were applied in two study areas in the United States for validation and comparison purposes. The first study area includes three agri-  (Figure 1). The CSP sites locate in humid continental climate zone with severe winter and hot summer, but without a dry season. The mean annual temperature in the CSP sites is around 10 °C and the mean annual precipitation is about 790 mm. Given the large field sizes, the 250 m observations from MODIS near the center of these sites do not suffer from the problem of mixed pixels. CSP#1 was continuously planted in corn since 2001, while CSP#2 was planted in a corn (odd years)-soybeans (even years) rotation before 2009 and continuously planted to corn since 2009. CSP#3 has been planted in a corn (odd years)-soybean (even years) rotation since 2001. During the 2001 to 2012 study period, ground-based leaf area index (LAI) data of corn in these three sites were collected. The second study area includes forty-six AmeriFlux (AF) sites (AmeriFlux website: https://ameriflux.lbl.gov/, accessed on 20 March 2021) evenly distributed across the continental United States (CONUS) from 19° N latitude to 50° N latitude ( Figure 2). The climate of CONUS varies with the latitude and a range of geographic features, including mountains and deserts, ranging from the frost-free tropical of southernmost Florida to the cold climate zone in northeastern CONUS, from the rain drenched northwest coast to the drought-ridden deserts of the southwestern CONUS [58]. These forty-six sites were selected by two steps: (1) dividing the CONUS into a 5° × 5° grid, and (2) selecting the first site of each grid to make sure that the sites are generally evenly distributed across CO-NUS. In total, forty-six AmeriFlux ( Figure 2) sites were selected with various vegetation types, climates, and topographies.

Data Requirement
The MOD09GQ, MOD09Q1 and MOD13Q1 products were downloaded using the Google Earth Engine (GEE) platform (https://code.earthengine.google.com/, accessed on 20 March 2021). NDVI time series data at three temporal scales (daily, 8-day and 16-day) were obtained from MODIS MOD09GQ, MOD09Q1 and MOD13Q1 products, respectively. The MOD09GQ product provides the daily surface spectral reflectance (SSR) of visible red and NIR bands at a spatial resolution of 250 m. MOD09Q1 and MOD13Q1 products are temporally composited data over 8-and 16-day periods, respectively, to minimize cloud cover and atmospheric effects over a given time interval. The daily and 8-day composite NDVI time series data (NDVIGQ and NDVI09Q1) were calculated from the SSRs of MOD09GR and MOD09Q1 products (Equation (1)). 16-day composite NDVI time series data (NDVI13Q1) were directly derived from the layer of the MOD13Q1 product. MODIS products used the first day of the composite window as the nominal observation date of the composite value, instead of the actual acquisition date of the selected NDVI value for the composite period.
where and refer to the visible red and near-infrared band of MODIS.

Cloud Mask Product
The MOD35_L2 product provides a daily cloud mask indicator at 250 m spatial resolution. In order to minimize the influence of cloud cover, the cloud mask product was used to filter the cloud-free observations in this study. The MOD35_L2 product was not available on GEE platform and downloaded from the NASA website (https://modis.gsfc.nasa.gov/, accessed on 20 March 2021).

Methodology
At CSP sites, the LAI data were used to analyze the relationship between groundmeasured LAI and reconstructed daily NDVI data by different methods. The application of NDVI reconstruction methods for the AF sites allowed the reconstruction methods to be tested over various surface conditions compared with the CSP sites. The comparison of the methods' performances in the AF sites will further validate their advantages and disadvantages. Two existing NDVI reconstruction methods, HANTS and the combination of SG-filtering and Piecewise Cubic Hermite Interpolating Polynomial (PHCIP) interpolation, were used and compared with the DAVIR-MUTCOP method in this study. The steps to apply these methods are introduced in detail in the following subsections.

HANTS
The HANTS method was widely used to reconstruct the missing and biased satellitederived NDVI time series data by decomposing periodic time-dependent data into sum of sinusoids (see [59]). One of the advantages of the HANTS method is that it can easily be applied. The MATLAB source code of HANTS can be freely downloaded from the website (http://gdsc.nlr.nl/gdsc/en/tools/hants, accessed on 20 March 2021). The application of the HANTS method only needs the original NDVI data calculated from the MOD09GQ product. The HANTS method has five parameters, including the number of periodic terms in the NDVI series (N), the maximum number observations in a time series (normally one year) (M), the fit error tolerance (FET), the damping factor (Delta), and the degree of over determinedness (DOD). The parameters must be carefully set when using the HANTS method. However, according to previous studies [10,59], there is no objective guideline to determine the parameters. By testing the HANTS method for the CSP sites, the parameters were set to the values when the bias between estimated LAI and ground measured LAI is minimal. In the open source code, N and M were set as 365 and 10, respectively; FET was set as 0.05; Delta was set as 0.5; DOD was set as 5. The valid range of NDVI was defined as "0 to 1". The "low" option was selected as the flag indicating the direction of outliers with respect to the current curve.

SG Filtering
SG filtering was used to filter the original temporal composite NDVI data. To reconstruct daily NDVI data, the PHCIP method was then used to interpolate the filtered temporal composite NDVI data. According to the source of original NDVI data, three strategies were built in this study: (1) using the NDVI data calculated from acquisition-datecorrected MOD09Q1 product only (named SG8A hereafter); (2) using the NDVI data from the acquisition-date-corrected MOD13Q1 product only (named SG16A hereafter); and (3) combining the original daily NDVI data calculated from the MOD09GQ product and MOD35_L2 cloud mask product by filtering the original daily NDVI data with cloud mask data (named SG35L2 hereafter).
Before the application of the methods, the data of MOD09GQ, MOD09Q1 and MOD13Q1 was pre-processed using their quality flag layers as well as the median moving method [60]. As for the SG8A method, the outliers of the MOD09Q1 product were identified by the rule that the value of the outlier differed from the median by more than one standard deviation in a 7-point moving window [60]. As for the SG16A method, the outliers of the MOD13Q1 product were identified by using a 5-point moving window [40] and one standard deviation. As for the SG35L2 method, the outliers of the cloud-free data were identified by using a 25-point moving window according to Eklundh and Jösson [61] and one standard deviation.
The SG filtering and PHCIP interpolation were implemented in MATLAB. Parameters of SG filtering (i.e., temporal frame size and polynomial order) were set and tested according to the temporally composite windows of the input data. The values of these parameters were determined by a trial-and-error method using the data from the CSP study sites. The temporal frame sizes (Fsc represented the frame size of composite product, Fse represented the frame size of selected daily observations) for the SG8A, SG16A and SG35L2 methods were set as 3 (Fsc = 3), 2 (Fsc = 2) and 5 (Fse = 5), respectively. The polynomial order for all these methods was set as 1. The same parameter setting was used when these methods were applied for the AF sites. Then the sensitiveness of the reconstruction methods to the frame sizes of Fsc and Fse, was tested and analyzed by changing their value in the CSP and AF sites.

DAVIR-MUTCOP Method
The DAVIR-MUTCOP method used original daily NDVI data calculated from the MOD09GQ product, as well as the original composite NDVI data calculated from the MOD09Q1 product or the MOD13Q1 product to reconstruct daily NDVI time series. DAVIR-MUTCOP method is applied and evaluated in this study using the two standard MODIS compositing products that included: (1) combining daily MOD09GQ product and 8-day MOD09Q1 composite product (named DAVIR-MUTCOP8 hereafter) and (2) combining daily MOD09GQ product and a 16-day MOD13Q1 composite product (named DAVIR-MUTCOP16 hereafter). The development of the DAVIR-MUTCOP method is based on two assumptions that: (1) the temporal variation of NDVI data calculated from composite product can reflect the general variation of daily NDVI, and (2) the reasonable daily NDVI data derived from the MOD09GQ product fluctuates around the filtered composite NDVI time series with the actual acquisition date, considering that the applied scientific composite algorithm selected the value with the nominal highest quality within the composite window and therefore minimized the noise of the composite data.
The flowchart to implement DAVIR-MUTCOP method are shown in Figure 3 and presented in detail as follows: Step 1: Calculating the actual acquisition date of NDVI from composite product (NDVIC) using Equation (2). Similar to the previous study conducted by Guindin-Garcia et al. [32], the actual acquisition dates within the compositing period changed without any predictable pattern in this study (not shown). The original composite NDVI time series was adjusted to its actual acquisition date obtaining a new time series NDVIAD, in which the observation dates were not equidistant. SG8A and SG16A methods were also implemented using the NDVIAD.
where m is the sequence number of composite NDVIC in a year, n is the sequence number of original daily NDVIGQ (derived from MOD09GQ product) in a year, tm is the actual acquisition date of mth NDVIc during a period of one year, tn is the date of n th NDVIGQ during a period of one year, D is the time interval of composite product. When the composite product was MOD09Q1, D was set as 8. When the composite product was MOD13Q1, D was set as 16.
Step 2: Denoising and filtering the NDVIAD data. It is assumed that a biased NDVIAD is usually lower than the true data. If the value of mth NDVIAD was lower than its immediate neighbors ((m − 1)th and (m + 1)th NDVIAD), the mth NDVIAD was temporarily deleted.
Step 3: The remaining NDVIAD data were filtered using SG filtering with frame size (Fsc) set as 3 (2) for MOD09Q1 (MOD13Q1) product and polynomial order set as 1. The filtered NDVIAD was named NDVIF. The sensitiveness of the method to Fsc in reconstructing the daily NDVI time series was discussed in Section 4.4.2.
Step 4: Interpolation for daily NDVI series (named NDVIP) from NDVIF using PHCIP interpolation algorithm. The normalized NDVIP (named NNDVIP) and the absolute value of the first derivative of NDVIP (named DNDVIP) were also computed. For the original NDVIAD obtained in Step 1, the NDVIAD of i th day (day of year) should be considered as the reasonable data if it satisfied Equation (3). After the new reasonable NDVIAD time series was obtained, we returned to Step 3. When the NDVIAD time series did not change, the loop was ended.
Step 5: Finding out the effective NDVIGQ (named ENDVIGQ). The effective NDVIGQ is expected to fluctuate around the trajectory of the reasonable NDVIAD. Based on the final NDVIP, DNDVIP, and NNDVIP in step 4, the ENDVIGQ was selected by using Equation (4).
Step 6: Reconstruction of the daily NDVI time series. SG filtering (frame size (Fse) = 5, polynomial order = 1) was used to filter the ENDVIGQ and filtered data were named NDVIFE. The NDVIFE was then interpolated by PHCIP algorithm to obtain the final reconstructed daily NDVI time series (named NDVIRD). The sensitiveness of the method to Fse in reconstructing the final NDVIRD was discussed in Section 4.4.3.

Evaluation of the Method's Performance
In CSP sites, the relationship between LAI and reconstructed NDVI, as well as the relationship between the NDVI data reconstructed by different methods, were used to quantitatively validate the accuracy and analyze the uncertainty and error of different reconstruction strategies in terms of the bias (Equation (5)), root mean square error (RMSE) (Equation (6)), slope and determination coefficient ( where K was the number of NDVI data, NDVIm1,i and NDVIm2,i were the NDVI reconstructed by method 1 and method 2.

The MODIS Products
In order to study the MODIS products, the NDVI time-series of the second CSP site in 2009 was taken as an example for analysis ( Figure 4). The original MOD09GQ product shows extensive, noise-contaminated NDVI data. It is difficult to interpret time-series data with excessive noise features and reconstruct the time-series data with high smoothness and spatial continuity [42,43]. Daily NDVI observations filtered by MOD35_L2 cloud mask data still included considerable noises (Figure 4a,d). This is probably explained by the facts that: (1) MOD35_L2 product still includes error and uncertainty considering that sub-pixel clouds, thin clouds and cloud shadows are difficult to be fully detected or removed [55,56]. Thus, SG35L2 method tends to underestimate NDVI values, especially during or approaching the NDVI peak period (Figure 4), and (2) in addition to cloud cover, other types of noise that exist in the dataset, including other atmospheric contamination and off-nadir viewing effects, can increase or decrease NDVI values [47,48]. 8-and 16-day composite data had much less noise than daily NDVI data (Figure 4b,c), because the CV-MVC method generally selects the representative NDVI with higher accuracy in the given time interval to minimize the influence of the noises in daily time series data, though there were still remaining data noises in composite data, especially for the 8-day composite data product (Figure 4e,f).
However, obvious temporal shifts were observed in 8-and 16-day NDVI values compared to the original daily observations, especially during the period when NDVI values changed quickly (Figure 4b,c), as the nominal dates of the composite products MOD09Q1 and MOD13Q1 is defined as the first day of the composite window. While during the green-up stage, the CV-MVC method is expected to select a higher NDVI value from late in the composite period than what would be observed on the first day of the composite period, as the NDVI is increasing across the 8 or 16-day window. During the senescence phase when the NDVI is decreasing, the CV-MVC method is expected to select a NDVI value towards the beginning of the composite period because the NDVI would gradually decrease as the days in the composite period progressed. However, usually, a decreased NDVI value observed after the first day is always selected considering that the factor, such as high-frequent noise or large sensor viewing angles, might affect the observations of the first day and even the nearby days as well. In addition, the "horizontal shift" (temporal shift) can introduce "vertical shift" (NDVI value shift). Thus, the composite data commonly tends to overestimate the NDVI values before the NDVI peak and underestimates the NDVI values after the NDVI peak. This is also the reason to temporally correct the composite data by using actual acquisition dates in the DAVIR-MUTCOP method, as shown in the flowchart in Figure 3.

The Temporal Resolution of Different Reconstruction Stratiges
Temporal composition of the NDVI time series reduced the data noises to some extent. However, it also significantly sacrificed the temporal resolution, especially for 16day composite data. For example, there might be few or even no clear observations in composite data during the period when NDVI increased rapidly from the bottom to the peak (Figure 4c,e,f). The effect of the data gap combined with the above-discussed temporal shift due to nominal observation data might introduce obvious changes and uncertainties in NDVI time-series trajectories and NDVI values when directly using MODIS composite data to reconstruct NDVI time-series curves.
As shown in Figure 4b,c,e,f, the DAVIR-MUTCOP method corrected the temporal shift of the composite data products and included more validate daily observations (black dots) to increase the temporal resolution as well, which were significantly helpful for accurate time-series NDVI reconstruction. It is shown from Figure 5 that the DAVIR-MUTCOP method significantly promotes the temporal resolution of NDVI observations in CSP sites, especially during the growing season from April to October. The 8-and 16day composite data had about 3.5 and 1.5 clear observation per month on average, respectively (i.e., 9-day and 20-day temporal resolution, respectively), while the DAVIR-MUTCOP method obtained at least 10 clear observations per month at the growing season (i.e., <3-day temporal resolution). Besides, it is shown in Figure 4 that the DAVIR-MUTCOP method provided a less biased and a more smoothed daily NDVI time series than SG35L2 method.
In addition, as shown as the averaged clear observations per month in Figure 6, the temporal resolution of the reconstructed NDVI time series was obviously improved by the DAVIR-MUTCOP method in the six selected AF sites with different landcover types. The improvement of the temporal resolution is particularly helpful to monitor the rapid vegetation changes during the periods when the onset or end of the growing season occurs. For example, the number of clear original observations included in DAVIR-MUTCOP method was about 3 and 6 times of that in MODIS 8-and 16-day composite data, respectively, in April, May, October and November at P17 site covered by deciduous broadleaf forests (Figure 6c). Leaf-expansion and leaf-fall of the deciduous broadleaf forest usually occur in these four months at the P17 site. It might be difficult to detect these two stages directly from the composite data, especially for 16-day composite data with < 2 clear observations per month, considering that the vegetation changed rapidly and the duration of these period was short. For example, Nagai et al. [7] reported that the duration of leaf-expansion and leaf-fall varied from 11 to 16 days during the period from 2004 to 2010 at the Takayama site.

The Comparison of the NDVI Time-Series Curves Reconstructed by Different Stratiges
The temporal variations of NDVI time-series data in the three CSP sites and forty-six AF sites reconstructed by SG8A, SG16A, HANTS, SG_35L2, DAVIR-MUTCOP8 and DAVIR-MUTCOP16 methods were presented in Figures S1 and S2, respectively. HANTS achieves high accuracy in reconstructing some of the NDVI time series. However, daily NDVI data in some years reconstructed by the HANTS method suffered from significant fluctuation. Especially, the missing data and noisy NDVI values lead to the failure of this method. For example, the reconstructed time series of the third CSP site (CSP#3) in 2001 and 2006 presented many local fluctuations and peaks not reflecting the terrestrial biotic dynamics ( Figure S1). In addition, the accuracy of HANTS method was relatively sensitive to its parameter settings. The accuracy might significantly decrease when it is applied in other areas without parameter optimization. It is difficult to automatically optimize the parameters in different areas and years [10,59] without ground-based reference data. Therefore, the HANTS method was not further discussed in this study.
The SG35L2 method used MOD35_L2 cloud mask data to screen cloudy free data. However, the reconstructed time series still included some local fluctuations and noises when compared to the composite data ( Figures S1 and S2). The 41th AF site (P41) was also selected to be presented in Figure 7 as an example for further analysis. The SG35L2 method also presented many local fluctuations and troughs even with Fse increased from 5 to 15. The time series reconstructed by the SG8A method were also jagged and exhibited by various temporally localized peaks and troughs due to the data noises. While both the DAVIR-MUTCOP8 and DAVIR-MUTCOP16 methods using the same Fsc with the SG8A/16A method and the same Fse with SG35L2 provided more smooth and consistent results. SG16A provided more smooth results than the SG8A method, which seemed to be similar to that of the DAVIR-MUTCOP method. 16-day composite data with a longer compositing window commonly includes less noise. However, compositing data sacrifices the temporal resolution. It is more degraded with the increase of the temporal compositing window.

The Relationship between Ground-Observed LAI and Reconstructed NDVI
In order to quantitatively validate the reconstructed NDVI time series, ground-observed plant biophysical parameter LAI were used to build its relationships with NDVI values reconstructed by SG8A/16A, DAVIR-MUTCOP8/16 and SG35L2 methods at the three CSP sites. As shown in Figure 8, DAVIR-MUTCOP methods obtained similar performance with SG8/16A methods in terms of coefficient of determination (R 2 ) at CSP sites. This is possibly explained by two facts: (1) the composite product in CSP sites has high data quality without frequent cloud contamination compared to most of the other areas across CONUS, (2) the time series of CSP sites do not suffer from mixed pixel problem, consequently regular shape of the time-series curves minimized the uncertainty and error introduced by the degradation of temporal resolution and data interpolation during the periods between two consecutive observations.
But it is interesting that the fitting lines for ground observed LAI and NDVI values reconstructed by DAVIR-MUTCOP8 and DAVIR-MUTCOP16 method were almost overlayed (Figure 8a). An obvious difference was observed between those of SG8A and SG16A as shown in Figure 8b. Similarly, the scatters plot (Figure 9b) of SG8A-reconstructed NDVI against SG16A-reconstructed NDVI was more spread than that (Figure 9a) of DAVIR-MUTCOP8-reconstructed NDVI against DAVIR-MUTCOP8-reconstructed NDVI (R 2 = 0.98 against 0.99, bias = 0.004 against 0.001, RMSE = 0.034 against 0.020), especially when NDVI changed rapidly (NDVI = 0.4~0.8). This indicates that the DAVIR-MUTCOP method can obtain more stable NDVI reconstruction unaffected by the size of compositing window in CSP sites.

The Sensitiveness of the Reconstructed NDVI to Fsc
In order to further validate the robustness of SG8A, SG16A and DAVIR-MUTCOP methods, the sensitiveness of the reconstruction results to Fsc in different methods was tested. Due to the degradation of the composite data, the rapid or unexpected changes might be confused with rapid drop caused by data noises. The reconstruction results might be sensitive to filtering parameters. For example, as shown in Figure 10a,b and Figure 11a,b, SG8A and SG16A methods were sensitive to the adopted Fsc in SG filter when applied in the 39 th and 10 th AF sites (P39 and P10). Some data noises might still remain when small Fsc was used, while bigger Fsc might even cause the temporal shift of the reconstructed curves for both the SG8A and the SG16A method. The sensitiveness of the reconstruction to the filtering parameters combined with the data noises in the composite data significantly increased the uncertainty of the reconstructed NDVI time-series data, especially for 16-day composite data or in the period when NDVI changed rapidly. In addition, the SG16A method with a longer temporal compositing window was shown to be more sensitive to the parameters, presenting less consistent NDVI time-series trajectories (Figures 10a,b and 11a,b) and more spread scatters of the NDVI values with the change of Fsc rather than the SG8A method (Figure 10e,f). The DAVIR-MUTCOP method, including both original daily observations and composite products, was shown to be insensitive to the adopted Fsc compared to the SG8A and SG16A methods, presenting consistent NDVI time series in the growing season with different compositing windows and different Fsc in an SG filter (Figures 10 and 11).

The Sensitiveness of the Reconstructed NDVI to Fse
The sensitiveness of the DAVIR-MUTCOP and SG8A/16A methods to Fsc was discussed above. The sensitiveness of DAVIR-MUTCOP and SG35L2 methods to Fse was then further analyzed. As shown in Figure 7, the SG35L2 method was significantly affected by Fse in the SG filter due to the existence of the remaining noises. The scatter plot in Figure 12 also showed that the results derived from DAVIR-MUTCOP method with different Fse were more consistent with each other than those from SG35L2.
Overall, the performance of the method that only used original daily or composited data varied with the temporal compositing window as well as the parameters adopted in the data smoothing method. For example, 16-day composite data had superior performance at the P41 site, but inferior performance at P39 and P10 site than 8-day composite data. As for the SG16A method, smaller Fsc was suggested in the P39 site where larger Fsc might result in temporal shift of the time-series trajectories (Figure 10b), while larger Fsc was suggested in the P10 site where smaller Fsc might fail to denoise the data ( Figure  11b). Therefore, usually, it is difficult to find a method or parameter that is generally applicable to remove the noise, retain the short biotic fluctuations and rebuild the accurate growth curves meanwhile in regional applications [2]. This is also the reason why many previous studies that evaluated the different NDVI time-series data smoothing and reconstruction methods stated that no single method always outperforms all others under all these different situations [2][3][4], considering the data noise, the variation and complexity of the vegetation changes under different landcover types at regional scale, as well as the limitation of each method. While the DAVIR-MUTCOP method is more robustness and appliable to reconstruct NDVI time series at the regional scale than other methods.

The Choice of Temporal Compositing Window
Usually, 8-day composite data with higher temporal resolution is more widely used and preferred to reconstruct NDVI time series in the previous studies [9,23]. However, in the cloud-prone areas, 8-day composite data might be still noisy when clear observations might not exist within the 8-day window. Thus, there might be a trade-off: higher temporal resolution or less data noise, when choosing the temporal compositing window of the composite data.
As for the SG8A and SG16A methods, on one hand, the NDVI values in the time series with actual acquisition dates were not equidistant and the period between two consecutive observations varied widely and reached up to 15 days for 8-day composites and 30 days for 16-day composites [32]. The longer interval between two consecutive observations makes it more difficult to correctly reconstruct daily NDVI time series with a continuous missing value. However, for the other hand, 16-day composite data commonly includes less noise, as it has a longer compositing window than 8-day composite data. Therefore, it has a greater possibility of obtaining a higher and clearer observation within the compositing window. It indicates that in practical application, it is hard to choose the composite product for SG8A or SG16A method, which might limit their applications at the regional to the global scale.
However, as shown in Figures 7-12, the DAVIR-MUTCOP16 method had similar a performance to the DAVIR-MUTCOP8 method in terms of accuracy and robustness. In addition, the temporal resolution of reconstructed data was not sacrificed by choosing 16day composite data instead of 8-day composite data, as the DAVIR-MUTCOP method used original daily data as well. As shown in Figures 5 and 6, the NDVI time series reconstructed by DAVIR-MUTCOP16 method generally obtained comparable temporal resolution with the DAVIR-MUTCOP8 method. While the 16-day composite data product with a longer compositing window is less influenced by data noises, it indicates that 16-day composite data MOD13Q1 rather than 8-day composite data MOD09Q1 is suggested when the growth and changes of the studied vegetation is slowly or smooth, which might favor its applications in cloud-prone areas. While 8-day composite data of MOD09Q1 is suggested when the changes of the studied vegetation might be rapid and unexpected.

The Choice of the Filtering/Fitting Method and Parameters
Many smoothing and denoising methods have been developed and applied to reconstruct the NDVI time series data [2,41], e.g., SG filtering method, HANTS methods, logistic method. The widely used SG filtering method was applied in the proposed DAVIR-MUTCOP method in this study to denoise the composite data and to reconstruct the continuous NDVI time series based on the selected daily NDVI observations as well. However, many other smoothing methods also have potential to be used to replace SG filtering in the DAVIR-MUTCOP method.
As for the parameters of different data smoothing methods, it is usually difficult to optimize a set of parameters for the regional applications. Both local methods are usually sensitive to parameters setting [2]. While the DAVIR-MUTCOP method was shown to be obviously less sensitive to the variation of the frame size than other reconstruction strategies when it varied in a reasonable range. The frame sizes of SG filter (Fsc and Fse) used in DAVIR-MUTCOP8 (Fsc = 3~5, Fse = 5~10) and DAVIR-MUTCOP16 (Fsc = 2~3, Fse = 5~10) methods are also suggested in the future applications considering the robust performance of the DAVIR-MUTCOP8/16 method when it was applied in different types of vegetation covers.

The Potential Application of the DAVIR-MUTCOP Method in the Future
NDVI time series data with high temporal resolution plays a very important role in many applications, e.g., remote sensing of phenology, agriculture, and forest disturbance. Compared to the existed approaches presented in this study, the DAVIR-MUTCOP method has been proven to be an effective and robust way to reconstruct NDVI time series under various conditions. The DAVIR-MUTCOP method achieved a high level of robustness in reconstructing daily data time series without changing the parameters over time and space, as illustrated by the results at the CSP and AF sites. In addition, although the DAVIR-MUTCOP method was only applied and validated for NDVI time series reconstruction based on MODIS data, it proposed a universal way to reconstruct daily time series of other VIs (e.g., Enhanced Vegetation Index (EVI), Wide Dynamic Range Vegetation Index (WDRVI)) from MODIS. It is also expected to be applied to VI time series datasets from other operational sensors, such as Advanced Very High Resolution Radiometer (AVHRR) and Visible Infrared Imaging Radiometer Suite (VIIRS). Especially, the land science of VIIRS is expected to build and expand on the heritage of land science from the NOAA AVHRR and EOS MODIS. VIIRS data will be used to expand upon the MODIS applications (NASA, https://earthdata.nasa.gov/, accessed on 20 March 2021). VIIRS has similar products with MODIS, such as daily surface reflectance products (e.g., VNP09GA), 8-day composite surface reflectance products (e.g., VNP09A1) and 16-day vegetation indices products (e.g., VNP13A1), etc.

The Limitation of the DAVIR-MUTCOP Method
The DAVIR-MUTCOP method can generally reconstruct MODIS NDVI time series with high accuracy and high temporal resolution. It still has limitations when the composite data selected from the temporally composite window still contains continuous noise (e.g., prolonged cloudy periods), i.e., the accuracy of temporally composite products might affect the accuracy of reconstruction to some extent. Fortunately, the use of composite data with a longer compositing window, e.g., 16-day, does not affect the application of the DAVIR-MUTCOP method and a longer compositing window can increase the reliability and availability of noise-free data in the given time interval.

Conclusions
Although the original daily observations derived from the satellites include massive incorrect values, there are still many observations with high quality that are not selected in the temporally composite product. The core of the proposed DAVIR-MUTCOP method is making full use of the all effective daily NDVI observations as possible by using the composite NDVI time series with actual acquisition date to screen the daily observations. The DAVIR-MUTCOP and existed HANTS, SG8A, SG16A, SG35L2 methods were applied and evaluated in two study areas in the USA.
SG8A and SG16A methods with temporal correction by the actual acquisition date minimized the temporal shift of observations. However, the unequal distances between the two observations with an actual acquisition date increased the difficulties to reconstruct the daily NDVI time series, considering the distance can vary from one to nearly double the compositing window. The HANTS method is sensitive to the outliers and continuous contamination, which results in the failure of NDVI reconstruction. The parameters' setting and optimization of the HANTS method is also a problem, especially for regional or global application. Cloud mask data products improve the performance of the SG35L2 method compared to the SG1 method. However, cloud mask data products also include uncertainty; for example, the undetected cloud. In addition, other existed factors except for cloud cover can also influence daily observations.
The DAVIR-MUTCOP method combined both daily and composite products and presented the best performance across time and space than other methods applied in this study in terms of accuracy, robustness and temporal resolution. The success of the DAVIR-MUTCOP method not only minimizes excessive noise of the original daily observation and corrects the temporal shift of the composite data by cross-checking the daily and composite data, but also reserves more daily observations with high quality to improve the temporal resolution. Note that the DAVIR-MUTCOP method assumes that the filtered composite data are relatively reliable. In the cloud-prone areas, the DAVIR-MUTCOP method might still fail to reconstruct the daily VI time series when the composite observations, e.g., MOD09Q1 or MOD13Q1 products, are still continuously noisy and unreliable. In these areas, longer temporal composite windows are suggested to be adopted to obtain more reliable composite data.
The DAVIR-MUTCOP method provided a universal way to reconstruct high-temporal-resolution time series of other VIs or from other operations, e.g., AVHRR and VIIRS. Our next work will include: (1) further quantificationally validating the DAVIR-MUTCOP method under different frequency of cloud cover and different change rates of underlying surface by comparing to the near surface optical sensors; (2) validating this method in related applications, e.g., phenology detection, forest disturbance, plant physiological parameter estimation at the reginal scale.