Quality Assessment of MODIS Time Series Images and the Effect on Drought Monitoring

Drought Monitoring by remotely sensed moisture vegetation indexes is being an active research subject as the vegetation spectral responses are showed to be highly correlated to water content. The MODIS (MODerate resolution Imaging Spectro-radiometer) sensor of the Terra satellite provides MOD09A1 product of BRDF (Bidirectional Reflectance Distribution Function) used in computing moisture vegetation indexes (MVI). The exploration of an MVI time-series in the Kroumirie forest in Northern Tunisia showed important noise due to both clouds contamination and sensor defaults that had to be removed. Amongst methods for removing these imperfections, TIMESAT tool was designed for correcting time-series of satellite data and also to retrieve seasonal parameters from smoothed vegetation indexes. The methodology of smoothing functions to fit the time series data is based on two stages. First, a least square fit to the upper envelope of the vegetation indices series is applied. The second stage is achieved by local and adaptive fitting functions. The corrections have been made by spikes removal due to abrupt change of MVI variations and by fitting the MVI time-series to the upper envelop to correct the negative biases of remote sensing vegetation indexes. The adaptive SavitskyGolay function filter compared to local filtering process produces variations that conserve local variations for all the tested MVI. Seasonal vegetation parameters were extracted for each year of the time-series analysis and compared to the Standardized Precipitation Index (SPI) calculated at meteorological station level and for different time scales. Positive relations were found between SPI and the seasonal parameters expressed by the length and the amplitude of the season, indicating MODIS derived MVI sensitivity to water deficit or surplus conditions. The 6-month SPI showed the best performance when related to water sensitive indexes suggesting that MODIS derived indexes are more correlated to the precipitation variations over seasons.


