Analysis of the Effects of Drought on Vegetation Cover in a Mediterranean Region through the Use of SPOT-VGT and TERRA-MODIS Long Time Series

Abstract

The analysis of vegetation dynamics and agricultural production is essential in semi-arid regions, in particular as a consequence of the frequent occurrence of periods of drought. In this paper, a multi-temporal series of the Normalized Difference of Vegetation Index (NDVI), derived from SPOT-VEGETATION (between September 1998 and August 2013) and TERRA-MODIS satellite data (between September 2000 and August 2013), was used to analyze the vegetation dynamics over the central region of Tunisia in North Africa, which is characterized by a semi-arid climate. Products derived from these two satellite sensors are generally found to be coherent. Our analysis of land use and NDVI anomalies, based on the Vegetation Anomaly Index (VAI), reveals a strong level of agreement between estimations made with the two satellites, but also some discrepancies related to the spatial resolution of these two products. The vegetation’s behavior is also analyzed during years affected by drought through the use of the Windowed Fourier Transform (WFT). Discussions of the dynamics of annual agricultural areas show that there is a combined effect between climate and farmers’ behavior, leading to an increase in the prevalence of bare soils during dry years.

1. Introduction

In semi-arid regions, drought is a frequent phenomenon leading to serious problems in agriculture and food safety, and increasing attention has been drawn to its high attendant economic and social costs [1]. Although it is essential to quantify the occurrence of droughts, this is generally found to be difficult since they are spatially variable and context-dependent [2]. In general, in situ meteorological data from weather stations is analyzed through the use of drought indices [3,4] such as the Palmer Drought Severity Index (PDSI; [5]) or the Standardized Precipitation Index (SPI; [6,7]). The PDSI is based on long-term historical precipitation and mean temperature data. The SPI is based on precipitation data only. In areas with a high density of weather stations, drought conditions can be interpolated with only small associated errors. In the case of regions with sparsely distributed weather stations, statistical techniques such as the Inverse Distance Weighted method or a stochastic model of Ordinary Kriging can be considered. These methods are not completely sufficient to ensure that drought quantification is achieved with high accuracy, in areas where limited ground measurements are made. The problem is more complicated in areas of the globe in which weather stations are completely absent. For this reason, in the last two decades scientists have proposed various new approaches for the estimation of drought conditions, based on the use of remotely sensed satellite data [8,9,10,11].

A number of different indices based on data recorded by remote optical sensors have been developed to quantify drought [12,13,14]. In fact, remote sensing allows soil and vegetation dynamics and its variations to be monitored over time [15,16,17]. Most studies are based on the use of the Normalized Difference Vegetation Index (NDVI), which has demonstrated good accuracies for the quantification of green vegetation cover or vegetation abundance [18]. This index is expressed by: NDVI = (RNIR − RRED)/(RNIR + RRED), where RNIR is the near-infrared (NIR) reflectance and RRED is the red reflectance. A high NDVI indicates a strong level of photosynthetic activity [19,20].

The Vegetation Condition Index (VCI) has been tested at different sites, revealing its high potential for the detection and monitoring of drought. Kogan et al. [12], Seiler et al. [21], and Quiring et al. [22] observed a high correlation between the VCI and agricultural production in different regions of the globe (Africa, America, Europe). Gitelsen et al. [23] showed that when tested at six different sites in Kazakhstan, the VCI could explain 76% of the observed variations in crop density. Amri et al. [14] proposed the use of the notion of the Vegetation Anomaly Index (VAI), derived from SPOT-VGT data, and observed a high correlation between this product and precipitation levels measured at the site studied in the present paper. Mu et al. [24] proposed the Drought Severity Index (DSI), combining the Vegetation Anomaly Index and an evapotranspiration anomaly index, which has shown considerable potential for the quantification of drought in different climatic regions.

In North Africa, droughts are the cause of major losses in non-irrigated agriculture, and crop yields are severely reduced during dry years [14]. Droughts are also an important factor in environmental degradation, because they limit the development of vegetation cover and make the soil more sensitive to erosion in the case of intense precipitation.

