forest: 0.1 < EVI < 0.2; Medium density forest: 0.2 < EVI < 0.3; High density forest: 0.3 < EVI < 0.38). Isohyets are computed from 30 years observations 1960-1990 from INM (National Institute of Meteorology).
3. 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 proce- ssed. Meteorological data were available in four stations (illustrated in Figure 1) between 2002 and 2009, whereas MODIS images were available since 2003.
3.1. Drought Assessment by SPI
The Standardized Precipitation Index (SPI)  , 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 (3 months, 6 months, 12 months and more) and for each month of the year. 3-month SPI 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 (Figure 2).
3.2. 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  .
Figure 2. Annual precipitations in meteorological stations of the Kroumirie forest. (Data source: DGRE: General Department of Water Resources).
provided as Hierarchical Data Format-Earth Observing System (HDF-EOS) downloaded from the online USGS Global Visualization Viewer (GloVis)  . Data were selected and ordered within a graphical interface. Tiles covering North Africa region were clipped to the extent of Tunisia. The images were rectified from sinusoidal projection to geographic projection by the use of MODIS Reprojection Tool, and next they were converted to Universal Transverse Mercator (UTM, zone 32 North). From visualization and profiles examination, we noticed important noise due to clouds contamination or sensor defaults that had to be removed. We also noticed that band 5 of the MOD09A1 product presents a regular strip effect.
Monitoring of vegetation moisture is possible through water sensitive indexes (MVI) based on BRDF combination in the various spectral wavelengths. The most used MVI reported in the literature are 1) NDWI: Normalized Difference Water Index  , 2) NDII: Normalized Difference Infrared Index  , and 3) GVMI: Global Moisture Vegetation Index  . 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  as shown in Figure 3.
In a previous research in the Kroumirie forest comparing field moisture data with MVI computed from MODIS  , 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  . The expression of the various MVI of the present study is given in Table 2.
4. Smoothing of MVI Time-Series by TIMESAT
TIMESAT software is dedicated to analyze time series satellite sensor data  . A free version of TIMESAT is available at  . 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. Figure 4 represents the general framework of MODIS time-series integration and processing stages in TIMESAT.
Figure 3. Leaf spectral reflectance simulated for MODIS and LANDSAT bands. Response in SWIR bandsis sensitive to relative water content (RWC)  .
Table 2. Moisture vegetation indexes based on MODIS spectral bands (MOD09A1).
4.1. 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  .
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
Figure 4. Flowchart of MODIS time-series processing in TIMESAT.
Figure 5. Fitted functions to the upper envelope of NDVI time series. Thin solid line reprsents the original vegetation index data and the thick line displays the fitted function. (a) Fitted function from the first weights; (b) Fitted function from decreased weights for low NDVI  .
procedure. The first step results in a fitted function affected by original lowest values of NDVI Figure 5(a), whereas Figure 5(b) shows the effect of low data weight decrease on the final fitted function.
The third stage of time series corrections is achieved by one of the three filtering methods: 1) adaptive Savitsky-Golay function, 2) local model Gaussian asymmetric function, and 3) local model double logistic function.
4.2. 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 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).
where (c1, c2) 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 “x1” determines the position of the maximum or minimum with respect to the independent time variable t, (x2, x3) determine the width and flatness of the right function half, (x4, x5) determine the width and flatness of the left function half.
where “x1” determines the position of the left inflection point, “x2” gives the rate of change, (x3, x4) give the same parameters for the right inflection point.
The second method is known as the adaptive Savitzky-Golay SG filter (   ) 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 yi, this yi 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 yi:
where yi is the original data series representing one date value of VI, Cj is the weight of the jth original data value yj 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.
4.3. 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.
4.4. 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.
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  .
Figure 7. An example of parameters setting for filtering an NDII6 eight years time-series in TIMESAT tool.
・ 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.
5. Results and Discussion
5.1. 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 time- series 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 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 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  and  and Leaf area index time-series  .
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.
The local and adaptive filtering functions show seasonal cycle of typical deciduous forests with one curve by year. Local fitting functions (the asymmetric Gaussian function and the double logistic function) are fit to data in intervals around maxima and minima, and therefore they capture less local variations. 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.
5.2. 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 (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.
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 R2 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  , 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  .
Figure 9. SPI values for the period 2002-2009 calculated for 3, 6 and 12 months.
Table 3. Regression coefficients R2 of the linear relation between 6-month SPI (except underlined values obtained for 3-month SPI) and L-INTEG-MVI.
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 R2 coefficients. Figure 10 illustrates for each station the SPI and the time-course of the L-INTEG-MVI having the highest R2 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
Table 4. Regression coefficients R2 of the second degree relation between 6-month SPI (except underlined values obtained for 3-month SPI) and L-INTEG-MVI.
Figure 10. SPI at multiple time scales and L-INTEG-MVI variation between 2003 and 2009.
integral function of this index for moderately dry years. The fact that NDWI vanishes when climate conditions are normal to dry is explained by the vegetation spectral response in bands B2 and B5 used for its computation which are relatively close in these bands (Table 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 R2 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.
This study presented moisture vegetation indexes time-series corresponding to 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.
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.
The author wishes to thank the Ministry of Agriculture and Water Resources of Tunisia for providing climatic data (especially Mrs H. Ben Mansour from the DGRE). The author also thanks NASA Earth Science Enterprise for providing free MODIS images on the web. The author presents a special thank to the developers of the TIMESAT tool from Lund University, Sweden.