Next Article in Journal
Shoreline Changes along Northern Ibaraki Coast after the Great East Japan Earthquake of 2011
Previous Article in Journal
Crustal Strain and Stress Fields in Egypt from Geodetic and Seismological Data
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

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

1
College of Resources and Environment, Huazhong Agricultural University, Wuhan 430070, China
2
Center for Advanced Land Management Information Technologies, University of Nebraska-Lincoln, 3310 Holdrege St, Lincoln, NE 68583, USA
3
School of Environmental Studies & State Key Laboratory of Biogeology and Environmental Geology, China University of Geosciences, Wuhan 430074, China
4
State Key Laboratory of Information Engineering in Surveying, Mapping and Remote Sensing, Wuhan University, Wuhan 430079, China
5
Guangxi Key Laboratory for Spatial Information and Geomatics, Guilin University of Technology, Guilin 541004, China
6
Changjiang River Scientific Research Institute, Changjiang River Water Resources Commission, Wuhan 430010, China
7
Changjiang Institute of Survey, Planning, Design and Research, Wuhan 430010, China
8
Guanxi Key Laboratory of Water Engineering Materials and Structures, Guanxi Institute of Water Resources Research, Nanning 530023, China
*
Author to whom correspondence should be addressed.
Submission received: 20 February 2021 / Revised: 31 March 2021 / Accepted: 31 March 2021 / Published: 5 April 2021
(This article belongs to the Section Remote Sensing in Agriculture and Vegetation)

Abstract

:
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 (R2 = 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 (R2 = 0.99 ~ 1.00, bias = 0.001, root mean square error (RMSE) = 0.020). Only MODIS data were used in this study; nevertheless, the DAVIR-MUTCOP method proposed a universal and potential way to reconstruct daily time series of other VIs or from other operational sensors, e.g., AVHRR and VIIRS.

1. 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 determined 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 satellite-based 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.

2. Study Area and Data Requirement

2.1. 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 agriculture experimental sites (#1: 41°9′54.2″ N, 96°28′35.9″ W; #2: 41°9′53.5″ N, 96°28′12.3″ W; #3: 41°10′46.8″ N, 96°26′22.7″ W, 49 to 65 ha) located near Mead, Nebraska at the University of Nebraska-Lincoln’s (UNL) Eastern Nebraska Research and Extension Center (ENREC) as part of UNL’s Carbon Sequestration Program (CSP, http://csp.unl.edu/public/, accessed on 20 March 2021) (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 CONUS. In total, forty-six AmeriFlux (Figure 2) sites were selected with various vegetation types, climates, and topographies.

2.2. Data Requirement

2.2.1. Multi-Temporal NDVI Data

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.
NDVI =   ρ n i r ρ r / ρ n i r + ρ r
where ρ r and ρ n i r refer to the visible red and near-infrared band of MODIS.

2.2.2. 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).

3. Methodology

At CSP sites, the LAI data were used to analyze the relationship between ground-measured 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.

3.1. HANTS

The HANTS method was widely used to reconstruct the missing and biased satellite-derived 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.

3.2. 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-date-corrected 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.

3.3. 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.
t m = t n , i f N D V I c , m = N D V I G Q , n , D ×   m 1   +   1 n D × m
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 ith 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.
N D V I P , i 0.05 C R a t d · D N D V I P , i <   N D V I A D , i < N D V I P , i + 0.05 + C M F R · N N D V I P , i + C R a t u · D N D V I P , i
where CMFR, CRatu, CRatd were the multiplying factor (CMFR = 0.05, CRatu = 20, CRatd = 20).
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).
N D V I P , i 0.025 G Q R a t d · D N D V I P , i <   E N D V I G Q , i < N D V I P , i + 0.05 + G Q M F R · N N D V I P , i + G Q R a t u · D N D V I P , i
where CMFR, CRatu, CRatd were the multiplying factor (CMFR = 0.05, CRatu = 10, CRatd = 5).
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.

3.4. 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 (R2).
b i a s = i = 1 K N D V I m 1 , i N D V I m 2 , i K
R M S E = i = 1 K N D V I m 1 , i N D V I m 2 , i 2 K 1
where K was the number of NDVI data, NDVIm1,i and NDVIm2,i were the NDVI reconstructed by method 1 and method 2.

4. Results

4.1. 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.

4.2. 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 16-day 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 16-day 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.

4.3. 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.

4.4. The Robustness of the Different Reconstruction Stratiges