Introduction
Global and regional monitoring of drought is becoming an active research subject to determine the impacts of drought on water management and socio-economic development, especially in limited water resources regions where drought episodes highly control water availability and functioning of forested and cultivated ecosystems.Various approaches for drought definitions and drought indexes were developed and synthesized in review studies such as [1].On the other hand, during the last two decades, the models used in the assessment of drought causes and manifestations integrated more and more indicators from multisensors satellite images to study their correlation with drought indexes.The monitoring of vegetation moisture by remote sensing became an active research subject in various studies within natural and cultivated areas where the vegetation spectral responses are showed to be highly correlated to physical indicators of water content and consequently to drought episodes.These studies that are based on linking spectral vegetation response to meteorological drought indexes, and are carried out to overcome sparse and inadequate network of weather stations.For example, [2] proposed a comprehensive analysis to evaluate the performance of remote sensing data towards drought conditions for fire prone vegetation in Australia.In [3] spectral vegetation response with thermal images and TRMM precipitation images (Tropical Rainfall Measuring Mission) was integrated in a unique drought spectral index that had been shown to be correlated with drought indexes such as SPI (Standardized Precipitation Index) and variation of crop yield in China.In the southern shore of the Mediterranean basin known as one of the most exposed regions to climate changes [4], there is a need of consistent drought monitoring to ensure a balanced ecohydrological functioning and to help mitigate climate changes effects.In [5] hydric stress and drought impacts on ecohydrologic balance were studied in a forested ecosystem of northern Tunisia.
Vegetation indexes (VI) derived from remote sensing, which are empirical combination of BRDF (Bidirectional Reflectance Distribution Function) in various wavelengths, contribute in detecting both biomass driven by chlorophyll activity and vegetation water content inducing various spectral responses in middle infra-red waves known as short wave infrared (SWIR).Time-series of remotely sensed VI are valuable sources of data for cultivated and natural vegetation landscapes environmental monitoring.In [6]  ).Many studies assessed the performance of MODIS (MOD erate resolution Imaging Spectro-radiometer) as a keystone satellite sensor system providing multiple products, such as MOD09A1 product, a near-daily coverage of the Earth surface, at 500 m in visible and infrared ranges.This product has also the advantage of two narrow discrete channels in the SWIR bands with a signal to noise ratio above 100 [7] which both could be useful for the monitoring of leaf water content [8].MODIS bands 5, 6 and 7 are well suited for canopy water monitoring because of the plant water sensitivity in these wavelengths combined with the high atmospheric water vapor transmittance [9].
Despite the wide use of MODIS sensors, many users reported that MODIS time-series can be subject to errors due to gaps, clouds and noise (e.g.[10] [11] [12]).Therefore an adequate use of this product requires a correction of remotely sensed time-series data.Amongst the various approaches of time-series images smoothing, local and global fitting methods are widely used.Local fitting methods are based on surrounding value of data in a time-series determined by median smoothing approaches [13] or by Savitzky-Golay filter approach [14] [15].Global approaches fit the data to long time scale of observations such as Fourier transform used in frequency analysis of a signal (for an extended review of filtering methods, refer to [11]).The elimination of noise by filtering process allow san estimation of vegetation seasonal patterns such as the length of the season and the seasonal amplitude.Therefore, most of time-series smoothing algorithms are compared for their phonologic metric detection such as start and end of the season [16] [17].They also could bring valuable information on drought occurrence as addressed in [18] where Fourier transform was used on a time-series of NDVI to identify specific frequency components corresponding to dry years.
In the present study, local approaches are tested within the tool TIMESAT developed by [15].This tool was used in various time-series filtering (e.g.[11] [12] [19]).Our main objective was to compare different smoothing algorithms of time-series moisture vegetation indexes computed from MODIS spectral bands in a forested Mediterranean region of northern Tunisia.These algorithms are embedded in the TIMESAT software that was used principally for anomalies corrections of the vegetation indexes time-series.Fitting functions used in the different algorithms allow the determination of the start and the end of vegetation season in the study region.We hypothesized that these comparative outputs of the season length and the amplitude reached by vegetation indexes could be adequate indicators of vegetation response to water deficits or surplus when compared to the SPI drought index.
The specific objectives of the present work are: 1) to make a time-series filtering of various moisture vegetation indexes (MVI) derived from the BRDF product MOD09A1 using the TIMESAT tool; 2) to evaluate the most performing filtering functions to make necessary corrections of the imperfections of the used data; 3) to determine the vegetation seasonality parameters for the years of the time-series of the study; 4) to investigate the relationship between seasonal parameters derived from filtered MVI time-series and the drought meteorological index SPI within the study ecosystem.

Study Area
The research focuses on the northern ecoregion of Tunisia known as the Kroumirie Forest, belonging to the Mediterranean North African forests (7˚52'E-37˚29'N to 9˚33'E-36˚22'N) (Figure 1).This ecoregion covering an area of 255,000 ha, is characterized by deciduous forest mostly represented by tree layer of cork oak (Quercus suber), shrub layer (e.g.Erica arborea and Pistacial entiscus) and herbaceous litter layer.This region has an average annual precipi- tation of 700 mm, but locally, precipitation can reach 1500 mm [20].Isohyets of annual average precipitation illustrated in Figure 1

Times-Series Precipitation Data and MODIS Images
The analysis period targeted in the present study covered the years between 2003 and 2009.Both 8-daystime-series of MODIS images and daily precipitations of four meteorological stations within the study region were acquired and processed.Meteorological data were available in four stations (illustrated in Figure 1) between 2002 and 2009, whereas MODIS images were available since 2003.

Drought Assessment by SPI
The Standardized Precipitation Index (SPI) [21], is a common measure of meteorological drought calculated from in-situ precipitation data registered at a given station.SPI calculation requires monthly precipitations at a local station to capture the cumulated deficit or surplus of precipitations over a given timescale.
SPI is based on fitting the total time-series precipitation to a gamma probability density function and, transforming the gamma distribution to a normal distribution with zero as mean and 1 as standard deviation.Values of SPI allow the characterization of the meteorological drought in a station as this index represents the number of standard deviations of an observed precipitation from the long term mean precipitation value calculated for a long term period.On the basis of computed SPI values, a drought category is associated to the station as shown in Table 1.SPI can be estimated at a station level for different time scales

MODIS Time-Series Data
The MOD09A1 product of the BRDF corresponds to 8 days with a spatial resolution of 500 m.Images corresponding to the period 2003-2010 were Table 1.Drought classification based on SPI ranges [21].Water Index [8], 2) NDII: Normalized Difference Infrared Index [23], and 3) GVMI: Global Moisture Vegetation Index [24].These MVI were developed from the combination of the near infrared (NIR) (B2 of MOD09A1) and the SWIR (B5, B6 and B7 of MOD09A1); the latter are particularly sensitive to the water content of vegetation [25] as shown in Figure 3.
In a previous research in the Kroumirie forest comparing field moisture data with MVI computed from MODIS [26], it has been shown that MVI are highly correlated to vegetation water content in this ecosystem.In other water limited ecosystems, a study led in Australia confirmed that MVI are sensitive to drought intensity and are highly correlated to SPI [27].The expression of the various MVI of the present study is given in Table 2.

Smoothing of MVI Time-Series by TIMESAT
TIMESAT software is dedicated to analyze time series satellite sensor data [15].
A free version of TIMESAT is available at

Methodological Stages of TIMESAT
The methodology of smoothing MVI time-series by TIMESAT is based on three stages to accomplish the correction of the synthesized as follows: The first step is based on a preliminary processing of data; it is applied to remove values of MVI corresponding to abrupt changes (spikes).This processing makes use of quality data available as additional band with the MODIS product.
Secondly, a least square fits to the upper envelope of the MVI is used to overcome the fact that vegetation indexes generated from remotely sensed data are negatively biased as proved for example for NDVI (Normalized Difference Vegetation Index) time series by [29].
For a time-series of vegetation index, data are organized in two-dimensional arrays where each array corresponds to a vegetation index at a given time.The fitting function is a weighted combination of the original vegetation index where weights could be adjusted during the processing.In the beginning of this stage, weights are derived from data quality product.Next, the weights are adjusted so that the highest data are influencing the correction process (values with small weights influence the fit less than values with large weights).Figure 5 shows the adaption the upper envelope through weights varying from a multi-step

Filtering Methods
The filters in TIMESAT tool are of two types: 1) local model fitting based on minima and maxima, and 2) adaptive filtering.
In ( ) ( ) where (c 1 , c 2 ) determine the base level and the amplitude of the fitting function f(t), "x" is the non linear parameter determining the shape of the basis function and "t" is the time variable.
( ) where "x 1 " determines the position of the maximum or minimum with respect to the independent time variable t, (x 2 , x 3 ) determine the width and flatness of the right function half, (x 4 , x 5 ) determine the width and flatness of the left function half.
( ) where "x 1 " determines the position of the left inflection point, "x 2 " gives the rate of change, (x 3 , x 4 ) give the same parameters for the right inflection point.
The second method is known as the adaptive Savitzky-Golay SG filter ( [31] [32]) and uses a linear combination of nearby values in a window to fit locally a polynomial function by a least square method.The filtered value is then calculated to fit to a polynomial function while varying locally the size of the window to avoid getting very important variations between original value and filtered value.For a VI time-series, if we assume that at a date "t", VI(t) is noted y i , this y i is replaced by a linear combination of nearby values of a window where "n" represents the size of the window.Equation (4) gives the SG filtered value of y i : where y i is the original data series representing one date value of VI, C j is the weight of the j th original data value y j in the smoothing window size equal to (2n + 1).For example, for n = 1, the window size convolution is equal to 3, meaning that VI (t) is replaced by the filtered value combining VI(t − 1), VI(t) and VI(t + 1).
The window size "n" can be adjusted to avoid obtaining large variations around a given VI value.

Seasonal Vegetation
Considering a time-series of VI with one growing season per year, seasonality parameters resulting from fitting functions can be extracted for each full season of the time-series.Various parameters could be computed from the filtered MVI as shown in Figure 6: the start, the end and the length of the season; the seasonal amplitude and the integral of the function describing the season.These seasonal parameters allow the comparison of the vegetation response to water availability offering the possibility to investigate and explain inter-annual variations of the ecosystem vegetation status over space and time.

Setting of TIMESAT Parameters
The parameters that need to be set in the beginning of time-series filtering by TIMESAT are illustrated by the interface example of Figure 7.The main parameters that have to be specified are: • Amplitude cutoff value: time-series with smaller amplitude than the cutoff will not be processed; this parameter could be set to 0 to process all data.• Spike method: the median filter method assigns a zero weight to values that are substantially different from both the left and right-hand neighbors and from the median in a window.The difference from the median is measured in units of the standard deviation of the time-series that we set to 2 ("Spike Parameters" =2).
• Number of envelope iterations: this setting allows two additional fits where the weights of the values below the fitted curve is decreased forcing the fitted function toward the upper envelope.• Adaptation strength: this number indicates also the strength of the upper envelope adaptation.After some trials, we set the value 3 for envelope iterations and the strength adaptation was set to 2.

Analysis of MVI Time-Series Filtering
The moisture vegetation indexes MVI were computed from the BRDF product covering the time-series period from January 2003 to December 2010.In each year, 46 images were acquired and processed on the basis of one image each 8 days.Thus the total images proceeded in this data set are 368.Next, each timeseries vegetation index had been processed in TIMESAT where local and adaptive filtering functions were tested.Results were first analyzed based on visual exploration of time-series profiles at various locations in the Kroumirie Forest region.Figure 8 shows original data and the three fitting functions of the NDII6 time-series for "Ain Draham" and "Ouchtata" sites representing maximum and local variations for all the tested MVI.This has also been reported in other works where the SG adaptive filter has been successfully used in minimizing noise in NDVI time-series [33] and [12] and Leaf area index time-series [11].This results in more obvious limits of start and end of the season.
The year 2003 presented the highest annual precipitations amongst the precipitation time-series as illustrated in Figure 2; this explains that computed MVI for this year show highest values compared to the other years.

Drought Analysis by SPI and MODIS Time Series
The SPI for four meteorological stations within the study region were computed for each month of the period where daily precipitations are available From smoothed MVI time-series by TIMESAT, the seasonal vegetation parameters were extracted around the positions of the meteorological stations and compared to the SPI variations.The length of the season and the amplitude of the smoothed MVI vary from one year to another suggesting that water deficit/ surplus conditions influence the seasonal vegetation parameters.In particular, the integral of the MVI-function describing the season from the start to the end L-INTEG-MVI (denoted "h" in Figure 6) was extracted as it integrates both the amplitude and the length of the vegetative season, and could be an interesting indicator of precipitation deficit or surplus effects on the canopy.
For the explored stations, we noticed the existence of positive relations between SPI and L-INTEG-MVI indicating that it is possible to make a monitoring of seasonal variation response to climate by remote sensing.As a first approximation, a linear relation between SPI and L-INTEG-MVI was evaluated and the R 2 coefficients were computed for each tested MVI (Table 3).Most significant results were obtained with 6-month SPI reflecting medium-term trends in precipitation and more rarely with 3-month SPI which reflects short-term moisture conditions.This suggests that water sensitive indices are more correlated to the precipitation variations over seasons.
Results show differential capacity of MVI to detect water surplus or deficit.
Indeed in [26], it has been shown that in this peculiar ecosystem characterized by a gradient of vegetation density, the models linking water-sensitive indexes and the hydric status of the canopy are highly dependent of the biomass expressed by the LAI (Leaf Area Index).This could explain the variable performance of MVI from one site to another as they had different sensitivities to the biomass as reported in other works [27].In a second phase we tested a second degree relation between SPI and L-INTEG-MVI.Results reported in (Table 4) show that we have substantially the same relations but with increasing R 2 coefficients.Figure 10 2, Figure 3).On the other hand, the NDII and the GVMI indexes, computed from MODIS B2, B6 and B7, generate low but not zero values of L-INTEG-MVI for these years.The amplitude variation of L-INTEG-NDWI is small compared to that computed from NDII6, NDII7, GVMI6 and GVMI7 (Figure 10).Therefore, despite the high R 2 coefficients values obtained with first and second degree relations between SPI and L-INTEG-NDWI, this B5-based MVI is not recommended for vegetation response to climate conditions in the study region.
From the previous analysis, our hypothesis that the season length and the amplitude reached by remotely sensed vegetation indexes could be adequate indicators of vegetation response to water deficits or surplus is validated in the sites where SPI were calculated.It will be possible to improve these preliminary results by spatial correlations between SPI extended to the whole region and the seasonal parameters extracted from the fitted MVI time-series fitted functions.

Conclusions
This study presented moisture vegetation indexes time-series corresponding to The contribution of this study is that the hypothesis of correlation between time-series vegetation indexes, seasonal parameters and drought conditions had been addressed and verified in the Mediterranean peculiar ecosystem.The interest of remote sensing indexes lies in the fact that the information is already spatialized and allows a monitoring of an ecosystem response to drought conditions via time-series moisture vegetation indices after their corrections by adequate filtering methods.These preliminary results should be improved by in-situ data on seasonal parameters as well as spatial correlations between SPI and spectral vegetation indexes extended to the whole region.
an exhaustive literature review assessment was presented for prominent applications based on time-series of vegetation indexes.The time-series NDVI products have been widely used as different sensors afforded these data such as the 8-km resolution AVHRR-GIMMS (Advanced Very High Resolution Radiometer-Global Inventory Modeling and Mapping Studies) launched in the end of the eighties and the 1-km spatial resolution SPOT-VEGETATION program launched in the nineties.Since the early 2000 s, the availability of MODIS products has brought a breakthrough in time-series data with the 250-m spatial resolution NDVI product (MOD-13Q1 show an important gradient of precipitation decrease in the direction NW-SE.The Mediterranean climate of the study region is characterized by rainy autumns and winters where precipitations are distributed between September and March.The deciduous forest of Kroumirie region is highly detected by MODIS sensors.The processing of the Enhanced Vegetation Index (EVI) (MOD13Q1 bi-weekly product: 250 m spatial resolution) allows the discrimination of agriculture (EVI < 0.1) and forested areas that we classify in 3 density classes as illustrated in Figure 1.

( 3
months, 6 months, 12 months and more) and for each month of the year.3-monthSPI provides short and medium-term moisture conditions whereas 6-month SPI shows effects over seasons.12-month SPI and more are used to analyze long term precipitations patterns.In the Kroumirie region, records of four meteorological stations were used to calculate the SPI from monthly precipitations of the period 2002-2009, except for Feija station where records of 2002-2003 were missing (Figure2).

Figure 2 .
Figure 2. Annual precipitations in meteorological stations of the Kroumirie forest.(Data source: DGRE: General Department of Water Resources).
[28].Thistool was designed primarily for analyzing time-series of satellite data and to extract seasonal information from any kind of time series.The program provides different smoothing functions to fit the time-series satellite images.Figure4represents the general framework of MODIS time-series integration and processing stages in TIMESAT.
the first case, data are fitted to local model functions defined in intervals around maxima and minima in the time-series data.The local model functions have the general form given by Equation (1).The local fitting function uses two types of functions which are: the asymmetric Gaussian function Equation (2) and the double logistic function Equation (3).

Figure 6 .
Figure 6.Seasonality parameters extraction.Points (a) and (b) mark, respectively, start and end of the season; Points (c)and (d) give the 80% levels; (e) displays the point with the largest value; (f ) displays the seasonal amplitude and (g) the seasonal length.Finally; (h) and (i) are integrals showing the cumulative effect of vegetation growth during the season [30].

Figure 7 .
Figure 7.An example of parameters setting for filtering an NDII6 eight years time-series in TIMESAT tool.
minimum time-series annual precipitation.The other MVI (non illustrated here) represent the same patterns of fitting functions with variations in the maximum reached by the MVI and the amplitude of the function over time.Compared to original values of NDII6, the fitting function process eliminated the extreme values (spikes) resulting in sensor imperfection or clouds.Moreover, data were fitted to the upper envelope to compensate for the negative biases of vegetation indexes.The SG adaptive filter Figure 8(a), Figure 8(c) compared to local filtering process Figure 8(b), Figure 8(d) produce variations in time-series that conserve

Figure 8 .
Figure 8. Time-series plot of NDII6 values.Adpative SG filter for (a) Ouchtata (c)Ain Draham.Local filters for (b) Ouchtata (d)Ain Draham.(FMT1 is SG sadaptive filter; FMT2 is the asymmetric Gaussian function; FMT3 is the double logistic function. (2002-2009), except for Feija station where records of 2002-2003 were missing.The SPI computed at the station level indicate that the study region experienced episodes of dry and wet conditions of different magnitude.To assess the possible relation between MVI derived from MODIS and drought intensity, the SPI values for the month of March were considered.For example, computation of 3-month SPI of March in a given year is based on the cumulative precipitations of January, February and March and it is compared to the same 3-month period of the time series years.Figure 9 represents the SPI variations for 3, 6 and 12 months time scales showing the levels of drought affected to each station during the analysis period.From this analysis, we can conclude that the seasons 2002-2003 and 2008-2009 are the wettest ones, whereas for the other seasons, conditions were rather normal or moderately dry.
illustrates for each station the SPI and the time-course of the L-INTEG-MVI having the highest R 2 coefficients.The strong relation between L-INTEG-NDWI, where NDWI is computed from MODIS B2 and B5, and the 6-month SPI is explained by nil values of the

Figure 10 .
Figure 10.SPI at multiple time scales and L-INTEG-MVI variation between 2003 and 2009.

8 -
days images acquired from 2003 to 2009 and calculated from MODIS spectral bands after anomalies corrections by filtering functions in TIMESAT tool.All the tested filtering functions had eliminated spikes and the time-series data were fitted to the upper envelope thus overcoming the negative biases of vegetation indices.Amongst the three methods tested in this tool, it has been shown that adaptive Savitzky-Golay fitting function is more appropriate because it simultaneously corrects the imperfections represented by spikes while keeping the variations due to a natural evolution of the moisture vegetation indices.Local filtering asymmetric Gaussian function and double logistic function capture less local variations; however they generate more obvious limits of start and end of vegetative season.Analysis of the SPI computed from an eight-year history (2002-2009) of precipitation records demonstrated that smoothed moisture vegetation indexes have the capacity to detect wet and drought conditions characterized either by a deficit or a surplus of precipitations.Seasonal functions integrating both the amplitude and the length of the vegetative season were analyzed in comparison to the SPI within each station at different time scales.Regression coefficients of the (SPI, L-INTEG-MVI) linear and polynomial positive relations show different performances for the tested MVI, but at least one MVI performs for each tested station.Moisture vegetation indexes computed from MODIS short wave infra red bands 6 and 7 outperformed the index computed from band 5.The best agreement between drought index SPI computed for three time scales (3, 6 and 12 months) and the seasonal parameters were with 6-month SPI suggesting that MODIS derived indices are more correlated to the precipitation variations over seasons.

Table 3 .
Regression coefficients R 2 of the linear relation between 6-month SPI (except underlined values obtained for 3-month SPI) and L-INTEG-MVI.

Table 4 .
Regression coefficients R 2 of the second degree relation between 6-month SPI (except underlined values obtained for 3-month SPI) and L-INTEG-MVI.