The aim of this paper is to propose an analysis of variations in vegetation cover, in particular for the case of annual agricultural land use in a semi-arid region of North Africa, in the context of drought events. It is based on two satellite products having different spatial resolutions, namely SPOT-VGT and TERRA-MODIS. Section 2 presents the studied site and the database used in this study. Section 3 compares the behavior of inter-annual NDVI time series, and provides a discussion of land use mapping and NDVI anomalies from TERRA-MODIS and SPOT-VGT observations. Section 4 describes the specificities of agricultural areas confronted by drought events. Our conclusions are presented in Section 5.

2. Database

2.1. Studied Site

The studied site is the Kairouan plain, which is situated in central Tunisia (35°–35°45′N; 9°30′–10°15′E) (Figure 1) and is characterized by a semi-arid climate [25]. The average annual rainfall is approximately 300 mm per year, with the rainy seasons lasting from October to May, followed by dry summers. The rainfall patterns in this semi-arid area are highly variable over time and space. As shown in the annual precipitation plot of Figure 2, this site is frequently affected by drought events: during the 15-year period of the present study, from 1998 to 2013, we identify five dry years with an accumulated precipitation of less than 200 mm. The mean daily temperature in Kairouan City ranges between a minimum of 10.7 °C in January and a maximum of 28.6 °C in August, with a mean value equal to 19.2 °C. The mean annual potential evapotranspiration (Penman) is close to 1600 mm. The landscape has no relief, and land use is dominated by agriculture, with two main types of vegetation cover: annual agriculture and olive trees.

2.2. Satellite Data

(a) SPOT-VGT

The 10-day synthesis (S10) products derived from SPOT-VGT data are available at full resolution (1 km), and include the 10-day NDVI data used in the present study. For these products, top-of-atmosphere corrections were applied using the SMAC algorithm, which corrects for molecular and aerosol scattering, water vapor, ozone, and other gas absorption effects. The parameters taken into account in the atmospheric corrections are the aerosol optical depth (AOD), atmospheric water vapor, and ozone, and a Digital Elevation Model used for atmospheric pressure estimation. The water vapor parameter is derived from data provided by Météo-France once every six hours, with a 1.5° × 1.5° grid cell resolution. Although the AOD is currently retrieved from B0 data (blue spectral band of SPOT-VGT, 0.43–0.47 µm), combined with the NDVI, prior to 2001 it was a static dataset, which varied as a function of latitude only. Different systematic errors (misregistration of the different channels, calibration of the linear array detectors for each spectral band) are corrected in the final P product, which is re-sampled to a plate carree geographic projection. The S10 products are available at http://free.vgt.vito.be. In the present study, we make use of this product for the period between September 1998 and August 2013.

(b) TERRA-MODIS

In this study we used 250-m spatial resolution, 16-day composites of MODIS NDVI data (MOD13Q1, collection 5). This product is retrieved from atmosphere-corrected, daily bidirectional surface reflectance observations, using a compositing technique based on product quality assurance metrics to remove low quality pixels. The 250-m spatial resolution is the finest available from the MODIS NDVI dataset, and the 16-day composite was selected in order to ensure a high probability of having the best quality pixel (reduced cloud effects) representing the NDVI within each 16-day period. Here, we make use of this product in order to analyze the influence of drought, during the period between September 2000 and August 2013.

2.3. Precipitation

Precipitation estimations were based on measurements recorded by a network of 30 rain gauges, distributed over the entire site. As our studied site is mainly flat, without mountains, the terrain was assumed to have no influence on the rainfall’s spatial distribution. The Inverse Distance Weighting (IDW) interpolation algorithm, which computes values at non-sampled points by estimating the weighted average of data observed at nearby points, was used to estimate daily precipitation maps [26]. The precipitation time series and NDVI data set are thus available with similar spatial resolutions. As a consequence of the small number of rain gauges and the assumption of rainfall with no particular spatial distribution, the difference between precipitation products at 1 km and 250 m resolutions is very small.

2.4. NDVI Temporal Series

Figure 3 shows the NDVI time series over the last 13 years (2000–2013), for three types of vegetation: pastures, annual agricultural areas, and olive groves, based on data provided by the SPOT-VGT and TERRA-MODIS sensors. For each sensor, the NDVI and precipitation data correspond to an average value, taken over a selected area characterized by homogenous vegetation cover. It is retrieved using high resolution land use maps and direct observations in the site. Among the three types of vegetation cover, the NDVI with the lowest dynamic range is observed for olive groves. These are characterized by a dispersed vegetation cover, with a spacing of approximately 15 m between two successive trees accompanied by a high percentage of bare soil. The highest temporal dynamics are observed in the areas characterized by annual agricultural vegetation, where the NDVI can be seen to vary strongly between wet and dry years. For example, a maximum NDVI equal to 0.6 is observed during 2010–2011, the wettest year (430 mm precipitation), whereas 0.1 is observed during the period 2000–2001 (176 mm precipitation), the driest hydrological year. In all land use cases, a high level of agreement is found between the NDVI dynamics observed with SPOT-VGT and TERRA-MODIS. However, during two different one-year periods (2004–2005 and 2006–2007) some discrepancies were observed in the case of the annual agriculture vegetation cover, for which the maximum value of NDVI derived from TERRA-MODIS data is clearly lower. This behavior could be attributed to differences in the data processing algorithms, in particular for cloudy conditions. It can also be noted that the precipitation recorded between November and January, which is essential for growth of annual agriculture, was relatively low (approximately 105 mm) during these two periods, when compared to that recorded during wet years (170 mm for 2003–2004 and 190 mm for 2005–2006). It is thus likely that the NDVI values retrieved from the TERRA-MODIS data are more suitable for the accurate determination of the vegetation dynamic.

3. Analysis of Vegetation Dynamics with SPOT-VGT and TERRA-MODIS NDVI Products

3.1. Temporal Dynamics of Land Use Mapping

In order to accurately analyze the vegetation dynamics corresponding to the different types of vegetation cover present at the studied site, it is essential to consider annual land use maps. As a consequence of the lack of high-resolution satellite images, suitable for the retrieval of land-use maps for all of the studied years between 1998 and 2013, we adopt an approach based on the use of low spatial resolution SPOT-VGT and TERRA-MODIS NDVI images. For this, we consider three characteristic classes, labelled: “olive trees and bare soil”, “annual agriculture” and “pastures”. As shown by the NDVI time series, olive groves are affected by quite small variations in NDVI, between dry and wet years, with the NDVI remaining very close to that of bare soil (between 0.15 and 0.2). The similarity of these values is due to the large average spacing between olive trees (between 15 and 20 m), as well as to the continuous use of soil tilling (up to five times per year), a practice that is commonly used to increase infiltration and eliminate the growth of grass and weeds following rainfall events. For this reason, we group olive groves and bare soils together in the same class of land use.

A linear mixing theory developed by Settle et al. [27] is considered, in which it is assumed that the reflectance (respectively NDVI) of a mixed pixel is given by the sum of the mean reflectance (respectively NDVI) values of the different land cover classes within the pixel, weighted by their respective fractional covers. The first step in this approach involves the identification of typical NDVI profiles, representative of each of the above land cover classes in the disaggregation methodology. For this step, we consider information related to the class composition of each pixel, retrieved from a high-resolution land-cover map. Disaggregation techniques are designed to estimate the annual proportion (between 0 and 1) of specific classes occurring within each pixel, by considering the NDVI time series acquired for each year. The result is a set of three fractional images, each corresponding to one specific class of land-cover. While this information describes the composition of each pixel, it does not provide any indication as to how the classes are spatially distributed within each pixel. Figure 4 illustrates the results of the three land use mapping using SPOT-VGT and MODIS data for the wettest (2010–2011) and driest (2000–2001) years. A qualitative comparison of the SPOT-VGT and MODIS land use maps shows that they are in good agreement. Approximately the same areas are retrieved, corresponding to a strong presence of “annual agriculture” or “pasture” or “olive groves and bare soil”. During dry years, a clear decrease is observed in the presence of annual agriculture. In order to verify the accuracy of the land use maps produced from low-resolution satellite images, we compared high-resolution SPOT-HRV maps of annual agricultural areas with the proposed low-resolution SPOT-VGT and TERRA-MODIS NDVI maps over a period of eight years, for which high-resolution satellite data exists.

Table 1. High resolution optical data (SPOT, LANDSAT) acquisitions.

Agricultural Year	Date of Acquisition	Satellites
1999/2000	5 January 2000	SPOT 2
	13 March 2000
	9 July 200