4.4.1. 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 (R2) 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 (R2 = 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.

4.4.2. 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 39th and 10th 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 (Figure 10a,b and Figure 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 (Figure 10 and Figure 11).

4.4.3. 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.

5. Discussion

5.1. 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 Figure 7, Figure 8, Figure 9, Figure 10, Figure 11 and Figure 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 16-day composite data instead of 8-day composite data, as the DAVIR-MUTCOP method used original daily data as well. As shown in Figure 5 and Figure 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.

5.2. 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.

5.3. 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.

5.4. 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.

6. 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.

Supplementary Materials

The following are available online at https://0-www-mdpi-com.brum.beds.ac.uk/article/10.3390/rs13071397/s1. Figure S1: The temporal variation of NDVI data originally derived from MOD09GQ product (dots) and reconstructed by SG8/16A, HANTS, SG35L2, and DAVIR-MUTCOP8/16 methods in the three CSP sites (lines). Figure S2: The temporal variation of NDVI data originally derived from MOD09GQ product (dots) and reconstructed by SG8/16A, SG35L2, and DAVIR-MUTCOP8/16 methods in the forty-six AF sites (lines). Code S1: DAVIR-MUTCOP method’s Matlab code: “DAVIR_MUTCOP_sourcecode.m”. Data S1: Example datasets: “OriginalCompositeData.txt” (from MOD09Q1 product) and “OringinalDailyData.txt” (from MOD09GQ product).

Author Contributions

Conceptualization, L.Z. and S.H.; methodology, L.Z. and S.H.; software, L.Z. and S.H.; validation, L.Z., B.D.W. and S.H.; formal analysis, L.Z. and S.H.; investigation, S.H., X.Z., G.P., D.X., R.W., R.M. and W.W.; resources, B.D.W. and G.Z.; data curation, L.Z.; writing—original draft preparation, L.Z.; writing—review and editing, L.Z., B.D.W. and S.H., X.Z., G.Z., D.X., R.M.; visualization, S.H.; supervision, B.D.W. and S.H.; project administration, L.Z. and S.H.; funding acquisition, L.Z. All authors have read and agreed to the published version of the manuscript.”

Funding

This work was supported by the National Nature Science Foundation of China program (Grant No. 41901353), the Fundamental Research Funds for the Central Universities, China University of Geosciences (Wuhan) (No. 162301202679), Huazhong Agricultural University (No. 2662019QD054), the National Key Research and Development Program of China (Grant No. 2017YFC1502406-03), the Open Research Fund of Guangxi Key Laboratory of Water Engineering Materials and Structures, Guangxi Institute of water resources research (Grant No. GXHRI-WEMS-2019-03) and National Key Research and Development Program of Guangxi (Grant No. 2019AB20009).

Institutional Review Board Statement

Not applicable.

Informed Consent Statement

Not applicable.

Data Availability Statement

The data presented in this study are available in the Supplementary Materials.

Acknowledgments

The authors would like to thank the anonymous reviewers for their valuable comments to improve the quality of the paper. We also acknowledge for the data support from Carbon Sequestration Program of University of Nebraska-Lincoln’s (UNL) (CSP, http://csp.unl.edu/public/, accessed on 20 March 2021) and AmeriFlux (AF) sites (AmeriFlux website: https://ameriflux.lbl.gov/, accessed on 20 March 2021).

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. Zhang, Y.; Song, C.; Band, L.E.; Sun, G.; Li, J. Reanalysis of global terrestrial vegetation trends from MODIS products: Browning or greening? Remote Sens. Environ. 2017, 191, 145–155. [Google Scholar] [CrossRef] [Green Version]
  2. Zeng, L.; Wardlow, B.D.; Xiang, D.; Hu, S.; Li, D. A review of vegetation phenological metrics extraction using time-series, multispectral satellite data. Remote Sens. Environ. 2020, 237, 111511. [Google Scholar] [CrossRef]
  3. Cai, Z.; Jönsson, P.; Jin, H.; Eklundh, L. Performance of Smoothing Methods for Reconstructing NDVI Time-Series and Estimating Vegetation Phenology from MODIS Data. Remote Sens. Basel 2017, 9, 1271. [Google Scholar] [CrossRef] [Green Version]
  4. Hird, J.N.; Mcdermid, G.J. Noise reduction of NDVI time series: An empirical comparison of selected techniques. Remote Sens. Environ. 2009, 113, 248–258. [Google Scholar] [CrossRef]
  5. De Beurs, K.M.; Townsend, P.A. Estimating the effect of gypsy moth defoliation using MODIS. Remote Sens. Environ. 2008, 112, 3983–3990. [Google Scholar] [CrossRef]
  6. Jin, S.; Sader, S.A. MODIS time-series imagery for forest disturbance detection and quantification of patch size effects. Remote Sens. Environ. 2005, 99, 462–470. [Google Scholar] [CrossRef]
  7. Nagai, S.; Saitoh, T.M.; Suzuki, R.; Nasahara, K.N.; Lee, W.; Son, Y.; Muraoka, H. The necessity and availability of noise-free daily satellite-observed NDVI during rapid phenological changes in terrestrial ecosystems in East Asia. For. Sci. Technol. 2011, 7, 174–183. [Google Scholar] [CrossRef]
  8. Narasimhan, R.; Stow, D. Daily MODIS products for analyzing early season vegetation dynamics across the North Slope of Alaska. Remote Sens. Environ. 2010, 114, 1251–1262. [Google Scholar] [CrossRef]
  9. Zeng, L.; Wardlow, B.D.; Wang, R.; Shan, J.; Tadesse, T.; Hayes, M.J.; Li, D. A hybrid approach for detecting corn and soybean phenology with time-series MODIS data. Remote Sens. Environ. 2016, 181, 237–250. [Google Scholar] [CrossRef]
  10. Zhou, J.; Jia, L.; Menenti, M. Reconstruction of global MODIS NDVI time series: Performance of Harmonic ANalysis of Time Series (HANTS). Remote Sens. Environ. 2015, 163, 217–228. [Google Scholar] [CrossRef]
  11. Tan, B.; Morisette, J.T.; Wolfe, R.E.; Gao, F.; Ederer, G.A.; Nightingale, J.; Pedelty, J.A. An Enhanced TIMESAT Algorithm for Estimating Vegetation Phenology Metrics from MODIS Data. IEEE J. Sel. Top. Appl. Earth Obs. Remote Sens. 2011, 4, 361–371. [Google Scholar] [CrossRef]
  12. Zhang, X.; Friedl, M.A.; Schaaf, C.B.; Strahler, A.H.; Hodges, J.C.F.; Gao, F.; Reed, B.C.; Huete, A. Monitoring Vegetation Phenology Using MODIS Time-Series Data. Remote Sens. Environ. 2003, 84, 471–475. [Google Scholar] [CrossRef]
  13. Justice, C.O.; Townshend, J.R.G.; Vermote, E.F.; Masuoka, E.; Wolfe, R.E.; Saleous, N.; Roy, D.P.; Morisette, J.T. An overview of MODIS Land data processing and product status. Remote Sens. Environ. 2002, 83, 3–15. [Google Scholar] [CrossRef]
  14. Spruce, J.P.; Sader, S.; Ryan, R.E.; Smoot, J.; Kuper, P.; Ross, K.; Prados, D.; Russell, J.; Gasser, G.; McKellip, R.; et al. Assessment of MODIS NDVI time series data products for detecting forest defoliation by gypsy moth outbreaks. Remote Sens. Environ. 2011, 115, 427–437. [Google Scholar] [CrossRef]
  15. Hansen, M.C.; DeFries, R.S.; Townshend, J.R.G.; Sohlberg, R.; Dimiceli, C.; Carroll, M. Towards an operational MODIS continuous field of percent tree cover algorithm: Examples using AVHRR and MODIS data. Remote Sens. Environ. 2002, 83, 303–319. [Google Scholar] [CrossRef]
  16. Gao, F.; Masek, J.; Schwaller, M.; Hall, F. On the blending of the Landsat and MODIS surface reflectance: Predicting daily Landsat surface reflectance. IEEE Trans. Geosci. Remote 2006, 44, 2207–2218. [Google Scholar]
  17. Hilker, T.; Wulder, M.A.; Coops, N.C.; Linke, J.; McDermid, G.; Masek, J.G.; Gao, F.; White, J.C. A new data fusion model for high spatial- and temporal-resolution mapping of forest disturbance based on Landsat and MODIS. Remote Sens. Environ. 2009, 113, 1613–1627. [Google Scholar] [CrossRef]
  18. Tian, F.; Fensholt, R.; Verbesselt, J.; Grogan, K.; Horion, S.; Wang, Y. Evaluating temporal consistency of long-term global NDVI datasets for trend analysis. Remote Sens. Environ. 2015, 163, 326–340. [Google Scholar] [CrossRef]
  19. Motohka, T.; Nasahara, K.N.; Murakami, K.; Nagai, S. Evaluation of Sub-Pixel Cloud Noises on MODIS Daily Spectral Indices Based on in situ Measurements. Remote Sens. Basel 2011, 3, 1644–1662. [Google Scholar] [CrossRef] [Green Version]
  20. Zhou, J.; Jia, L.; Menenti, M.; Gorte, B. On the performance of remote sensing time series reconstruction methods—A spatial comparison. Remote Sens. Environ. 2016, 187, 367–384. [Google Scholar] [CrossRef]
  21. Nagai, S.; Saigusa, N.; Muraoka, H.; Nasahara, K.N. What makes the satellite-based EVI–GPP relationship unclear in a deciduous broad-leaved forest? Ecol. Res. 2010, 25, 359–365. [Google Scholar] [CrossRef]
  22. Sakamoto, T.; Yokozawa, M.; Toritani, H.; Shibayama, M.; Ishitsuka, N.; Ohno, H. A crop phenology detection method using time-series MODIS data. Remote Sens. Environ. 2005, 96, 366–374. [Google Scholar] [CrossRef]
  23. Peng, Y.; Gitelson, A.A.; Sakamoto, T. Remote estimation of gross primary productivity in crops using MODIS 250m data. Remote Sens. Environ. 2013, 128, 186–196. [Google Scholar] [CrossRef]
  24. Housman, I.; Chastain, R.; Finco, M. An Evaluation of Forest Health Insect and Disease Survey Data and Satellite-Based Remote Sensing Forest Change Detection Methods: Case Studies in the United States. Remote Sens. Basel 2018, 10, 1184. [Google Scholar] [CrossRef] [Green Version]
  25. Van Leeuwen, W.J.; Huete, A.R.; Laing, T.W. MODIS vegetation index compositing approach: A prototype with AVHRR data. Remote Sens. Environ. 1999, 69, 264–280. [Google Scholar] [CrossRef]
  26. Vermote, E.F.; Roger, J.C.; Ray, J.P. MODIS Surface Reflectance User’s Guide, Collection 6th ed. 2015. Available online: https://modis-land.gsfc.nasa.gov/pdf/MOD09_UserGuide_v1.4.pdf (accessed on 20 March 2021).
  27. Gao, X.; Huete, A.R.; Didan, K. Multisensor comparisons and validation of MODIS vegetation indices at the semiarid Jornada experimental range. IEEE Trans. Geosci. Remote. 2003, 41, 2368–2381. [Google Scholar]
  28. Miura, T.; Didan, K.; Huete, A.R.; Rodriguez, E.P. A performance evaluation of the MODIS vegetation index compositing algorithm. In Proceedings of the IEEE 2001 International Geoscience and Remote Sensing Symposium, Sydney, NSW, Australia, 9–13 July 2001; pp. 1812–1814. [Google Scholar]
  29. Miura, T.; Huete, A.R.; Didan, K.; van Leeuwen, W.J.; Yoshioka, H. An assessment of the MODIS vegetation index compositing algorithm using quality assurance flags and sun/view angles. In Proceedings of the IEEE 2000 International Geoscience and Remote Sensing Symposium, Honolulu, HI, USA, 24–28 July 2000; pp. 545–547. [Google Scholar]
  30. Wessels, K.J.; Bachoo, A.; Archibald, S. Influence of composite period and date of observation on phenological metrics extracted from MODIS data. In Proceedings of the 33rd International Symposium on Remote Sensing of Environment: Sustaining the Millennium Development Goals, Stresa, Lago Magglore, Italy, 4–8 May 2009. [Google Scholar]
  31. McKellip, R.; Ryan, R.E.; Blonski, S.; Prados, D. Crop Surveillance Demonstration Using a Near-Daily MODIS Derived Vegetation Index Time Series. In Proceedings of the Third International Workshop on the Analysis of Multitemporal Remote Sensing Images, Biloxi, MS, USA, 16 May 2005; NASA Stennis Space Center: Biloxi, MS, USA. [Google Scholar]
  32. Guindin-Garcia, N.; Gitelson, A.A.; Arkebauer, T.J.; Shanahan, J.; Weiss, A. An evaluation of MODIS 8-and 16-day composite products for monitoring maize green leaf area index. Agric. For. Meteorol. 2012, 161, 15–25. [Google Scholar] [CrossRef] [Green Version]
  33. Colditz, R.R.; Ressl, R.A. The impact of the day of observation of image composites on adequate time series generation. In Proceedings of the SPIE Earth Resources and Environmental Remote Sensing/GIS Applications IV, Dresden, Germany, 24 October 2013. [Google Scholar]
  34. Testa, S.; Mondino EC, B.; Pedroli, C. Correcting MODIS 16-day composite NDVI time-series with actual acquisition dates. Eur. J. Remote Sens. 2014, 47, 285–305. [Google Scholar] [CrossRef]
  35. Thayn, J.B.; Price, K.P. Julian dates and introduced temporal error in remote sensing vegetation phenology studies. Int. J. Remote Sens. 2008, 29, 6045–6049. [Google Scholar] [CrossRef]
  36. Gatis, N.; Anderson, K.; Clement, E.G.; Luscombe, D.J.; Hartley, I.P.; Smith, D.; Brazier, R.E. Evaluating MODIS vegetation products using digital images for quantifying local peatland CO2 gas fluxes. Remote Sens. Ecol. Conserv. 2017, 3, 217–231. [Google Scholar] [CrossRef]
  37. Fisher, J.I.; Mustard, J.F.; Vadeboncoeur, M.A. Green leaf phenology at Landsat resolution: Scaling from the field to the satellite. Remote Sens. Environ. 2006, 100, 265–279. [Google Scholar] [CrossRef]
  38. Wang, C.; Zhu, K. Misestimation of Growing Season Length Due to Inaccurate Construction of Satellite Vegetation Index Time Series. IEEE Geosci. Remote. Sens. Lett. 2019, 16, 1185–1189. [Google Scholar] [CrossRef]
  39. Ahl, D.E.; Gower, S.T.; Burrows, S.N.; Shabanov, N.V.; Myneni, R.B.; Knyazikhin, Y. Monitoring spring canopy phenology of a deciduous broadleaf forest using MODIS. Remote Sens. Environ. 2006, 104, 88–95. [Google Scholar] [CrossRef] [Green Version]
  40. Testa, S.; Soudani, K.; Boschetti, L.; Borgogno Mondino, E. MODIS-derived EVI, NDVI and WDRVI time series to estimate phenological metrics in French deciduous forests. Int. J. Appl. Earth Obs. 2018, 64, 132–144. [Google Scholar] [CrossRef]
  41. Xu, X.; Conrad, C.; Doktor, D. Optimising Phenological Metrics Extraction for Different Crop Types in Germany Using the Moderate Resolution Imaging Spectrometer (MODIS). Remote Sens. Basel 2017, 9, 254. [Google Scholar] [CrossRef] [Green Version]
  42. Qiu, B.; Feng, M.; Tang, Z. A simple smoother based on continuous wavelet transform: Comparative evaluation based on the fidelity, smoothness and efficiency in phenological estimation. Int. J. Appl. Earth Obs. 2016, 47, 91–101. [Google Scholar] [CrossRef]
  43. Chen, J.; Jönsson, P.; Tamura, M.; Gu, Z.; Matsushita, B.; Eklundh, L. A simple method for reconstructing a high-quality NDVI time-series data set based on the Savitzky–Golay filter. Remote Sens. Environ. 2004, 91, 332–344. [Google Scholar] [CrossRef]
  44. Pan, Z.; Hu, Y.; Cao, B. Construction of smooth daily remote sensing time series data: A higher spatiotemporal resolution perspective. Open Geospat. Data Softw. Stand. 2017, 2, 25. [Google Scholar] [CrossRef] [Green Version]
  45. Jönsson, P.; Eklundh, L. TIMESAT—A program for analyzing time-series of satellite sensor data. Comput. Geosci. 2004, 30, 833–845. [Google Scholar] [CrossRef] [Green Version]
  46. Van Leeuwen, W.J.; Huete, A.R.; Laing, T.W.; Didan, K. Vegetation change monitoring with spectral indices: The importance of view and sun angle standardized data. In Proceedings of the SPIE The International Society for Optical Engineering, Florence, Italy, 17 December 1999. [Google Scholar]
  47. Bhandari, S.; Phinn, S.; Gill, T. Assessing viewing and illumination geometry effects on the MODIS vegetation index (MOD13Q1) time series: Implications for monitoring phenology and disturbances in forest communities in Queensland, Australia. Int. J. Remote Sens. 2011, 32, 7513–7538. [Google Scholar] [CrossRef]
  48. Tittebrand, A.; Spank, U.; Bernhofer, C.H. Comparison of satellite-and ground-based NDVI above different land-use types. Theor. Appl. Climatol. 2009, 98, 171–186. [Google Scholar] [CrossRef]
  49. Disney, M.; Lewis, P.; Thackrah, G.; Quaife, T.; Barnsley, M. Comparison of MODIS broadband albedo over an agricultural site with ground measurements and values derived from Earth observation data at a range of spatial scales. Int. J. Remote Sens. 2004, 25, 5297–5317. [Google Scholar] [CrossRef]
  50. Gu, J.; Li, X.; Huang, C.; Okin, G.S. A simplified data assimilation method for reconstructing time-series MODIS NDVI data. Adv. Space Res. 2009, 44, 501–509. [Google Scholar] [CrossRef]
  51. Zhao, W.; Duan, S. Reconstruction of daytime land surface temperatures under cloud-covered conditions using integrated MODIS/Terra land products and MSG geostationary satellite data. Remote Sens. Environ. 2020, 247, 111931. [Google Scholar] [CrossRef]
  52. Gao, Y.; Xie, H.; Lu, N.; Yao, T.; Liang, T. Toward advanced daily cloud-free snow cover and snow water equivalent products from Terra–Aqua MODIS and Aqua AMSR-E measurements. J. Hydrol. 2010, 385, 23–35. [Google Scholar] [CrossRef]
  53. Van Leeuwen, W.J.D.; Orr, B.J.; Marsh, S.E.; Herrmann, S.M. Multi-sensor NDVI data continuity: Uncertainties and implications for vegetation monitoring applications. Remote Sens. Environ. 2006, 100, 67–81. [Google Scholar] [CrossRef]
  54. Luo, Y.; Trishchenko, A.P.; Khlopenkov, K.V. Developing clear-sky, cloud and cloud shadow mask for producing clear-sky composites at 250-m spatial resolution for the seven MODIS land bands over Canada and North America. Remote Sens. Environ. 2008, 112, 4167–4185. [Google Scholar] [CrossRef]
  55. Wilson, A.M.; Parmentier, B.; Jetz, W. Systematic land cover bias in Collection 5 MODIS cloud mask and derived products—A global overview. Remote Sens. Environ. 2014, 141, 149–154. [Google Scholar] [CrossRef]
  56. Wang, T.; Fetzer, E.J.; Wong, S.; Kahn, B.H.; Yue, Q. Validation of MODIS cloud mask and multilayer flag using CloudSat-CALIPSO cloud profiles and a cross-reference of their cloud classifications. J. Geophys. Res. Atmos. 2016, 121, 611–620. [Google Scholar] [CrossRef]
  57. Sun, L.; Gao, F.; Anderson, M.C.; Kustas, W.P.; Thenkabail, P.S. Remote Sensing Daily Mapping of 30 m LAI and NDVI for Grape Yield Prediction in California Vineyards. Remote Sens. Basel 2017, 9, 317. [Google Scholar] [CrossRef] [Green Version]
  58. Rumney, G.R. North America, climate of. In Climatology. Encyclopedia of Earth Science; Springer: Boston, MA, USA, 1987; pp. 612–624. [Google Scholar] [CrossRef]
  59. Roerink, G.J.; Menenti, M.; Verhoef, W. Reconstructing cloudfree NDVI composites using Fourier analysis of time series. Int. J. Remote Sens. 2000, 21, 1911–1917. [Google Scholar] [CrossRef]
  60. Xin, Q.; Li, J.; Li, Z.; Li, Y.; Zhou, X. Evaluations and comparisons of rule-based and machine-learning-based methods to retrieve satellite-based vegetation phenology using MODIS and USA National Phenology Network data. Int. J. Appl. Earth Obs. 2020, 93, 102189. [Google Scholar] [CrossRef]
  61. Eklundh, L.; Jönsson, P. TIMESAT 3.3 with seasonal trend decomposition and parallel processing Software Manual. 2017. Available online: http://web.nateko.lu.se/timesat/docs/TIMESAT33_SoftwareManual.pdf (accessed on 20 March 2021).
Figure 1. The location of CSP with three agriculture experimental sites (the satellite image was from CSP website: http://csp.unl.edu/public/, accessed on 20 March 2021).
Figure 1. The location of CSP with three agriculture experimental sites (the satellite image was from CSP website: http://csp.unl.edu/public/, accessed on 20 March 2021).
Remotesensing 13 01397 g001
Figure 2. The distribution of the selected forty-six evenly distributed AmeriFlux (AF) sites.
Figure 2. The distribution of the selected forty-six evenly distributed AmeriFlux (AF) sites.
Remotesensing 13 01397 g002
Figure 3. The flowchart of the DAVIR-MUTCOP method.
Figure 3. The flowchart of the DAVIR-MUTCOP method.
Remotesensing 13 01397 g003
Figure 4. The temporal variation of original daily, 8-day and 16-day NDVI data derived from MOD09GQ, MOD09Q1 and MOD13Q1 products, respectively (named 09GQ, 9Q1O and 13Q1O, respectively); temporally corrected 8-day and 16-day NDVI data using actual acquisition dates (named 9Q1A and 13Q1A, respectively); daily clear NDVI data filter by MOD35_L2 cloud mask data (named 35L2); and the clear NDVI values selected by the DAVIR-MUTCOP8 and the DAVIR-MUTCOP16 method (named DM8 and DM16) at the second CSP site (ac) and the 41th AF site (df).
Figure 4. The temporal variation of original daily, 8-day and 16-day NDVI data derived from MOD09GQ, MOD09Q1 and MOD13Q1 products, respectively (named 09GQ, 9Q1O and 13Q1O, respectively); temporally corrected 8-day and 16-day NDVI data using actual acquisition dates (named 9Q1A and 13Q1A, respectively); daily clear NDVI data filter by MOD35_L2 cloud mask data (named 35L2); and the clear NDVI values selected by the DAVIR-MUTCOP8 and the DAVIR-MUTCOP16 method (named DM8 and DM16) at the second CSP site (ac) and the 41th AF site (df).
Remotesensing 13 01397 g004
Figure 5. The temporal variation of the NDVI time series curve from 2001 to 2012 and the available averaged clear observations per month in the three CSP sites for NDVI time-series reconstruction by the SG8A, SG16A, DAVIR-MUTCOP8 (DM8) and DAVIR-MUTCOP16 (DM16) methods. The averaged clear observations per month (0 ~ 30) from January to December during the period from 2001 to 2012 in three CSP sites was presented in different colors which varied from blue to red with corresponding values (averaged clear observations per month) recorded on the color pixels.
Figure 5. The temporal variation of the NDVI time series curve from 2001 to 2012 and the available averaged clear observations per month in the three CSP sites for NDVI time-series reconstruction by the SG8A, SG16A, DAVIR-MUTCOP8 (DM8) and DAVIR-MUTCOP16 (DM16) methods. The averaged clear observations per month (0 ~ 30) from January to December during the period from 2001 to 2012 in three CSP sites was presented in different colors which varied from blue to red with corresponding values (averaged clear observations per month) recorded on the color pixels.
Remotesensing 13 01397 g005
Figure 6. The temporal variation of NDVI time series curve from 2001 to 2012 and the averaged clear observations per month for NDVI time-series reconstruction by the SG8A, SG16A, DAVIR-MUTCOP8 (DM8) and DAVIR-MUTCOP16 (DM16) methods. The averaged clear observations per month (0~30) from January to December during the period from 2006 to 2010 in six selected AF sites was presented in different colors varying from blue to red with corresponding values (averaged clear observations per month) recorded on the color pixels. The six selected AF sites under different land cover types include P10 (a), P15 (b), P17 (c), P24 (d), P40 (e) and P41 (f).
Figure 6. The temporal variation of NDVI time series curve from 2001 to 2012 and the averaged clear observations per month for NDVI time-series reconstruction by the SG8A, SG16A, DAVIR-MUTCOP8 (DM8) and DAVIR-MUTCOP16 (DM16) methods. The averaged clear observations per month (0~30) from January to December during the period from 2006 to 2010 in six selected AF sites was presented in different colors varying from blue to red with corresponding values (averaged clear observations per month) recorded on the color pixels. The six selected AF sites under different land cover types include P10 (a), P15 (b), P17 (c), P24 (d), P40 (e) and P41 (f).
Remotesensing 13 01397 g006
Figure 7. The NDVI time-series curves of the 41th AF sites (P41) from 2006 to 2010 reconstructed by SG35L2 (a), SG8A/16A and DAVIR-MUTCOP8/16 (DM8/16) methods (b) using different Fsc or/and Fse in the SG filters.
Figure 7. The NDVI time-series curves of the 41th AF sites (P41) from 2006 to 2010 reconstructed by SG35L2 (a), SG8A/16A and DAVIR-MUTCOP8/16 (DM8/16) methods (b) using different Fsc or/and Fse in the SG filters.
Remotesensing 13 01397 g007
Figure 8. The scatter plot of ground measured LAI against NDVI values reconstructed DAVIR-MUTCOP8/16 method (a), SG8A/16A (b) and SG35L2 methods (c).
Figure 8. The scatter plot of ground measured LAI against NDVI values reconstructed DAVIR-MUTCOP8/16 method (a), SG8A/16A (b) and SG35L2 methods (c).
Remotesensing 13 01397 g008
Figure 9. The scatter plot of NDVI values reconstructed by DAVIR-MUTCOP8 against those by the DAVIR-MUTCOP8 method (a) and the scatter plot of NDVI values reconstructed by SG8A against those by the SG16A method (b).
Figure 9. The scatter plot of NDVI values reconstructed by DAVIR-MUTCOP8 against those by the DAVIR-MUTCOP8 method (a) and the scatter plot of NDVI values reconstructed by SG8A against those by the SG16A method (b).
Remotesensing 13 01397 g009
Figure 10. The NDVI time-series curves as well as the scatters plots of the NDVI values of the 39th AF sites (P39) from 2006 to 2010 reconstructed by SG8A (a,e), SG16A (b,f), DAVIR-MUTCOP8 (c,g) and DAVIR-MUTCOP16 (d,h) methods using different Fsc.
Figure 10. The NDVI time-series curves as well as the scatters plots of the NDVI values of the 39th AF sites (P39) from 2006 to 2010 reconstructed by SG8A (a,e), SG16A (b,f), DAVIR-MUTCOP8 (c,g) and DAVIR-MUTCOP16 (d,h) methods using different Fsc.
Remotesensing 13 01397 g010
Figure 11. The NDVI time-series curves of the 10th AF sites (P10) from 2006 to 2010 reconstructed by SG8A (a), SG16A (b), DAVIR-MUTCOP8 (c) and DAVIR-MUTCOP16 (d) methods using different Fsc.
Figure 11. The NDVI time-series curves of the 10th AF sites (P10) from 2006 to 2010 reconstructed by SG8A (a), SG16A (b), DAVIR-MUTCOP8 (c) and DAVIR-MUTCOP16 (d) methods using different Fsc.
Remotesensing 13 01397 g011
Figure 12. The scatter plot of ground-observed LAI against NDVI values reconstructed by DAVIR-MUTCOP8 (a), DAVIR-MUTCOP16 (b) and SG35L2 (c) methods with Fse varying from 5 to 15.
Figure 12. The scatter plot of ground-observed LAI against NDVI values reconstructed by DAVIR-MUTCOP8 (a), DAVIR-MUTCOP16 (b) and SG35L2 (c) methods with Fse varying from 5 to 15.
Remotesensing 13 01397 g012
Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Share and Cite

MDPI and ACS Style

Zeng, L.; Wardlow, B.D.; Hu, S.; Zhang, X.; Zhou, G.; Peng, G.; Xiang, D.; Wang, R.; Meng, R.; Wu, W. A Novel Strategy to Reconstruct NDVI Time-Series with High Temporal Resolution from MODIS Multi-Temporal Composite Products. Remote Sens. 2021, 13, 1397. https://0-doi-org.brum.beds.ac.uk/10.3390/rs13071397

AMA Style

Zeng L, Wardlow BD, Hu S, Zhang X, Zhou G, Peng G, Xiang D, Wang R, Meng R, Wu W. A Novel Strategy to Reconstruct NDVI Time-Series with High Temporal Resolution from MODIS Multi-Temporal Composite Products. Remote Sensing. 2021; 13(7):1397. https://0-doi-org.brum.beds.ac.uk/10.3390/rs13071397

Chicago/Turabian Style

Zeng, Linglin, Brian D. Wardlow, Shun Hu, Xiang Zhang, Guoqing Zhou, Guozhang Peng, Daxiang Xiang, Rui Wang, Ran Meng, and Weixiong Wu. 2021. "A Novel Strategy to Reconstruct NDVI Time-Series with High Temporal Resolution from MODIS Multi-Temporal Composite Products" Remote Sensing 13, no. 7: 1397. https://0-doi-org.brum.beds.ac.uk/10.3390/rs13071397

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop