A comparison of methods for smoothing and gap filling time series of remote sensing observations – application to MODIS LAI products

Moderate resolution satellite sensors including MODIS (Moderate Resolution Imaging Spectroradiometer) already provide more than 10 yr of observations well suited to describe and understand the dynamics of earth’s surface. However, these time series are associated with significant uncertainties and incomplete because of cloud cover. This study compares eight methods designed to improve the continuity by filling gaps and consistency by smoothing the time course. It includes methods exploiting the time series as a whole (iterative caterpillar singular spectrum analysis (ICSSA), empirical mode decomposition (EMD), low pass filtering (LPF) and Whittaker smoother (Whit)) as well as methods working on limited temporal windows of a few weeks to few months (adaptive Savitzky–Golay filter (SGF), temporal smoothing and gap filling (TSGF), and asymmetric Gaussian function (AGF)), in addition to the simple climatological LAI yearly profile (Clim). Methods were applied to the MODIS leaf area index product for the period 2000–2008 and over 25 sites showed a large range of seasonal patterns. Performances were discussed with emphasis on the balance achieved by each method between accuracy and roughness depending on the fraction of missing observations and the length of the gaps. Results demonstrate that the EMD, LPF and AGF methods were failing because of a significant fraction of gaps (more than 20 %), while ICSSA, Whit and SGF were always providing estimates for dates with missing data. TSGF (Clim) was able to fill more than 50 % of the gaps for sites with more than 60 % (80 %) fraction of gaps. However, investigation of the accuracy of the reconstructed values shows that it degrades rapidly for sites with more than 20 % missing data, particularly for ICSSA, Whit and SGF. In these conditions, TSGF provides the best performances that are significantly better than the simple Clim for gaps shorter than about 100 days. The roughness of the reconstructed temporal profiles shows large diff rences between th various methods, with a decrease of the roughness with the fraction of missing data, except for ICSSA. TSGF provides the smoothest temporal profiles for sites with a % gap > 30 %. Conversely, ICSSA, LPF, Whit, AGF and Clim provide smoother profiles than TSGF for sites with a % gap < 30 %. Impact of the accuracy and smoothness of the reconstructed time series were evaluated on the timing of phenological stages. The dates of start, maximum and end of the season are estimated with an accuracy of about 10 days for the sites with a % gap < 10 % and increases rapidly with the % gap. TSGF provides more accurate estimates of phenologi al timing up to a % gap < 60 %.


Introduction
Leaf area index (LAI) is recognized as an essential climate variable (GCOS, 2006) since it plays a key role in vegetation functioning and exchanges of mass and energy between the atmosphere, the plant and the soil.LAI is defined as half the area of the green elements per unit horizontal surface (Chen and Black, 1992).Satellite observations in the reflective solar domain have been used intensively for more than a decade to monitor LAI dynamics over the globe using medium resolution sensors such as MODIS (Moderate Resolution Imaging Spectroradiometer) (Myneni et al., 2002), VEGETATION (Deng et al., 2006;Baret et al., 2007), MERIS (MEdium Resolution Imaging Spectrometer) (Bacour et al., 2006a) or AVHRR (Advanced Very High Resolution Radiometer) Published by Copernicus Publications on behalf of the European Geosciences Union.
The current products show significant discontinuities mainly due to cloud or snow cover (Weiss et al., 2007) as well as system failure.Further, because of residual cloud contamination, imperfect atmospheric or directional corrections as well as instability of the inversion process, products may be characterized by significant temporal noise not expected with the actual LAI values.LAI is the result of incremental processes in vegetation such as leaf development and senescence.Therefore, LAI should be relatively smooth with time.However, this temporal smoothness property is not always observed over current LAI products as demonstrated by several authors (Weiss et al., 2007;Verger et al., 2008;Camacho et al., 2013).The MODIS LAI product was recognized as less smooth mainly because of a combination of factors including a short (8 days) compositing window, the maximum value compositing algorithm used and instability in the retrieval algorithm, though MODIS collection 5 products are an improvement over the collection 4 products (De Kauwe et al., 2011).Several investigations have been focusing on the improvement of the MODIS products using specific mathematical filters which use either temporal or spatial techniques to get temporally smoothed and spatially continuous products (Gao et al., 2008;Borak and Jasinski, 2009;Jiang et al., 2010;Verger et al., 2011;Yuan et al., 2011).Spatial filters using pixel-level or regional ecosystem statistical data include geostatistical and regression methods (Goovaerts, 1997;Berterretche et al., 2005;Wang et al., 2012).Nevertheless, spatial filters may fail for LAI products derived from coarse resolution satellites to represent the complexity of real landscapes mainly over mixed pixels where LAI could vary widely within a short distance.To overcome this limitation some studies tried to combine both temporal and spatial methods by using historical high-quality data and temporal curves from neighbor pixels (e.g.Moody et al., 2005;Fang et al., 2008;Gao et al., 2008).Our study refers only to temporal methods.
Time series processing is thus an important ingredient of a biophysical algorithm in order to get the expected continuous and smooth dynamics required by many applications.The earliest methods used in remote sensing, often called compositing, were reviewed by Qi and Kerr (1997) They mostly operate over a local temporal window, focusing on minimizing artifacts due to cloud or snow contamination, atmospheric or directional residual effects.They included the well-known MVC (maximum value compositing) (Holben, 1986) and BISE (Best Index Slope Extraction from Viovy et al., 1992).More recently, Hird and McDermid (2009a) reviewed the abundant literature on time series processing in remote sensing, mostly focusing on NDVI (normalized difference vegetation index) (Rouse et al., 1974).They further compared six methods operating over a restricted temporal window and demonstrated that under their conditions, logistic (Jönsson and Eklundh, 2004) or asymmetric Gaussian function (Beck et al., 2006) curve fitting methods were outperforming more simple local filtering methods (Jiang et al., 2010), compared three statistical methods both to smooth the time series and to provide forecast over a season.These methods are based on the decomposition of the time series into noise, seasonal variability and trend and require, therefore, relatively long time series of observations.Fourier transforms (Azzali and Menenti, 1999), or wavelet decomposition (Martinez and Gilabert, 2009) have also been used to characterize the phenology of vegetation from medium resolution observations.However, several studies have demonstrated the superiority of local methods, i.e. based on a restricted temporal window, as compared to the Fourier transform methods applied to the whole time series (Jönsson and Eklundh, 2002;Beck et al., 2006;Ma and Veroustraete, 2006;Hird and McDermid, 2009a).Physically based corrections were also proposed in order to correct for the known factors of variability, resulting in the GIMMS data set (Tucker et al., 2005).However orbital drift and directionality were rather corrected using Empirical Mode Decomposition techniques (Pinzon et al., 2005).More recently, the long-term data record (LTDR) series derived from AVHRR sensors (Vermote et al., 2009) and CYCLOPES (Baret et al., 2007) and GEOV1 (Meroni et al., 2013) derived from VEG-ETATION also proposed global time series based on physical principles.Alcaraz-Segura et al. (2010) and Meroni et al. (2013) showed that significant differences were observed between these several NDVI time series, making the identification of anomalies and trends more complex.The choice of the smoothing gap filling or compositing method may have a large impact on the accuracy of the phenology extracted from the reconstructed time series (Hird and McDermid, 2009b;Atkinson et al., 2012).
This brief review of studies focusing on satellite time series from medium resolution sensors shows that a number of methods are available.It is however still difficult to identify the potentials and limitations associated since no comprehensive evaluation is available.Comparison is often qualitative, or when quantitative, it is mostly centered on a small sample of global conditions.Most of them have been applied to NDVI rather than on true biophysical variable such as LAI.Further, very little attention was paid to the missing data structure: as a matter of fact satellite observations present missing data mostly because of cloud masking which creates irregular time steps between actual observations of the surface.Gap filling is therefore an important aspect of the process.Finally, only a small fraction of the studies were employing methods capable of processing the time series as a whole such as the decomposition methods.
The objective of this study is to evaluate the capacity of several methods to provide faithful reconstruction of time series in the presence of significant amount of missing observations as well as observations contaminated by uncertainties.The methods will therefore be compared using several criterions including the ability to run over periods without observations of variable length, the fidelity of reconstructed values with the actual ones and the smoothness of the reconstructed temporal profiles.In addition, consequences on the ability to capture phenological stages will be also quantified since this is a usual application of the time series Eight methods were selected because they were well referenced while being based either on local curve fitting techniques, or decomposition techniques working on the time series as a whole.Except ICSSA and EMD, all other methods are commonly used for processing biophysical time series data.ICSSA and EMD, commonly used in other subject areas, were considered in this study due to their potential ability to recover seasonal trends and reduce noise.The study is based on the MODIS LAI collection 5 product (Shabanov et al., 2005) at 1 km spatial and 8 days temporal sampling over a 9 yr period.The MODIS LAI products were demonstrated to get relatively good accuracy (closeness to actual ground observations) but were suffering from lack of precision (temporal or spatial consistency), resulting in shaky temporal profiles as stated earlier.Processing of such time series is therefore expected to result in a significant improvement of its consistency.A sample of sites selected to represent the range of variability expected over the globe was considered.

Approach, data and methods
The MODIS LAI products will first be described.Then the 8 methods selected will be briefly presented.Finally, the approach for evaluating the methods and the associated metrics used will be described.

The MODIS data and preprocessing
The data used in this study are the MODIS Collection 5 LAI products (MOD15A2) derived from TERRA and AQUA platforms.The products were downloaded from the land processes distributed active archive center (https://lpdaac.usgs.gov) for the 2000-2008 9 yr period.They correspond to 1 km spatial sampling interval using a sinusoidal projection system.The temporal sampling is 8 days based on a daily composition: all observations available in the 8-day compositing window are accumulated, and the one getting the maximum FaPAR (Fraction of Absorbed Photosynthetically Active Radiation) value is selected.The main MODIS LAI retrieval algorithm relies on the inversion of a 3-D radiative transfer model using the red and near-infrared bidirectional reflectance factor values, associated uncertainties, the viewillumination geometry and biome type (within eight types based on MOD12Q1 land-cover map) as inputs (Myneni et al., 2002;Shabanov et al., 2005).If the main algorithm fails, a back-up procedure is triggered to estimate LAI from biome specific NDVI-based relationships.However, the LAI estimates using the back-up algorithm are of lower quality mostly due to residual clouds and poor atmospheric correction (Yang et al., 2006a, b).Hence, these estimates are not used in this study and are considered as missing observations.Although the MODIS LAI product has been extensively validated (e.g.De Kauwe et al., 2011;Ganguly et al., 2008), a high level of noise was inducing shaky temporal profiles and unrealistic seasonality (Kobayashi et al., 2010), which justifies the interest of using MODIS LAI products for smoothing and gap-filling investigations.
A first preprocessing step was applied to remove unexpected abrupt variations in the time series: values that are substantially different from both their left-and right-hand neighbors and from the median in a 72 day length local window are considered as missing values as proposed (Jönsson and Eklundh, 2004) in the TIMESAT toolbox.Further, for Evergreen broadleaf forests presenting reduced seasonality and high level of variability in the time series because of frequent occurrence of residual cloud contamination, any value lower than the first decile are eliminated since these usually low LAI values are not expected in this canopy type.

The methods investigated
Eight candidate methods (Table 1) were selected, including both decomposition techniques generally applied to the whole time series and curve fitting methods working on a limited temporal window.Decomposition methods split the signal into additive components.The time series are then reconstructed using only the components of interest, usually removing the high frequency components considered as noise.Decomposition methods should capture the seasonality and the trend signals observed over the whole time series, which may be exploited in the reconstruction phase to replace missing values.Curve fitting techniques adjust the parameters of a functional by minimizing a cost function that is usually the sum of quadratic differences between observations and simulations.Because the adjustment is operated over a limited temporal window, only a limited amount of information is used when filling gaps.
Iterative Caterpillar Singular Spectrum Analysis Method (ICSSA) is a modification of the CSSA (Golyandina and Osipov, 2007) method developed to describe time series and fill missing data by decomposing the time series into empirical orthogonal functions (EOF).This modified version was proposed by Kandasamy et al. (2012) to correct for overestimation of seasonal valleys and better fitting to the peaks as compared to the original CSSA formulation.The method  2011) for each 8-day period * The maximum gap length applies on the yearly reconstructed time series, not on the original time series.requires 2 parameters: the window length and the number of eigenvectors (orthogonal functions) used for the reconstruction.Better reconstruction can be obtained for large number of eigenvectors but at the cost of a decrease in the smoothness.After trial and error, the number of eigenvector was set to 1, and the window length was set to 40 days.This method allows filling gaps and forecasting data at the extremities of the time series.
Empirical mode decomposition method (EMD).This method proposed by Huang et al. (1998) consists in decomposing the time series into a small number of intrinsic mode frequencies (IMFs) derived directly from the time series itself using an adaptive iterative process where the data are represented by intrinsic mode functions, to which the Hilbert transform can be applied.The method requires setting 2 parameters: the threshold for convergence and the maximum number of IMFs.The threshold for convergence is set to 0.3 according to Huang et al. (1998) and the maximum number of IMFs was restricted to 10 after trial and error.The first IMF, mostly affected by noise, was smoothed using a uniform mean kernel to remove the high frequency fluctuations at the expense of a loss in the amplitude (Demir and Erturk, 2008).Note that the EMD method requires the time series to be continuous.To allow the application of EMD to MODIS time series, the missing data within 128 days were filled by linear interpolation as proposed by Verger et al. (2011).However, when the time series contains gaps longer than 128 days, the whole series was not reconstructed.As a matter of fact, linear interpolation provides generally poor performances in case of long periods without observations.
Low pass filtering (LPF).This method originally proposed by Thoning et al. (1989) was adapted by Bacour et al. (2006b) for better retrieving the seasonality from AVHRR time series.A time dependent function with 10 terms (2 polynomial and 8 harmonic terms) is first adjusted to the data.Then, the residuals of this first fitting are filtered with a low pass filter using two cut-off frequencies defined to separate the intra-annual and inter-annual variations.The final reconstruction is obtained by summing the polynomial and harmonic terms with the filtered residuals.Although this method may also be considered to be based on curve fitting, it applies over the time series as a whole.This method requires the data to be continuous.Hence, similarly to the EMD method, the gaps within 128 days were filled by linear interpolation.For gaps longer than 128 days, LPF was considered unsuccessful and, in this case, it results in missing reconstructed values for the whole time series.
Whittaker smoother (Whit).This method proposed by Whittaker (1923) is based on the minimization of a cost function describing the balance between fidelity expressed as the quadratic difference between estimates and actual observations and roughness expressed as the quadratic difference between successive estimates.This balance is controlled by a smoothing parameter.The higher this value is, the smoother the result but at the expense of fidelity.Finding an appropriate value of the smoothing parameter is difficult, as it depends on the data.After trial and error the smoothing parameter was set to 100.The smoothness is also controlled by the order of differentiation which is fixed to 3 as proposed by Eilers (2003) for time series data sampled at regular intervals but with missing observations.Adaptive Savitzky-Golay filter (SGF).This method proposed by Chen et al. (2004) iteratively applies the Savitzky-Golay filter (Savitzky and Golay, 1964) to match the upper envelope of the time series.This specific adaptation was designed to minimize the effects of cloud and snow contamination that generally decreases the estimates of vegetation indices such as NDVI as well as biophysical variables such as LAI.The original Savitzky-Golay filter consists in a local polynomial fitting with two parameters: the length of the temporal window used and the order of the polynomial.As proposed by Chen et al. (2004), the values of these parameters were optimized for each case to get the best match between observations and reconstructed values.However, the range of variation of the window length was restricted to 72-112 days, and that of the polynomial order to 2-4.Missing data in the time series were filled by linear interpolation independently of the size of the gaps as proposed in the original version.
Temporal smoothing and gap filling (TSGF).This method is another adaptation of the Savitzky-Golay filter where the polynomial degree was fixed to 2 but the temporal window may be asymmetric and variable in length.It was designed by Verger et al. (2011) to better handle time series with missing observations.The temporal window corresponding to a nominal date is adjusted to include at least 3 observations within a maximum 64 days period on each side of the nominal date.If less than 6 observations are available within the maximum ±64 days temporal window, the polynomial fitting is not applied.The gaps in the reconstructed data are filled by an iterative (2 iterations) linear interpolation within 128 days window.For periods with missing data longer than 128 days, the interpolation is not applied which results in missing data.The possible flattening of the reconstructed time series observed over peaks was further corrected by scaling the smoothed series to the observations in a local window (Verger et al., 2011).
Asymmetric Gaussian fitting (AGF).This method has been proposed by Jönsson andEklundh (2002, 2004) within the TIMESAT toolbox.A Gaussian function is adjusted locally over the growing and senescing parts of each season.The functions are finally merged to get a smooth transition from one season to another.This method can handle small gaps.The original TIMESAT implementation of this method included 3 conditions preventing the near-constant and noisy data from being processed.Two of these conditions (minimum seasonality in the data and maximum fraction of missing data of 25 %) were not considered here to enlarge its domain of validity in case of missing data.This thus allows more rigorous comparison with the other methods.However, the last condition was kept: the method was not applied over seasons with gaps longer than 72 days and, in these cases, reconstructions for the entire time series are not provided.
Climatology (Clim).The climatology describes the typical yearly time course.It was included within the set of methods investigated since it may provide smooth and complete time series with, however, no changes from year to year.The climatology was computed every 8 days during a year by averaging the values available over a ±12 day window across all the years of the time series.The climatology was then corrected to provide more continuous and consistent temporal patterns as proposed in Baret et al. (2011).To provide smoother values, a simple Savitzky-Golay filter was applied (Savitzky and Golay, 1964).Note that the climatology may present missing values when no observations in the 24 day temporal window centered on a given date in the year were available across, the years of the time series.In such situations, linear interpolation was used to fill gaps shorter than 128 days.Gaps longer than 128 days will result in missing data.Once the average yearly time course was computed, it was replicated across all the years considered to provide a reconstructed time series.

Evaluation approach
The approach proposed to evaluate the methods is based on two steps as sketched in Fig. 1.The first one is dedicated to the preparation of reference time series over a limited number of representative sites.As a result of this step, two time series will be output: (i) a time series with no gaps and small uncertainties considered as the reference (LAI ref ), and (ii) a time series with no gaps with realistic uncertainties (LAI comp ).In the second step, time series with variable occurrence of missing data will be simulated (LAI sim ) from the previously LAI comp time series.Each method will be applied on this LAI sim data to get the corresponding reconstructed time series (LAI rec ).Finally, the LAI rec obtained with each of the 8 methods will be evaluated based on a range of metrics describing the fidelity, the roughness of time series and the accuracy of phenological stages that can be derived from the time series.

Generation of the reference and completed
LAI time series.
In the first step, few sites were selected with the objective to show a wide variability of seasonal patterns while having a minimal number of missing data.For this purpose, the 420 BELMANIP2 sites identified by Baret et al. (2006)  The same process was applied to SB, GC and EB classes, resulting in a total of 25 sites (5 sites for each of the 5 biome classes) (Table 2, Fig. 2) The eight methods presented earlier were applied to each of the 25 sites and show very good agreement.The median across all methods is a good approximation of the expected LAI product values (Fig. 3) with 4 out of the 8 methods investigated very close (RMSE (root mean square error) lower than 0.05) to the median across all methods.The time series made with the median across the 8 methods will therefore be considered as the reference values, LAI ref , in the following.This LAI ref does not show any missing data since the gaps in the original time series were filled by the reconstructed values of the 8 methods.LAI ref constitutes a good reference with minimal uncertainties attached to the LAI values because of the temporal smoothing coming from each method and the computation of the median across the 8 methods.A second set of time series was generated to provide realistic LAI values: the LAI ori were complemented at the location of missing data by LAI ref values contaminated by a noise that was randomly drawn within the distribution of residuals (LAI ref -LAI ori ) for each site.These realistic but continuous temporal profiles with no gaps (LAI comp ) will be used in the second step of the approach for simulating time series with gaps.