2005/2006	16 November 2005	SPOT 2
	26 March 2006
	22 July 2006
2006/2007	21 December 2006	SPOT 2/LANDSAT
	25 March 2007
	27 July 2007
2008/2009	21 December 2008	SPOT 5
	14 April 2009
	20 May 2009
	17 July 2009
2009/2010	8 October 2009	SPOT 5
	30 December 2009
	28 March 2010
	24 June 2010
2010/2011	24 December 2010	SPOT 5
	29 January 2011
	17 March 2011
	03 July 2011
2011/2012	6 November 2011	SPOT 5
	13 January 2012
	31 March 2012
	6 July 2012
2012/2013	26 December 2012	SPOT 5
	21 January 2013
	19 March 2013
	11 July 2013

shows the SPOT-HRV and LANDSAT high-resolution images used to produce the high-resolution land-use maps, for the eight agricultural seasons.

Figure 5 shows the high-resolution land use map of the studied site for the 2012–2013 agricultural season. An approach making use of decision tree learning was applied, based on the use of four multi-temporal data acquisitions for each year [25]. This allowed 11 land use classes to be identified: rivers, dams, sebkhas, reliefs, urban areas, dry olive, summer and winter vegetables, cereals, bare soils, and tree crops. Our annual agriculture class derived from SPOT-VGT and MODIS maps includes several classes derived from the SPOT-HRV land use maps, such as cereals and summer and winter vegetables, which were merged in order to compare high- and low-resolution estimations. The comparison was made according to the SPOT-HRV scene data only, and other regions of the Kairouan plain were masked in the SPOT-VGT and MODIS land use maps to allow comparisons to be made correctly. As shown in Figure 6, in the case of the total estimated area of annual agriculture land use, a better agreement is found between SPOT-HRV and TERRA-MODIS data than between the SPOT-HRV and SPOT-VGT data, with correlation coefficients equal to 0.89 and 0.59, respectively. This result could be explained by differences in spatial and temporal resolutions. As a result of the relatively small size of the agricultural fields (generally between 1 and 5 ha), the TERRA-MODIS outputs are found to be in better agreement with the high-resolution reality of land use at the study site.

3.2. Discussion of NDVI Anomalies with TERRA-MODIS and SPOT-VGT Data

In this section, we propose to use an anomaly index applied in [14], which provides a quantitative illustration of vegetation stress and the influence of drought on the vegetation cover. This index is based on statistics derived from the NDVI time series, and is referred to as the “Vegetation Anomaly Index” (VAI), written as:

$$VA{I}_i=\frac{NDV{I}_i-{\left(NDV{I}_i\right)}_{mean}}{{\sigma }_i}$$

(1)

where $NDV{I}_i$ is the NDVI estimate for a given month i, ${\left(NDV{I}_i\right)}_{mean}$ is the mean value of the NDVI during month i, derived from the previously described 15 years of NDVI time series, and ${\sigma }_i$ corresponds to the standard deviation of the NDVI values estimated for month i over the same 15 year period.

Figure 7 illustrates the mean level of this index for the studied site, using data from both SPOT-VGT and TERRA-MODIS for the period of interest, between 2000 and 2013. The VAI is computed for the hydrological years beginning in September, over the studied region. We first retrieve approximately the same VAI behavior for these two products and positive or negative levels for approximately the same periods. However, various discrepancies are encountered due to the different spatial resolutions of these two sensors. The VAI index can be seen to have a positive trend for both satellite products, which can be firstly explained by a slightly positive trend in autumn rainfall (between November and January) observed during the 13 years of the studied period. In addition, a small increase in the surface of irrigated areas was observed during the last 10 years of the studied period, which may also have contributed to this positive trend [28]. This outcome confirms the trend reported in a global study of drought by [24,29], based on the interpretation of the NDVI and DSI (Drought Severity Index) indices, showing that the studied region (and North Africa in general) is becoming wetter. We find a more significant trend for the SPOT-VGT time series, with a slope equal to 0.13, whereas the TERRA-MODIS data leads to a slope of 0.06. This difference is due mainly to the discrepancy between the annual maximum values of NDVI determined for these two products, during the one-year periods between 2004–2005 and 2006–2007 (see Figure 3).

Figure 8 compares the mean annual VAI indices for the two products, for the 13 years considered in this study. The two products are found to be in very good agreement, and an approximately linear relationship can be established between the two indexes:

$$VA{I}_{SPOT-VGT}=\alpha VA{I}_{MODIS}+\beta$$

(2)

Table 2. Parameters derived from linear regression between VAI_SPOT-VGT and VAI_MODIS over the studied site, for the period between 2000 and 2013.

Month	α	β
September	0.78	0.0075
October	1.05	0.01
November	0.88	0.03
December	0.96	0.06
January	1.02	0.02
February	0.99	−0.02
March	1.08	−0.05
April	1.07	0.0028
May	0.87	0.1
June	1.04	0.066
July	1.07	0.07
August	0.96	0.06

lists the linear regression parameters retrieved for each month of the year, showing that for most months α is close to 1, and α is close to zero.

4. Dynamics of the Area Used for Annual Agriculture

In the previous sections we propose a method for the quantification of drought, through the use of satellite time series. The area used for annual agriculture is certainly the most sensitive to drought events, as shown in Figure 3. It is essential to gain a clear understanding of the influence of drought in this area, in order to improve its agricultural productivity. In this section, following the analysis of the temporal variations observed in the NDVI data, we discuss the behavior of the NDVI signals in the spectral domain through the use of the windowed Fourier transform, thus allowing the specificities of extremely dry years to be analyzed.

4.1. Analysis of the Vegetation Using a Windowed Fourier Transform

The Fourier Transform (FT) can be used to analyze the frequency content of a signal in the time domain, by decomposing the original signal into a set of superimposed sine and cosine basis functions. By definition, these basis functions are not local—i.e., they have an infinite span and are globally uniform over time. Thus, the FT does not represent any temporal variability that may be present in the signal, and cannot provide any information about the time of occurrence of any specific frequential event. For such applications, the Windowed Fourier Transform (WFT) allows frequency components present in a time series to be localized in the time domain. This technique applies a fixed-width moving window over the data, and can be used to identify specific frequency components corresponding to dry years.

The short-time Fourier transform of the NDVI signal, X(t) is defined as:

$$X\left(t,f\right)={\displaystyle {\int }_{-\infty }^{\infty }x\left(t_1\right)}{h}^{*}\left(t_1-t\right)e^{-i2\pi f{t}_1}d{t}_1$$

(3)

where the function h(t) is a Hamming window centred at time t and multiplied by the signal X(t) prior to the Fourier transformation. The window function allows the signal to be viewed close to the time t only, and the Fourier transform provides a local estimation of the spectrum centred around time t. The width of this window was set to 90 days. Similarly to the ordinary Fourier transform and spectrum, the spectrogram can be written as:

$$S_x\left(t,f\right)={\left|X\left(t,f\right)\right|}^2$$

(4)

S_x can be expressed in decibels as:

$$P=10\log \left(S_x\left(t,f\right)\right)$$

(5)

Figure 9 shows the resulting WFT spectrograms (P values), for the annual agricultural cover produced from TERRA-MODIS time series. For this vegetation cover, a clear difference can be observed between dry and wet seasons, with the spectrogram revealing a greater intensity at high frequencies during wet years. Figure 10a shows the spectrograms corresponding to the relatively dry year from 2001 to 2002, and the wet year from 2010 to 2011, as a function of frequency. Higher levels of spectral power occur at high frequencies, during the wettest year (more than 4 dB). This observation is confirmed in Figure 10b, which plots the values of the spectrogram at 30 Hz (a high frequency). A similar behavior is observed for the NDVI temporal signal, with higher power occurring during wet seasons, and lower power occurring in the dry years, during which the dynamics of the vegetation cover are low. This leads to a relatively low intensity in the temporal variations of the NDVI signals. Conversely, during the wet seasons, NDVI variations resulting from fluctuations in the vegetation cover and rainfall events, etc. occur with a much higher intensity. Although in this section we present the results obtained with TERRA-MODIS products only, it is interesting to note that the results obtained with the SPOT-VGT products are very similar.

4.2. Monitoring of Bare Soils in Agricultural Areas

In the cultivated zones analyzed here, the bare soils represented a greater surface area during dry years, as a result of the farmers’ tendency to avoid sowing in the absence of rainfall events, during the period between November and January. The aim of the next section is to illustrate this behavior, by analyzing the NDVI time series.

In order to analyze the dynamics of bare soils, together with their relationship with drought events occurring during the studied period, we studied yearly remotely-sensed land use maps, recorded between 1998 and 2013, and broke them down into three classes (olive groves and bare soil (class 1), annual agriculture (class 2) and pastures (class 3)). In order to evaluate the dynamics of bare soils selected from class 1, during the course of wet and dry years, it is proposed to accurately determine the surface area covered by the olive groves only, corresponding to that identified during the wettest year (2010–2011). Under these climatic conditions, a very small percentage of bare soils remain uncultivated. Even in natural areas, pastures can be found everywhere. As this olive grove class has remained approximately constant in recent years, any increase in class 1 coverage during a dry year can be attributed to an increase in the surface area of bare soil.

Using this hypothesis, Figure 11 illustrates the temporal variation of surface areas covered by bare soils, determined from data recorded by the SPOT-VGT and TERRA-MODIS sensors. These two estimations are in good agreement, with the highest surface areas occurring during the driest years. As shown in Figure 11, only the TERRA-MODIS sensor was able to detect a strong increase in the surface area of bare soil during the dry years of 2000, 2001 and 2007. This can be explained by comparing the spatial resolutions of the two sensors: since that of SPOT-VGT is equal to 1 km, the probability of retrieving bare soil covering a complete pixel is lower than in the case of the TERRA-MODIS observations, which have a resolution of 250 m. These results clearly illustrate effects that are related to drought, but which may also be influenced by the farmers’ agricultural practices. In general, farmers do not sow cereals or vegetables whenever there are no significant rainfall events at the beginning of the agricultural season (end of autumn, beginning of winter). This result is confirmed by the correlation between precipitation and the surface area covered by bare soil, during the period between November and January. The correlations between bare soil surface area and precipitation levels are equal to 0.36 and 0.8, for the SPOT-VGT and TERRA-MODIS data, respectively. For this reason, the influence of drought on vegetation cover could not be accurately identified in these agricultural areas, since the farmers’ agricultural practices were not taken into account.

5. Conclusions

This paper proposes to interpret SPOT-VGT and TERRA-MODIS NDVI time series for the purposes of analyzing the behavior of vegetation as a consequence of drought. In the first step, land use mapping is implemented using data from both of these satellites. Three agricultural classes are considered: olive groves and bare soil, pastures, and annual agricultural areas. The results obtained from SPOT-VGT and TERRA-MODIS observations are found to be qualitatively coherent in the case of the latter class of agricultural land use. These results are validated through the use of high-resolution SPOT/HRV land use maps recorded over a period of eight years. For the total surface area covered by annual agricultural fields, a better agreement is found between SPOT-HRV and TERRA-MODIS data than between SPOT-HRV and SPOT-VGT data, with correlation coefficients equal to 0.89 and 0.59, respectively. VAI (Vegetation Anomaly Index) estimations are proposed for a period of more than 13 years, through the use of SPOT-VGT and TERRA-MODIS NDVI products. The VAI estimations are found to be coherent, despite differences in spatial resolution between these two sensors, and a linear relationship is established on a monthly basis, between the two sets of VAI estimations. A positive trend is observed for the VAI index determined with these two products. For the case of annual agricultural areas, a Windowed Fourier Transform spectrogram is computed from the NDVI time series, to simplify the analysis of the temporal variations of this product, and to determine when certain cyclic events occur. Drought years are qualitatively characterized by spectrograms having low levels of high frequency spectral power. The dynamics of the VAI and NDVI characterizing of agricultural areas are strongly related to variations in bare soil surface coverage. This relationship was derived from the inter-annual dynamics of land use mapping. Using the two satellites, SPOT-VGT and TERRA-MODIS, we are able to show that an increase in bare soil coverage occurs during dry years. This outcome is related mainly to the decision, made by most farmers, not to sow new crops during agricultural years that are preceded by an autumn with low levels of precipitation. This aspect is clearly revealed by the TERRA-MODIS data, due to its higher spatial resolution and greater sensitivity for the detection of fields with bare soils.