Simulation of time series with gaps
In the second step, emphasis was put on the occurrence of missing data (% gap).The gap structure observed over each one of the 420 sites was applied to the completed time series (LAI comp ).This allows the gap structure to be more realistic as compared to other strategies that would consist in randomly simulating gaps.However, vegetation type and the associated climate experienced, hence the cloud occurrence and corresponding gap structure, are probably correlated.To account for such possible dependency, the gap structure applied to one of the 25 sites was selected within the gap sites belonging to the same vegetation class (Table 3).Note that the balance amongst vegetation classes in BELMANIP2 was preserved (Table 3) providing approximate representativeness of global scale conditions regarding the occurrence of missing data: SB and CG represent about 2/3 of the land area, associated with relatively low fraction of missing data (gap percentage, Fig. 4).Forests represent about 1/3 of the global land area with relatively high fraction of missing data.However, sites with less than 9 observations for the whole 9 yr period (i.e. less than one observation per year in average) were not considered since none of the methods will be able to provide a fair reconstruction of the LAI time course.2) and the 384 sites used to simulate gaps in the data (×).SB: shrub savannah bare, GC: crop grassland, DB: deciduous broadleaf forests, EB: evergreen broadleaf forests, NF: needleleaf forests.A total of 384 sites were used to simulate the gaps in the data (Table 3).Because each vegetation class is represented by 5 sites used for the LAI ref and LAI comp values, a total of 1920 cases (384 × 5) with realistic LAI uncertainties and gap structure were finally available.

Metrics used to quantify performances
The performances of the 8 methods considered in this study were evaluated based both on LAI values as well as phenology.Note that when the reconstructed LAI values (LAI rec Fig. 1) were outside the definition domain (0 < LAI rec < 7), the reconstructed value was systematically set to the closest bound (0 or 7).Note that in several situations, the methods may fail to reconstruct the whole time series due to long periods of gaps.This will be quantified by the reconstruction fraction (% reconstructions), i.e the fraction of dates with reconstructed values in LAI rec time series, and the fraction of successful gap filling (% success), i.e. the fraction of gaps that were able to be filled.

Metrics based on LAI values
The RMSE (root mean square error) computed over all cases quantifies the fidelity of the reconstruction of the time series: where LAI j ref (t) and LAI j rec (t) are, respectively, the reference and the reconstructed values for date t and case j , n j is the number of dates with observations for case j and N is the number of cases considered.
Finally, the metrics proposed by Whittaker (1923) called here roughness will be used to quantify the shaky nature of the reconstructed time series: (2)

Metrics based on phenology
The 5 evergreen broadleaf forest sites were excluded from these metrics since the identification of seasonality was questionable at the single pixel scale considered in this study, and would result in large uncertainties in the timing of phenological stages if they exist (some sites do not show obvious seasonality).Three phenological events were considered: the start of season (SoS), maximum of season (MoS) and end of season (EoS).SoS and EoS were defined similarly to Jönsson and Eklundh (2002) as the timing when LAI reaches 20 % where P j ref (s) and P j rec (s) are, respectively, the reference and the reconstructed dates for the phenological events and case j , m j is the number of phenological events for case j (i.e. the number of seasons in the time series j and N is the number of cases considered).

Results
The methods will first be evaluated with regards to fidelity and roughness of the reconstructed time series.Then, they will be evaluated with regard to their ability to describe the phenology.In both the cases, the impact of the occurrence of missing data (% gap) and the length of gaps (LoG), i.e. the number of days between two consecutive valid observations, will be analyzed.

Performances for LAI reconstruction
Before investigating quantitatively the performances through the several metrics envisioned, the main features of each method will be qualitatively assessed.Five cases of LAI rec within the 1920 ones have been selected (Fig. 5) and their temporal profiles plotted against LAI ref .They represent the 5 typical vegetation types under a medium occurrence of missing data.Visual inspection shows that -The climatology is often shifted from the reference temporal profile, highlighting the inter-annual fluctuations, particularly for non-forest vegetation types (SB and CG in Fig. 5).
-In presence of periods with long and continuous missing data, several methods were not able to reconstruct the time series over these periods, particularly TSGF and Clim, while AGF, EMD, LPF fail for the entire time series (SB in Fig. 5).However, the other methods (ICSSA, Whit, SGF) showing continuity in LAI rec do not always provide realistic (as compared with LAI ref ) reconstructions in such cases.
-When observations show a significant level of temporal noise (the forest sites in  both in terms of fidelity (closeness to LAI ref ) and roughness, particularly for SGF and EMD.

Capacity to reconstruct the temporal series in presence of missing data
All the methods were not able to fill all of the gaps, i.e. to provide an estimated value in gaps.This was quantified by the % success, i.e. the fraction of gaps that were able to be filled.Whit, SGF, ICSSA allow to fill most of the gaps even if they are very long (Fig. 6a).Conversely, EMD, LPF and AGF show a rapid decrease of %Success with the length of the gaps, with no reconstructions for gaps longer than 88 (AGF) to 136 days (EMD, LPF).Even for small gaps, only 50 % of them were filled.This is due to the fact that a specific gap may be associated to other ones in a close vicinity in the time series.TSGF is able to fill gaps up to gap length of 128 days as expected by its definition.The climatology shows also a progressive decrease of % success with gaps of 128 days length being filled in 80 % of cases because of the accumulation of observations over the 9 yr period.The capacity to fill individual gaps has consequences on the reconstruction fraction (% reconstructions), i.e. the fraction of dates with reconstructed LAI values, LAI rec , relative to the total number of dates in the LAI comp time series.Only three methods (Whit, SGF, ICSSA) were able to provide a continuous reconstructed time series over all the cases investigated (Fig. 6b) even for large occurrence of missing data in agreement with Fig. 6a.In contrast, AGF is characterized by the smallest %Success (Fig. 6a) and % reconstruction (Fig. 6b): only 50 % of the dates are reconstructed for cases with more than 25 % of missing data in their time series (Fig. 6b).LPF and EMD that do not accept gaps longer than 128 days (Table 1) show also a similar drastic decrease of the reconstruction fraction with the occurrence of missing data in the cases considered (Fig. 6b).The TSGF method, although also not filling gaps longer than 128 days is more resilient to the occurrence of gaps: TSGF was able to reconstruct more than 50 % of the data for cases with more than 60 % of missing data.When a gap longer than 128 days appears in a time series, the remaining parts of the time series are reconstructed.This was not the case for LPF and EMD for which the whole time series was not reconstructed for cases having a gap longer than 128 days.Clim allows for the reconstruction of most time series, even for cases with large amounts of missing data, benefiting from the replications between years, cloudy days being not always the same day of the year.
To improve the robustness of the metrics used to characterize the performances on LAI and phenology reconstruction they will be computed only when % reconstructions > 50 % (Fig. 6b).As a consequence, all the methods will be compared for cases with less than % gap < 20 %; TSGF, Clim, ICSSA, SGF and Whit will be compared for 20 % < % gap < 60 %; and for % gap > 60 %, only Clim, IC-SSA, SGF and Whit will be compared (Fig. 6b).

Fidelity to LAI ref
Fidelity is quantified by RMSE.To better highlight the reconstruction capacity of the methods, RMSE were computed by comparing LAI rec with LAI ref either over actual dates with observations or over dates with missing data in LAI sim .Results show that, except for Clim and SGF, all the other methods show generally good fidelity with LAI ref for the dates where observations are available (Fig. 7).These good performances are observed almost independently from the occurrence of missing data (Fig. 7).The higher RMSE values observed for SGF are due to a positive bias induced by the fitting of the upper envelope of the observations.Clim shows a RMSE value close to that of SGF (Fig. 7) that mostly refers to the inter-annual variability of LAI seasonality.Over dates of missing observations, the reconstruction capacity degrades rapidly as a function of the length of gaps for all methods except Clim that keeps a RMSE value around 0.35 independent from the expected gap length (Fig. 8a).LPF and TSGF provide the best performances up to gap length around 100 days when Clim starts to be the best method.AGF, Whit EMD and SGF show similar performances with RMSE lower than the climatology up to gap length around 70 days, while ICSSA performances rapidly degrade with the length of gaps.The fidelity of reconstructions in gaps as a function of the fraction of missing observations in the time series (Fig. 8b) derives logically from the reconstruction performances as a function of the gap length (Fig. 8a).The RMSE values computed over dates of missing observations are relatively low for all methods up to % gap < 20 %.Then SGF, Whit and ICSSA show a rapid increase of the RMSE with % gap with poorer performances as compared to Clim for % gap > 30 % (Fig. 8b).These three methods show obvious artifacts when reconstructing long gaps (Fig. 5, nonforest sites 5 and 338).TSGF shows relatively low RMSE values up to 60 % gap (Figs. 7,8b).Clim shows similar performances over dates with missing data (Fig. 8b) and dates with observations (Fig. 7) as expected since it is not dependent on the local observations.

Roughness
The roughness was computed over the whole reconstructed time series and is presented as a function of % gap (Fig. 9).Results show that for % gap < 30 %, EMD and SGF show the highest roughness values in agreement with the previous qualitative observations (Fig. 5).The behavior of SGF is con-trolled by its iterative nature that puts emphasis on fidelity (relative to the upper envelope).For EMD, the 10 modes selected were showing variable patterns and it was difficult to find a better compromise between smoothness and fidelity.ICSSA shows a roughness value close to that of Clim for % gap < 20 %.However, the roughness of ICSSA strongly increases when % gap > 20 %.This is partially due to inconsistencies observed in its temporal pattern with abrupt variations in the periods with high discontinuities in the data (Fig. 5, jumps observed between the lowest and the highest values when data are missing for non-forest sites 5 and 338).AGF and LPF show the smoothest temporal profiles however limited to cases with % gap < 20 %.Whit provides always smooth reconstructed profiles, even for large amount of missing data.This is obviously controlled by the smoothing parameter.Whit is just slightly rougher than LPF and AGF.TSGF shows a slightly higher roughness values than Whit for the % gap < 30 %, with a significant decrease when % gap increases.This is due to the linear interpolation used to fill the gaps that explains also the decrease of roughness for EMD, SGF and Clim when % gap increases.

Performances for describing the phenology
The capacity of the several methods to identify the main phenological stages (SoS, MoS and EoS) was evaluated using the dates derived from LAI ref as a reference.The performances (RMSE in days) were analyzed as a function of the occurrence of missing data (% gap).Results show a general degradation of RMSE when % gap increases for the three stages considered.
Closer inspection of performances in terms of RMSE for SoS estimates shows large differences between methods (Fig. 10a).For complete time series (% gap = 0), RMSE values are between 3 days (LPF) and 15 days (AGF), with the exception of the climatology with a RMSE around 25 days, indicating a significant inter-annual variability in the timing of SoS.EMD, TSGF, Whit and SGF have RMSE around 10 days.For discontinuous time series, Whit, SGF and IC-SSA show a continuous and steep increase of RMSE with % gap.Conversely, the RMSE values of Clim and, in a lesser way, TSGF increase moderately with % gap.Similar patterns are observed for MoS (Fig. 10b) with however smaller differences between methods for % gap < 20 % and a slightly lower rate of increase of RMSE with % gap except for Clim.The performances for EoS (Fig. 10c) appear to be very similar to what is observed for SoS (Fig. 10a).The Climatology (Clim) performs better than ICSSA, SGF and Whit for % gap > 40 % for SoS and EoS, and for % gap > 50 % for MoS.TSGF yields the smallest RMSE for % gap > 20 % for SoS, MoS and EoS with however only small differences as compared to Clim for EoS (Fig. 10a-c).

Discussion and conclusion
This study compares 8 methods designed to improve the continuity and consistency of time series by filling gaps created by missing observations and smoothing the temporal profiles to reduce local uncertainties.However, the performances of the different methods for processing time series depend on their implementation (e.g.very different results with two variants of Savitzky-Golay filter: SGF and TSGF).The selected methods were applied here as close as possible to their standard implementation including the original parameterization as proposed by their authors.When the parameters for each method were not known they were adjusted using a trial and error approach.Other techniques based on systematic cross validation could have been implemented.This was not considered here since it would lead to significant increase in computation time not compatible with cur-rent operational processing lines capabilities.The time series considered correspond to actual MODIS LAI products over a sample of sites that were selected to be both representative of the diversity of seasonal patterns and of the distribution of the missing observations.This approach was expected to improve the realism of the context of the analysis that accounts for the implicit links between the vegetation type and the distribution of missing observations.This may be critical for filling gaps or smoothing the time series.The approach allowed defining a set of reference time series used to quantify the accuracy of each of the 8 methods as a function of the fraction of missing observations.
Results clearly show that some methods including LPF, AGF and EMD were failing in about 50 % of the situations when the fraction of missing observations was larger than 20 % which represents about 60 % of the situations investigated here.This is partly due to the principles on which these methods are based, but also partly to their implementation.Consequently, great care should be taken with the implementation of such methods to improve their rate of applicability in case of significant periods with missing observations.Conversely, ICSSA, Whit and Clim methods were applicable in almost all situations while TSGF shows intermediate behavior.
For the methods resilient to periods of missing observations of significant length, their capacity to provide realistic interpolation between actual observations was challenged in cases corresponding to medium to high fraction of missing data.SGF, designed to fit the upper envelope of observations, performs poorly (large RMSE and positive Bias) over MODIS LAI time series.Better filtering principles are thus required to reject outliers possibly contaminated by residual clouds.ICSSA and Whit show unreliable interpolated values in the medium (few weeks) to large (few months) periods of missing data although these methods are adjusted over the whole time series.The TSGF method appears to provide the most reliable interpolation capacity due to its  adaptive temporal window, although limited to gaps smaller than 128 days.For longer periods without observations, the Clim method appears to be the best one provided that enough data are available over the time series of years used to build the climatology.Note that the reconstruction performances for the best methods and for gaps shorter than 100 days fulfills the GCOS criterion on LAI uncertainties (RMSE ≤ 0.5) (GCOS, 2010) although the reconstruction uncertainty is only part of the error budget Each method is based explicitly (Whit) or implicitly (the other methods) on a balance between fidelity and roughness.This is clearly demonstrated when plotting Roughness and RMSE performances for each of the 25 selected sites (Table 2, Fig. 2) for a class of occurrence of missing data (Fig. 11).For each method, all the 25 sites are approximately organized around a line passing through the origin.The slope of this line indicates the balance between fidelity and roughness.For relatively continuous time series (0 % < % gap < 15 %), TSGF, ICSSA and Whit focus more on fidelity than smoothness (Fig. 11a).Conversely, LPF, EMD and AGF are focusing more on smoothness than fidelity.SGF constitutes a particular case because the fidelity is targeting the upper envelope of the points, resulting in larger RMSE values, while roughness is also quite important as described previously.Clim provides the steepest slope, with smooth temporal profiles but a loose match with observations.Note that the slope of Clim is in between that of LAI comp and LAI ref (for which RMSE was replaced by the standard deviation between observations).For the larger occurrence of missing data (Fig. 11b), the slopes increase significantly due to an increase of RMSE mostly due to inaccurate reconstructions in the gaps, and a decrease of roughness due to more simple patterns in the observations, except for ICSSA as noticed earlier.
The slope between RMSE and roughness (Fig. 11) appears thus a good indicator of the balance between fidelity (RMSE) and roughness of each method and its associated sets of parameters.The overall performances may be described by the distance to the ideal case (RMSE = Roughness = 0) in the (roughness, RMSE) feature plane averaged over the 25 sites considered: the closer to the origin (0,0), the smoother and the better match with LAI ref (low RMSE).The behavior of each method as a function of the occurrence of missing data is well sketched in the (performances, slopes) feature plane (Fig. 12).For low amounts of missing data, all the methods provide good performances except SGF and Clim for the reasons exposed previously.When the fraction of missing data increases, each method follows a particular pattern (the black arrow in Fig. 12) with a degradation of the performances and an increase in the slope indicating more emphasis on the smoothness of the temporal profiles through a reduction of the roughness.For medium and high occurrence of missing data, TSGF provides clearly the best overall performances although restricted to gaps smaller than 128 days, followed by Whit.SGF and ICSSA show poor performances.
The consequences of the application of the several time series processing methods on their capacity to describe phenology characteristics were finally evaluated.As expected, the methods providing the best accuracy on LAI estimation were also more accurate for dating specific phenological events such as start, maximum and end of season (Fig. 13).
The effect of gaps on the derivation of time series appears as a major limitation of the accuracy of the reconstructed temporal profiles.Techniques based on the processing of the time series as a whole (ICSSA, EMD, LPF, Whit and Clim) were not demonstrated to perform systematically better than techniques based on a limited temporal window (AGF, SGF, TSGF) although they were expected to fill long gaps with the "experience" gained across the several years available in the time series.Local methods were generally more faithful but were lacking capacity to fill long gaps.Most methods were performing poorer than Clim to gaps longer than about 100 days.Future works should therefore be dedicated to develop methods where the features derived from the exploitation of the several years available in the time series including the climatology, could be injected more explicitly as a background information for improving the reliability of methods working over a limited time window, such as a season or part of it, with emphasis on the capacity to provide accurate phenological timing as proposed in Verger et al. (2013).

Fig. 1 .
Fig. 1.Flow chart describing the approach used.M i corresponds to method i within the 8 investigated.LAI rec (M i ) corresponds to the reconstructed time series based on method M i .

Fig. 3 .
Fig. 3. Comparison between the original LAI values (LAI ori ) and the median of the reconstructed values (LAI ref ) based on the 8 methods considered over the 25 selected sites (Table 2, Fig. 2) (N = 7561; R 2 = 0.90; RMSE = 0.231; Bias = −0.008).The colors of this density plot correspond to the frequency of data in each of the 0.1 × 0.1 cells LAI values.

Fig. 4 .
Fig. 4. Cumulated distribution of the fraction of missing data (% gap) in the simulated time series (LAI sim ) for each of the 5 vegetation classes (SB, CG, DB, EB, NF).

Fig. 5 .
Fig. 5. Time series (LAI rec ) reconstructed by the several methods in presence of medium occurrence of missing data (25 % < % gap < 35 %).Black dots correspond to LAI comp at the location of observations.Empty circles on the x-axis correspond to the dates of missing data.The dashed black curve corresponds to LAI ref .Note that because of the structure of missing data, EMD, LPF and AGF were not reconstructed for sites 5 and 65, as well as AGF for site 338.
Fig.5: DB, EB and NF), significant differences are observed between the methods,

Fig. 6 .
Fig. 6.(a) Fraction of gaps reconstructed (% success) as a function of the length of the gaps (LoG).(b) Fraction of dates reconstructed (% reconstructions) as a function of the % gap.The horizontal dashed line represents the 50 % threshold of % reconstructions.The several methods are represented by different colors.Some values were slightly shifted vertically to ease the reading when curves were overlapping.

Fig. 7 .
Fig. 7. RMSE as a function of % gap.The RMSE is computed between LAI ref and the reconstructed LAI rec time series over dates with actual observations in LAI sim .The values were slightly shifted vertically to ease the reading when curves were overlapping.

Fig. 8 .Fig. 9 .
Fig. 8. RMSE as a function of the length of gaps (a) and fraction of missing observations (% gap ) (b).The RMSE is computed between LAI ref and the reconstructed LAI rec time series over dates with missing observations in LAI sim .

Fig. 10 .
Fig. 10.RMSE relative to the timing (in days) of the start of season (a), maximum of season (b) and end of season (c).The RMSE is evaluated between the phenological dates computed with LAI ref and those derived from the reconstructed LAI rec time series using the several methods investigated.

Fig. 11 .
Fig. 11.RMSE as a function of Roughness observed over the reconstructed times series using the several methods.Color dots correspond to the values for the different methods over the 25 sites for cases with % gap in the selected classes of occurrence of missing data: 0-15 % (a) and 45-55 % (b).Note that EMD, LPF and AGF are only displayed for the lowest occurrence of missing data (0 % < % gap < 15 %).Black circles correspond to LAI comp and black diamonds to LAI ref .Lines correspond to the zero-offset linear regressions.

Fig. 13 .
Fig. 13.Accuracy of the start of season retrieval expressed in RMSE (days) as a function of the accuracy of LAI estimated expressed in RMSE.The same colors (corresponding to methods) and markers (corresponding to classes of % gap) as in Fig. 12 are used.

Table 1 .
List of the methods investigated.Length of processing window (whole means that the processing window is the whole time series) and maximum gap length tolerated are indicated.

Table 2 .
Sites selected for the study for the 5 biomes.The site number in BELMANIP2 ensemble of 420 sites, latitude and longitude and fraction of missing data (% gap) are indicated.

Table 3 .
Number of sites per vegetation class in BELMANIP2 set of sites, and number of cases considered in the gap simulation experiment.