A Dynamic Vegetation Senescence Indicator for Near-Real-Time Desert Locust Habitat Monitoring with MODIS

Abstract

Desert locusts (Schistocerca gregaria) represent a major threat for agro-pastoral resources and food security over almost 30 million km² from northern Africa to the Arabian peninsula and India. Given the differential food preferences of this insect pest and the extent and remoteness of the their distribution area, near-real-time remotely-sensed information on potential habitats support control operations by narrowing down field surveys to areas favorable for their development and prone to gregarization and outbreaks. The development of dynamic greenness maps, which detect the onset of photosynthetic vegetation, allowed national control centers to identify potential habitats to survey, as locusts prefer green and fresh vegetation. Their successful integration into the daily control operations led to a new need: the near-real-time identification of the onset of dryness, a synonym for the loss of habitat attractiveness, likely to be abandoned by locusts. The timely availability of this information would enable control centers to focus their surveys on areas more prone to gregarization, leading to more efficiency in the allocation of resources and in decision making. In this context, this work developed an original method to detect in near-real-time the onset of vegetation senescence. The design of the detection relies on the temporal behavior of two indices: the Normalized Difference Vegetation Index, depending on the green vegetation, and the Normalized Difference Tillage Index, sensitive to both green and dry vegetation. The method is demonstrated in Mauritania, an ever-affected country, with 10-day MODIS mean composites for the years 2010 and 2011. The discrimination performance of three classes (“growth”, “density reduction” and “drying”) were analyzed for three classification methods: maximum likelihood (61.4% of overall accuracy), decision tree (71.5%) and support vector machine (72.3%). The classification accuracy is heterogeneous in both time and space and is affected by several factors, such as vegetation density, the north-south climatic gradient and the relief. Smoothing the vegetation time series resulted in an increase of the overall accuracy of about 5% at the expense of a loss in timeliness of ten days. To simulate near-real-time monitoring conditions, the decision tree was applied to the decade of 2010. Overall, the seasonal vegetation cycle appeared clear and consistent. The results obtained pave the way for an operational implementation of the senescence dynamic mapping and, consequently, to further strengthen the capacity of the locust control management.

1. Introduction

Desert locust (Schistocerca gregaria (Forksål, 1775)) is a well-known pest that represents a major threat for agro-pastoral resources, because of its voracity and ability to travel long distances [1,2]. Locust infestations and invasions jeopardize repetitively the food security in more than 60 countries from northern Africa to the Arabian peninsula and covering over 30 million km². During invasion periods, locusts swarms count millions of individuals capable of destroying hectares of crops in only a few hours [3]. The 2003–2005 invasion in West Africa hit more than eight million people and destroyed 80%–100% of the expected harvest, leading to a dramatic increase of food insecurity in the affected countries. Control operations, including the treatment of more than 13 millions hectares, reached a total cost of 400–500 million U.S. dollars [4]. During the 1987–1989 invasion, the school enrollment rate in Malian affected villages dropped by 25%, falling under 18% [5].

Desert locusts exhibit a phase polyphenism in which behavioral, physiological and morphological traits change due to variations in the local population density [6]. Three mutually-interdependent and concurrent phenomena lead to that transformation [7,8]: multiplication, population concentration and gregarization. Under these conditions, solitarious phase individuals shift to gregarious phase individuals, which are more voracious and mobile. The gregarization density threshold varies according to development stages and is estimated at 250–500 adults locusts per hectare [1].

Desert locusts have clear food preferences: some plants are consumed before others, and some species remain untouched [9–12]. Food plant preferences vary from one developmental stage to another [10,13] and depend on plant association [9]. This differential feeding directly influences the behavior of locusts, e.g., for habitat selection and marching speed [14,15]. For instance, locust bands move faster through shrubby habitats consisting of moderately palatable perennials than they do in sand dune communities, where locusts feed extensively on the highly palatable annuals [9].

To protect human populations from the disastrous consequences of locust invasions, the Food and Agriculture Organisation (FAO) promotes a preventative control strategy and coordinates national locust control centers that regionally operate surveys and control actions [16–18]. That preventative control strategy aims at avoiding the development of generalized invasions and is ensured by three essential activities [1,19,20]: (i) monitoring the meteorological and environmental conditions of outbreak areas by means of meteorological data and satellite imagery; (ii) surveying in areas potentially favorable for locusts; and (iii) carrying out control measures and treatment in areas that present a locust density above the critical density threshold.

Remote sensing plays a major part in this strategy, as it allows monitoring the locust habitats in near-real time and over its entire distribution area, especially in remote or unsecured regions. Given locusts’ differential food preferences and the extent and remoteness of the desert locust distribution area, this near-real-time satellite-derived information on the locust potential habitat supports control operations by narrowing down field surveys to areas prone to gregarization and outbreaks. Research has long been directed toward mapping locusts’ habitats with satellite remote sensing imagery, as the current spatial resolution of satellites remains insufficient for direct detection of locusts. Numerous studies tackled the issues of monitoring the desert locust breeding areas visually and by means of digital image processing [21–23]. In particular, the Normalized Difference Vegetation Index (NDVI) has been intensively used to link vegetation with locust occurrence. For instance, [24] adopted a visual approach to investigate the relationship of high NDVI areas with locust activity. Later, [22,25] demonstrated that an index based on the NDVI was correlated with rainfall and locust infestations. [26] concluded that the occurrences of high NDVI values correspond to the emergence of vegetation, which can be further exploited to monitor the reproduction areas of locusts. [27] related abundance and the spatial distribution of resources at the landscape scale to historical records of outbreaks thanks to satellite data. They also showed that the responses were different for vegetation growing near the Red Sea compared to Mauritania. [28] suggested an NDVI threshold of 0.14 to discriminate vegetation from bare soil. This threshold was applied on a mosaic of SPOT-VGT images covering the entire desert locust distribution area. This technique yielded an acceptable result with the detection of most homogeneous events. However, in certain cases with sparse vegetation, favorable areas are omitted. In order to reduce those omissions, 250-m MODIS images have been coupled to 1-km SPOT-VEGETATION images [29]. This combination allows reducing the omission errors thanks to MODIS’ higher spatial resolution. Furthermore, commission errors were reduced thanks to the MODIS SWIR bands. Since 2010, dynamic greenness maps (i.e., the green area, developed by [30]) are provided every ten days to the Food and Agriculture Organization and the affected countries to help them analyze the current situation and plan the control operations [30–32]. These maps highlight vegetated areas close to the onset of greenness, which correspond to green and fresh vegetation that is highly palatable to desert locust. These areas highlight priority areas for survey (i.e., locust potential habitats), increasing consequently the efficiency and the rapidity of field operations, while simultaneously reducing their costs.

Conversely, detecting the areas where conditions are becoming less favorable for locusts would also substantially increase the survey efficiency. Hence, mapping the senescence of vegetation would highlight the areas likely to be abandoned by locusts, as they prefer fresh and green vegetation. Such dryness maps would enable control centers to assess if they can concentrate their efforts in priority areas. Since recently, few studies that have addressed the detection of dry vegetation by remote sensing have tested monitoring methods efficient at a large scale [33,34], since field or airborne spectroscopy or high-resolution data from Landsat [35–37] or Hyperion [38–41] were used in most cases. [33,34] have shown the value of MODIS SWIR bands to monitor dry vegetation in a representative panel of Sahelian region. These recent studies laid the foundations for a better locust monitoring using remote sensing

In order to support a cost-effective and efficient strategy, the objective of this study is to develop a methodology, based on the latest advances in the field, for mapping the senescence of vegetation in arid and semi-arid areas compliant with the near-real-time requirements of the application. The tested hypothesis is that a reliable discrimination of the onset of vegetation senescence can be achieved thanks to the joint temporal trajectory of two spectral indices: the Normalized Difference Vegetation Index (NDVI), sensitive to green vegetation, and the Normalized Difference Tillage Index (NDTI), sensitive to both green and dry vegetation. Special emphasis was put on the study of the temporal and spatial variability of the accuracy of the resulting dryness maps. Finally, a case study demonstrates how the proposed method can support desert locust habitat monitoring by simulating a near-real-time monitoring, i.e., as soon as images are available, of an area in Mauritania during the 2010 growing season.

2. Material

2.1. Study Site

This study focuses on the Islamic Republic of Mauritania, a front-line country spreading over more than a million square kilometers between latitudes 14;50;N and 27;20;N and longitudes 5;40;W and 16;5;W (Figure 1). As a front-line country, its vast and diverse territory is a near-permanent habitat for desert locusts with complementary biotopes along the seasons. Sixty five to ninety percent of the Mauritanian active population is employed by the primary sector; therefore, the locust threat is high [2].

Sahelian and Saharan vegetation is dominated by annual plants, associated with sparse stands of woody plants. Annual plants germinate after the first rainfall (July) of the rainy season and complete their reproductive cycle rapidly, with fruiting generally occurring at the end of the rainy season (September). Then, annual vegetation withers to give litter or straw, which are grazing reserves for the dry season (a period of 8–9 months) [2,42]. During the dry season, herbaceous vegetation mass decreases by 50%–80%, even without being subjected to grazing, as a result of microbiological degradation [43].

Rainfall plays a key role for desert locusts: (1) it determines the edaphic humidity, which is essential to egg survival and development; and (2) it conditions the growth of vegetation that provides both shelter and food. Therefore, the locust cyclical dynamic is tightly tied to that of vegetation, which itself depends on rainfall. In Mauritania, rainfall decreases from south to north and from west to east, because of the Sahara-Sahel eco-climatic gradient and the increasing continental influences, respectively. The rainy season starts in July in the southern part of the country and spreads north. Locust populations over one year are first located mostly in the summer breeding area of the south. When vegetation dries out, they migrate to the north, where vegetation is still green, to complete their development [2].

2.2. Data Sources

In total, three kinds of data were utilized: field data, remote sensing time series and satellite-derived greenness maps.

The RAMSES database (Recognition and Monitoring System of the Environment of Schistocerca) developed by FAO compiles field observations collected by locust survey teams since 1988. Each geo-/chrono-referenced observation concerns: (i) the locust information; (ii) the meteorology; and (iii) the floristic information (

Table 1. RAMSES database fields and corresponding values.

Attribute	Possible Values
Date	11 January 2009–20 June 2011
Phenological stage	Sprout/Green/Greening/Drying/Dry
Total cover	Low/Moderate/Dense
Annual crop cover	0%–100%
Perennial crop cover	0%–100%
Prospected surface	0–100,000 hectares
Infested surface	0–100,000 hectares
Control action	Yes/No
Habitat	Wadi/Interdune/Plains/Basin

). The vegetation development stage and the percentage of coverage and greenness are evaluated for annual and perennial plants, first separately and then jointly. The vegetation development stage is established following five qualitative classes: sprout, growth, green, drying and dry. Annual and perennial species are listed in order of decreasing dominance [2]. The accuracy of the RAMSES database was evaluated at 88% [2], even if some inconsistencies have been documented [44].

Operational monitoring in the Sahelian and Saharan regions requires Earth observation data with large spatial coverage, a high revisit frequency and free access. With a spatial resolution of 250-m for the red in near-infra-red channels and 500-m for the short wave infra-red channels and a daily revisit frequency, the MODIS instrument aboard the Terra and Aqua satellites appears as an appropriate solution to cope with those requirements. Daily MOD09GA and MYD09GA products were downloaded for the period from January 2009 to April 2011.

For the corresponding period, the green area product [30] was downloaded. This product provides dynamic greenness maps for near-real-time operational monitoring of the desert locust habitat. The methodology behind the green area is based on a colorimetric transformation (R, NIR, SWIR) from the RGB color space to HSV (hue, saturation, value). In this color space, H appears to be a qualitative spectral index [30], which decouples information on land cover. Its temporal variation can be interpreted as a land cover change. Changes in S and V may be interpreted as changes in observation conditions (atmospheric effects, geometry of acquisition) and changes in vegetation (phenological stage, floristic composition). A set of thresholds were defined by statistical sampling to distinguish vegetation from non-vegetation on 10-day MODIS mean composite images [45]. The final product integrates a temporal component thanks to a time meter that accounts for the number of 10-day intervals a pixel has been tagged as vegetation. A color table allows an intuitive and straightforward identification of ephemeral (red), seasonal (light green) and perennial vegetation (dark green). This provides operationally and in near-real time a user-friendly vegetation dynamic map updated every 10 days over the whole desert locust area.

3. Methodology

3.1. Rationale

Several studies investigated dry vegetation detection from satellite imagery. Different objectives were pursued such as assessment of soil tillage intensity and soil conservation [38,46,47], evaluation of soil erosion risk and runoff [48–51], evaluation of the risk of wildfire in relation to dead fuel proportion [41,52,53], and improvement in land cover mapping [39,54]. Recently the study of Jacques et al. [33], demonstrated the value of SWIR (Short Wave Infra Red) MODIS bands to monitor dry vegetation in a Sahelian region. In particular, it has been shown that combinations (simple ratio and normalized difference ratio) of MODIS band 6 (1.628–1.652 µm) and 7 (2.105–2.155 µm) are sensitive to both dry and green vegetation cover and masses that cover brighter and rougher soil background. One of those indices is a MODIS adaptation (Equation (1)) of the Normalized Difference Tillage Index (NDTI) from Van Deventer et al. [55], initially developed for Landsat.

$$NDTI=\frac{{\rho }_{MOD6}-{\rho }_{MOD7}}{{\rho }_{MOD6}+{\rho }_{MOD7}}$$

(1)

where ρ_MOD6, ρ_MOD7 correspond to MODIS band 6 (1.628–1.652 µm), MODIS band 7 (2.105–2.155 µm) respectively.The effect of water absorption is strong in the SWIR range and is more important for the band 7 (2200 m⁻¹) than for band 6 (498 m⁻¹), which could potentially affect the NDTI. However in semi-arid area, the water content of the vegetation rapidly decreases during the growing season to a value (<40% at the peak of biomass) that has a low impact on the index. Special attention has to be paid over wet soils and during the early growing season (water content close to 80%).

Along the growing season, when vegetation is still green, NDVI and NDTI are linearly related, because they share similar sensitivity to green vegetation [33]. The relationship changes when the vegetation starts drying. The NDTI curve slowly decreases during senescence (as green vegetation dries, it remains detected by the NDTI). Comparatively, the NDVI decreases rapidly. Indeed, while green vegetation dries out, it is no longer detected by the NDVI and only sensitive to photosynthetically-active vegetation (Figure 2a). The moment at which a separation between the NDVI and NDTI curves appears, i.e., the stall point, characterizes the onset of senescence. In addition, the vegetation cover density diminishes due to the combined effects of grazing, wind erosion and microbiological degradation. When significant, NDVI and NDTI decrease similarly (Figure 2b). In this case, the indices do not allow a direct identification of the onset of the senescence. However, this description of the joint temporal behavior of NDVI and NDTI is not always observed. Vegetation indices are influenced by several other factors than the photosynthetic state, such as soil type and humidity, relief and the geometry of observation (view and Sun angles), which add noise into the relationship. Furthermore, these effects do not affect equally the NIR and the SWIR wavelengths, which makes a proper intercomparison of the two indices much more difficult. To overcome some of these constraints, better than considering the difference between the two indices, the slope of their time series could be a pertinent alternative. During the growing season, both indices increase. The slope of the time series is positive. During the dry season, the NDVI signal drops more rapidly than that of NDTI if the vegetation density decreases slowly. Therefore, drying could be identified by a higher negative slope observed for the NDVI than for NDTI. A negative slope similar for the two indices would be more difficult to interpret and would mean a density loss without certainty of drying.

For locusts, vegetation density reduction is an important factor that can foster gregarization by concentrating the population on small surfaces [1,27,56]. Therefore, information on density reduction remains valuable, even if the actual state of the vegetation (green or dry) is unknown. Areas identified as “drying” are associated with a low locust threat, from which the survey teams can be withdrawn. However, areas described as “growth” or “density reduction” indicate a higher locust threat because of either the presence of green vegetation (i.e., locust potential habitats) or risk of gregarization by population concentration.

3.2. Data Preprocessing

3.2.1. Field Observation Preprocessing

Prior further analysis, the RAMSES observations needed to be converted from its specific nomenclature into the three classes of interest, namely “growth”, “density reduction” and “drying”. First, the five initial vegetation development classes were merged into (i) growth (sprout, growth and green) and (ii) drying (drying and dry). Second, the samples (11,441 from May 2010–June 2011) were related to their corresponding dynamic greenness map value, i.e., the number of 10-day intervals during which the corresponding pixel is detected as vegetation. To ensure consistency with the green area product, the analysis focused on pixels of vegetation potentially entering senescence: only pixels occurring as fresh and green vegetation, during the current 10-day interval and the three previous were kept. Third, the vegetation index time series of the remaining samples (1349 observations) and the corresponding phenological stage of the development were analyzed by an expert. Each sample was classified by visual interpretation in one of the three following classes according to the joint temporal behavior of their NDVI and NDTI curves:

• Growth: the observation was located on the increasing part of the NDVI curve and was designated as “green” in the database RAMSES;
• Drying: the observation was reported as “drying” in RAMSES and was located on the decreasing part of the NDVI curve while the NDTI curve remained stable (i.e., presence of a stall point). There was no apparent decrease of the vegetation cover and the vegetation shifted from a “green” state to a “drying” state. The NDTI remained constant before entering a decrease phase due to the vegetation degradation;
• Density reduction: the observation was described as “green” or “drying” in RAMSES and was on the decreasing part of the NDVI and NDTI curves. As NDTI is sensitive to both green and dry vegetation, a decrease in its temporal trajectory indicates a reduction of vegetation. This class does not exclude the two others, but in this case, the actual state of the vegetation remained undetermined (likely drying). Caution should be taken in the interpretation of the errors of the classification of this class.

Only 94 points were found unclassifiable (absence of seasonal cycle or too noisy) and set aside. Besides, 126 observations were removed from the database, as they showed contradictory information: NDVI and NDTI indicated a “drying” status, whereas RAMSES listed them as “green”. These points can result from survey, encoding or geolocation errors. The final database contained 1129 samples (532 in “growth”, 308 in “density reduction” and 289 in “drying”) well distributed over the study area (Figure 1) and available for further calibration and validation purposes.

3.2.2. Remote Sensing Data Preprocessing

In order to avoid punctual artifacts and improve temporal consistency, 10-days mean composites have been produced from daily Terra and Aqua images according to the procedure detailed in Vancutsem et al. [45]. All of the quality-controlled reflectance values (without clouds and hazes) were averaged for each pixel and each used band (1,2,6,7) over a 10-day period. NDVI and NDTI were computed from these averaged reflectance values.

The noise inherent to the remotely-sensed signal may negatively affect the computation of indices and temporal metrics (slope computation in particular), decreasing de facto the classification accuracy. Smoothing may avoid such spurious effects and enhanced class separability. Smoothing methods have been intensively devised in the literature (e.g., Fourier analysis, asymmetric Gaussian model, double logistic model). In particular, Atkinson et al. [57] highlighted the performance of the Whittaker filter for smoothing of satellite image time series. The non-parametric Whittaker method fits a smoothed curve to discrete values by finding a trade-off between fidelity to the original data and the regularity of the smoothed time series by penalizing the roughness of the smoothed curve [58].

3.3. Senescence Detection and Dynamic Mapping

The overall methodology consists of two steps: metric selection and classifier benchmarking: first, in defining vegetation indices’ temporal metrics, from which the most fitted to the classification purpose were selected. Feature selection is limited in order to keep the method simple, reproducible and computationally economical. Furthermore, using few features allows one to link their spectral behaviors with the reality on the ground and, thus, benefit from expert knowledge to interpret the results. Three classification methods were then trained on the selected features and benchmarked to discriminate the three phenological classes “growth”, “density reduction” and “drying”. After the assessment of their respective accuracy, the best performing method was applied on a subset of Mauritania to produce 10-day interval vegetation status maps. A temporal combination of those provides dynamic dryness maps.

3.3.1. Metric Definition and Selection

Eleven combinations of NDVI and NDTI, hereafter referred to as metrics, have been tested (

Table 2. Definition of the metrics tested.

Metrics	Definition
NDVI − NDTI	difference between NDVI − NDTI
$\mathrm{\Delta }NDV{I}_t^{t-1}$	NDVI slope over one past 10-day interval
$\mathrm{\Delta }NDT{I}_t^{t-1}$	NDTI slope over one past 10-day interval
$\mathrm{\Delta }NDV{I}_t^{t-2}$	NDVI slope over two past 10-day intervals
$\mathrm{\Delta }NDT{I}_t^{t-2}$	NDTI slope over two past 10-day intervals
$\mathrm{\Delta }NDV{I}_{t+1}^{t-1}$	NDVI slope over the past 10-day interval to the next 10-day interval
$\mathrm{\Delta }NDT{I}_{t+1}^{t-1}$	NDTI slope over the past 10-day interval to the next 10-day interval
$\mathrm{\Delta }NDV{I}_t^{t-1}+\mathrm{\Delta }NDV{I}_t^{t-2}$	the cumulative sum of the NDVI slopes over each two past 10-day intervals
$\mathrm{\Delta }NDT{I}_t^{t-1}+\mathrm{\Delta }NDT{I}_t^{t-2}$	the cumulative sum of the NDTI slopes over each two past 10-day intervals
$\mathrm{\Delta }NDV{I}_t^{t-2}-\mathrm{\Delta }NDT{I}_t^{t-2}$	difference of NDVI and NDTI slope over two past 10-day intervals
$\mathrm{\Delta }NDV{I}_t^{t-2}\pm \mathrm{\Delta }NDT{I}_t^{t-2}$	sum of NDVI and NDTI slope over two past 10-day intervals

): (i) the difference between the two indices NDVI−NDTI; (ii) the slope over one past 10-day interval $\mathrm{\Delta }Inde{x}_t^{t-1}$; (iii) two past 10-day intervals $\mathrm{\Delta }Inde{x}_t^{t-2}$; (iv) the past 10-day interval to the next 10-day interval $\mathrm{\Delta }Inde{x}_{t+1}^{t-1}$; (v) the cumulative sum of the slopes over each two past 10-day intervals $\mathrm{\Delta }Inde{x}_t^{t-1}+\mathrm{\Delta }Inde{x}_t^{t-2}$; and (vi) the sum and the difference of slope over two past 10-day intervals of the two indices $\mathrm{\Delta }NDV{I}_t^{t-2}\pm \mathrm{\Delta }NDT{I}_t^{t-2}$. The cumulative sum over several 10-day intervals has also been assessed, as it should help strengthen the signal and thus contribute to better detection.

The performance of the eleven metrics was analyzed in a hierarchical fashion by two sets of receiver operating characteristic (ROC) curves. This method has long been recognized as a convenient and valuable tool for algorithms evaluation and comparison, but also for its specific properties, which handle distribution asymmetry [59]. The first set is for the discrimination of classes “growth” and “decrease” and the second for the distinction of classes “density reduction” and “drying” (Figure 3). The first set was derived from the preprocessed RAMSESdata combining the classes “density reduction” and “drying” into a single category, hereafter referred to as “decrease”. The second set of ROC curves included only the samples from the “decrease” category to further assess the metrics’ ability to discriminate “density reduction” from “drying”. This hierarchical approach is necessary due to the binary character of the ROC curve method.

In the present case, ROC curves allowed one to select metrics based on their classification performances. Two accuracy figures were compared: the sensitivity (or success rate) and the specificity. The sensitivity corresponds to the probability for a pixel of class ω₁ to be well classified by the algorithm (i.e., producer’s accuracy of class ω₁). The specificity provides the probability of a pixel of class ω₂ to be included in class ω₂ (i.e., producer’s accuracy of class ω₂). Specificity is equal to 1 − fp, where fp is the false alarm rate of class ω₁ (i.e., the probability of a ω₂ pixel to be incorrectly assigned to class ω₁) [59]. The ROC curve represents the trade-off between these two accuracy measures for different metrics thresholds. The area under the curve (AUC) is among the most commonly-chosen [60] ROC curve accuracy measure for method comparison. The AUC ranges from 0%–100%, 100% representing an error-free classification. As a random classification yields an AUC of 50%, no realistic classification should have an inferior AUC [59]. Overall accuracies corresponding to the highest AUC were also derived for each metric (

Table 3. Maximum accuracy parameters for the classification using the selected metrics for raw and smoothed data (in percent). The notation “growth-decrease” refers to the discrimination of the categories “growth” and “decrease” (Figure 3a) and “density-drying” to the discrimination of the classes “density reduction” and “drying” (Figure 3b).

Metrics	Raw Data				Smoothed Data
	Growth-Decrease		Density-Drying		Growth-Decrease		Density-Drying
	AUC	OA	AUC	OA	AUC	OA	AUC	OA
$\mathrm{\Delta }NDV{I}_t^{t-1}$	88.38	82.70	59.35	59.09	93.81	87.78	52.94	57.96
$\mathrm{\Delta }NDV{I}_t^{t-2}$	92.50	87.04	55.73	60.78	93.21	88.04	48.83	55.61
$\mathrm{\Delta }NDV{I}_t^{t-1}+\mathrm{\Delta }NDV{I}_t^{t-2}$	93.11	87.72	55.63	60.92	93.87	88.40	50.26	55.95
$\mathrm{\Delta }NDV{I}_{t+1}^{t-1}$	87.59	83.13	64.17	65.15	96.07	91.39	56.81	59.52
$\mathrm{\Delta }NDT{I}_t^{t-2}$	70.14	69.84	69.68	66.16	73.18	70.77	71.55	67.67
$\mathrm{\Delta }NDT{I}_t^{t-2}$	74.82	73.82	64.98	64.86	78.63	75.20	62.22	64.49
$\mathrm{\Delta }NDT{I}_t^{t-1}+\mathrm{\Delta }NDT{I}_t^{t-2}$	73.48	72.87	70.71	66.38	78.43	73.87	68.12	66.83
$\mathrm{\Delta }NDT{I}_{t+1}^{t-1}$	76.25	72.45	72.80	70.37	81.45	78.12	75.12	71.36
NDVI − NDTI	70.32	64.66	66.99	51.59	86.18	80.69	51.20	51.59
$\mathrm{\Delta }NDV{I}_t^{t-2}+\mathrm{\Delta }NDT{I}_t^{t-2}$	85.04	78.44	52.63	54.33	91.43	85.12	54.73	59.13
$\mathrm{\Delta }NDV{I}_t^{t-2}+\mathrm{\Delta }NDT{I}_t^{t-2}$	89.22	85.27	61.56	64.35	66.29	62.09	66.35	51.59

).

3.3.2. Classification Methods and Accuracy Assessment

Based on the two best features, selected thanks to the ROC curve accuracy assessment, three classification algorithms associated with various degrees of complexity have been benchmarked for the discrimination of the three classes of interest (“growth”,“density reduction” and “drying”): a decision tree (DT), a support vector machine (SVM) and a maximum likelihood (ML) classifier. The ML method is a parametric supervised classification derived from Bayes’ theorem. It uses the frequency histogram of each class directly as the discriminant function. A given pixel x is assigned to the class whose probability a posteriori (P(ω_i|x)) is the highest [61]. The decision threshold between ω_i and ω_j classes is defined by g_i (x) = g_j (x) with x ϵ ω_i if g_i (x) > g_j (x) for all j ≠ i [62]. The underlying assumption of this method is that the distribution of the vectors x values for each class is Gaussian and can be determined only by the average µ_i and the standard deviation σ_i of the considered class. Unlike other classification methods that use a set of measures to establish the decision rule in a single step, DT relies on a hierarchical approach with several steps. Decision tree operates on a so-called top-down approach: it begins with a primary node, containing all data, that is then divided to form all intermediate nodes downstream, which lead to terminal nodes. Each node corresponds to a binary decision rule that separates either a single class or several classes from all remaining classes [63]. Furthermore, the decision tree has other advantages, such as no assumption about the data distribution in each class and rapid training and implementation [64,65]. SVM belongs to the category of machine learning algorithms and allows discrimination of different classes in a multidimensional space. The basic principle is to locate the optimal boundaries between classes by maximizing the margin between the separating hyperplanes. Initially, the SVM was designed for binary classification, but has been expanded to multi-class classification [66]. In addition, the SVM has a low sensitivity to the training sample size and shows a capacity to handle high dimensional spaces [67].

Using two thirds of the preprocessed field data, each classifier was trained twice: the first time with metrics directly computed on the signal and second time on metrics derived from smoothed time series. The performance of the six classification (three algorithms and two sets of metrics) was assessed by means of the error matrix and derived statistics, such as: the overall accuracy (OA), the Kappa coefficient (κ), the omission and commission errors (OE and CE) and the F₁-score. In the framework of desert locust habitat monitoring, [32] showed that accuracy varies from one region to another and also along seasons within the same region. Therefore, the spatial, temporal and thematic variability of the classification accuracy has also been analyzed: (1) the RAMSES points were gathered by macro-region (i.e., homogeneous ecological units defined by [2]), and the global classification error (1 − OA) by macro-region has been represented spatially; (2) the RAMSES points were grouped by month, and the global error for each month was computed; and (3) the classification errors were calculated for each three vegetation density categories of the RAMSES database (“low” (25 points), “middle” (625 points), “dense” (368 points)).

3.4. Near-Real-Time Simulation Study Case

To assess the near-real-time capabilities of the proposed method, the most accurate classification method was then applied to a region of Mauritania from July 2010–June 2011 on a 10-day basis. To mimic the operational context, metrics are computed along the season with all images that would be available at the date of interest. To ensure consistency with the green area product, the method was only applied to pixels flagged as vegetation by the green area, as only vegetation can dry. One of the strengths of the dynamic greenness maps of the green area was its color code, allowing users to rapidly interpret the vegetation status and relate it to locust development [32]. Similarly, a color code was established for the dry season thanks to a simple time meter accumulating the detections for 1, 2, 3 and 4 or more 10-day intervals of each class. To the three classes discriminated by the algorithm (growth, density reduction and drying), the “dry” class was added for pixels not detected as vegetation anymore. In practice, this color code allows the interpreter to identify quickly, on the one hand, the areas from where the locusts are flying away (“dry” and “drying”) and, on the other hand, the areas still to be prospected, as density reduction might favor gregarization (“density reduction”).

4. Results and Discussion

4.1. Metric Evaluation

Regarding the discrimination of the “growth” and “decrease” classes (Figure 3a), the accumulated slopes $\mathrm{\Delta }NDV{I}_t^{t-1}+\mathrm{\Delta }NDV{I}_t^{t-2}$ present the maximum accuracy with an area under the curve of 93.11% and an overall accuracy of 87.72%. Discrimination is strengthened when the slope is computed over a longer period and then accumulated as it smooths the local variability (Figure 4a). Other metrics have an AUC ranging from 70% to 93% and an OA of 64%–87%. NDVI features (AUC from 88% to 93%), systematically outperforming their NDTI counterparts (70%–76%). Features combining the two features have an AUC in between these ranges.

Regarding the distinction of the “density reduction” and “drying” classes (Figure 3b), the slope of NDTI calculated with a delay of one 10-day interval $(\mathrm{\Delta }NDT{I}_{t+1}^{t-1})$ shows the highest accuracy measures. With an AUC of 72.80% and an OA of 70.37%, these values are lower than those found $\mathrm{\Delta }NDT{I}_t^{t-1}+\mathrm{\Delta }NDT{I}_t^{t-2}$ for the first step of discrimination. Despite and $\mathrm{\Delta }NDT{I}_{t+1}^{t-1}$, all metrics discriminate the first level more accurately than the second, highlighting the difficulty in the classes “density reduction” and “drying”. Even with their coarser spatial resolution, NDTI features consistently discriminate better than the NDVI features (OA of 64%–70% and 55%–59%, respectively). The difference between NDVI and NDTI yields poor result. This confirms that noise is added due to different factors, such as soil type and humidity, relief and the geometry of observation (view and Sun angles). As the “density reduction” class is not exclusive, misclassification errors between the “density reduction” and the “drying” classes might be “false alarm”, as the vegetation may still dry (or less likely, remain green) and see its density reduced simultaneously. Again, these NDTI and NDVI metrics sandwich the accuracy of the metrics that combine them. This tends to confirm the rationale of the study that postulates that the joint use of NDVI and NDTI enhance the recognition capabilities of vegetation dynamics.

With smoothing, the results of the ROC curves displayed an improvement in the discrimination of classes, especially for the first level of classification. The slope of NDVI with a delay of one 10-day interval $(\mathrm{\Delta }NDV{I}_{t+1}^{t-1})$ shows the best result for the distinction of categories “growth” and “decrease” (AUC of 96.07% and OA of 91.39%), while the slope of NDTI with a delay of one 10-day interval $(\mathrm{\Delta }NDT{I}_{t+1}^{t-1})$ separates best the classes “density reduction” and “drying” (AUC = 75.12% and OA = 71.36%). In particular, smoothing the data increases the accuracy of slope indices computed using data following the date of interest $(\mathrm{\Delta }NDV{I}_{t+1}^{t-1})$ by 10%. This can be partly explained by the reduced influence of the signal drop following the growth peak due to the smoothing.

4.2. Classification Benchmarking

The performance of the three classification methods was evaluated by global and class-specific accuracy measures (

Table 4. Accuracy parameters of the three classifiers with original and smoothed data. DT, decision tree.

Classifiers	Omission Errors			Commission Errors			F1-score			OA	κ
	Growth	Dens.	Sen.	Growth	Dens.	Sen.	Growth	Dens.	Sen.
DT	16	46	29	13	33	42	85	59	63	71	56
SVM	14	49	25	15	30	41	85	58	65	72	57
ML	38	62	15	10	40	53	72	46	60	61	42
DTSmoothing	11	36	31	7	32	40	90	65	63	76	63
SVMSmoothing	9	43	27	10	29	38	89	62	66	76	63
MLSmoothing	31	51	13	5	27	54	79	57	59	67	51

). The classifiers were trained with the RAMSES data and with the two best performing metrics of the unsmoothed ( $\mathrm{\Delta }NDV{I}_t^{t-1}+\mathrm{\Delta }NDV{I}_t^{t-2}$ and $\mathrm{\Delta }NDT{I}_{t+1}^{t-1}$) and smoothed case ( $\mathrm{\Delta }NDV{I}_{t+1}^{t-1}$ and $\mathrm{\Delta }NDT{I}_{t+1}^{t-1}$). The ML had consistently the poorest performance, while SVM yielded the best overall results with an overall accuracy of 72.33% and a Kappa of 57.74%. The difference of SVM with the DT was marginal: only 1% lower for the OA and the Kappa coefficient. The standard deviation of the overall accuracy was estimated at 2% for all three methods. Overall, the kappavalues range between 42% and 57%, which qualifies the accuracy as moderate according to Landis’ and Koch’s scale.

The omission and commission errors, as well as the F₁-score reveal all the complexity to discriminate the density reduction from the dry vegetation. The highest omission errors are obtained for the “density reduction” (46% for the DT and 49% for the SVM), while “drying” shows the highest commission error (42% for the DT and 41% for the SVM). Samples belonging to “density reduction” are often omitted and erroneously attributed to the class “drying”, which contributes to its commission error. In comparison, omission and commission errors of the “growth” class remain limited, with 16% and 13% for the DT and 14% and 15% for the SVM, respectively. The F₁-score reflects these trends and is therefore higher for “growth” (85% for the DT and 85% for the SVM) and then lower for “drying” (63% and 65%) and, finally, for “density reduction” (59% and 58%).

Smoothing the time series prior to metric extraction was found to improve the accuracy (Table 4): the overall accuracy increases by 6% for the ML (OA of 67%, σ_OA of = 1.8%), by 5% for the DT (OA of 76%, σ_OA of 1.9%) and by 4% for the SVM (OA of 76% σ_OA of 1.8%). However, smoothing might appear as a constraint in an operational context, as smoothed metrics of 10-day interval t will require the image of 10-day interval t + 1. Thus, smoothing delays the timing of monitoring by 10 days, and the optimal trade-off between timeliness and accuracy that meets the operational needs has yet to be identified.

All errors do not affect desert locust habitat monitoring similarly. Omission errors for drying vegetation (25%–29%) result in a loss of efficiency, as surveys in those unfavorable areas are not required. On the contrary, commission errors for this class (41%–42%) impact operations more negatively, as those areas could still be potential habitat for locusts. Errors for the “density reduction” class strongly affect the operations, as these areas still require surveying, because vegetation might still be green. When vegetation contracts into smaller patches, migrating insects gather to feed [7]. Mutual stimulation resulting from aggregation of solitary individuals in resources area leads to gregarization [12]. Thus, habitat fragmentation concentrates the locust population and increases the chances of contact (and gregarization), which favors recrudescence [12].

4.3. Analysis of the Error

As vegetation density has a strong influence on the measured surface reflectance, its impact on the detection performances was investigated by grouping the RAMSES samples into three density classes: “low” (25 points), “moderate” (625 points) “high” (368 points). The classification error decreases as the density increases for both detection DT and SVM (Figure 5). The decrease of ten percent between the “low” and the “moderate” density categories demonstrates the limit of the classifiers: a low vegetation density increases the risk of misclassification. The sensor’s spatial resolution (500 m) might also be a limiting factor for low vegetation density patch recognition.

The spatial distribution of the error was analyzed by stratifying the overall error by macroregion. This underlines the importance of the north-south climatic gradient on the detection accuracy (Figure 6). The error ranges from 18% to 34% for the DT and from 9% to 38% for the SVM classifier. However, these figures should be carefully interpreted with respect to the number of samples contributing to their computation (e.g., the error in the southwest region derived from only 11 points). The comparison of macroregions totaling a statistically sufficient number of samples leads to the following conclusion: the error rate increase in northern regions associated with more arid conditions. For instance, the southern region of Rguiba-Hodhein (216 samples) presents an error of 24% for the DT and of 20% for the SVM, as opposed to the northern region of Aftout (173 samples) where the error reaches 32% for both methods. It should be noted that even if located further north, the northwest region (202 points) displays a lower error rate than the Aftout region, especially for the SVM. This observation stresses the importance of the intrinsic characteristics of the ecological homogeneous regions within the climatic gradient. Indeed, regs of the northwest are covered by large vegetation patches when favorable rains occur, e.g., in 2010–2011, whereas the interdune landscape of Aftout favors discontinuous and limited vegetation patches.

Finally, the error repartition along time was investigated. Indeed, climatic variations along the seasons play an important part in plant development: the rainy season is associated with large and rather continuous and dense vegetation patches as opposed to the dry season, when vegetation is sparse. Ground data were clustered by month from August 2010–May 2011, and the overall error was assessed (Figure 7). In August, a typical month of the rainy season, the overall error is low for both classification methods (8% for the DT and 4% for the SVM). The error increases continuously during the following months with two peaks in October and January. High error values in April and May, the driest months of the covered period, stress the difficulty of accurate detection in the dry season. The observed error peak in October is explained looking at the spatial distribution of observation: ground observations shift to the north in October. In fact, survey teams follow the bio-geographical dynamics of locust populations that move up north to still vegetated areas. From an operational point of view, the interest of the dynamic dryness map lies at the end of the rainy season, at the onset of senescence. This onset appears in southern regions between September and October, but around January in northern regions. Hence, the onset of dryness in northern regions should be interpreted more cautiously than in southern areas.

Vegetation density appears thus as the driving factor for the error explanation, resulting in differential accuracies in both space and time according to the vegetation dynamics. Overall, the two classifiers are affected similarly by the density and the temporal and spatial factors. It is worth noting the interdependence of the factors: vegetation density depends on the climate and the season, as well as on the type of habitat. It is also worth mentioning that monitoring dry vegetation with the NDTI derived from MODIS uses spectral bands at 500-m spatial resolution. The monitoring is thus achieved at the cost of a loss in spatial details compared to the green area product (250-m). This would be overcome when considering higher spatial resolution satellites, such as Landsat-8 or the future Sentinel-2, which have coarser, but similar MODIS SWIR bands, allowing one to compute NDTI-like indices at 30 and 20 meters of spatial resolution, respectively. The potential of future monitoring based on 20/30-m rather than 250/500-m is clearly promising, particularly because the capability to identify vegetation density seems to play a critical role. This kind of spatial resolution would considerably reduce the resolution bias resulting from mixed coarse-resolution pixels [68] in fragmented areas, such as inter-dune vegetation patches or when vegetation patches cluster. This would also improve substantially the accuracy in the more arid northern areas. Nevertheless, the operational production of such 20/30-m dynamic maps represents a technical challenge, as the timeliness of their delivery is critical for planning the activities of national control centers and survey teams. Besides, providing an update on a 10-day basis would require the use of both Sentinel-2 A and B and most certainly their combination with Landsat-8 in order to ensure a sufficient temporal coverage. In the case of large cloud contamination, the fusion of high-resolution with low-resolution imagery ought to be considered.

Still, these findings have to be related to the bio-geographical cycle of desert locust in Mauritania. The senescence detection allows clearly a reliable mapping in summer breeding areas located in the south of the country. According to Babah and Sword [69], the south and southwest represent the main, if not the only breeding areas. Detections are less accurate in winter and autumn breeding areas dominated by fragmented vegetation patches where the vegetation density appears as the key limiting factor for detection at 500-m. These northern areas are mainly survival biotopes, but they remain fundamental for desert locust reproduction and, thus, for monitoring. Even if the assessment is limited to Mauritania, similar findings are expected in other areas of the distribution area, especially in those with similar climates.

4.4. Near-Real-Time Dynamic Dryness Mapping Case Study

In the previous section, it was found that the SVM and the DT perform equivalently. However, the DT looks more efficient for larger-scale implementation, because its simpler decision rules require less computing time, which is particularly valuable for near-real-time monitoring. In order to demonstrate the potential of the method to capture vegetation senescence dynamics, DT was applied every ten days on images covering a subset region of Mauritania from the beginning of July 2010 to the end June 2011, which simulates an operational context well. This area covers the intersection of the “southwest” and “Aftout” regions (Figure 8). The north-south climatic gradient is well represented, as shown by the small number of pixels detected by the green area in the northern part of the area labeled as bare soil in the land cover map (Figure 8a). The diagonal shown as bare soil corresponds to the highway linking Nouakchott to Nema. The presence of the town of Boutilimit and neighboring villages, with a high grazing pressure, explains the low canopy cover along the highway. The bare small circles enclosed in the vegetation (Figure 8b) correspond to villages. The landscape is mainly composed of narrow dunes with a NE-SW orientation and separated by rolling plains.

Consequently, the three classes (“growth”, “density reduction” and “drying”) were mapped every ten days (Figure 9). Grey areas are pixels labeled as grassland on the land cover map that are not detected as green vegetation during the considered 10-day interval. The vegetation growth appears gradually and first in the southeastern part and then spreads to northern regions, according to the rain progression along the rainy season. The vegetation density reduction first appears mainly in the south during the second 10-day interval of September 2010 (Figure 9f), then generalizes in October and, finally, leads to senescence in early November. However, this cycle does not occur homogeneously throughout the area. Drying patches appear during the first 10-day interval of September 2010 (Figure 9e), especially around the villages and the highway, while these pixels were not detected by the green area before and are classified as growing in late September (Figure 9g). Two areas identified as “growth” appear after a drying of the vegetation; one west near the highway on 11 October 2010 (Figure 9i) and one east during the next 10-day interval (Figure 9j). Some vegetation patches, predominantly in the north, become first senescent and loose their density afterwards as the transition from 01 October 2010 (Figure 9h) to 21 October 2010 shows (Figure 9j).

As these observations concern large areas, it limits the probability of misclassifications due to noise. Several reasons can explain the temporal evolution of these patches. First, sporadic rains heavily influences vegetation: a rainfall event, even limited spatially, may lead to a re-greening of the vegetation. Second, the vegetation density reduction can occur at any time by the combined action of wind and grazing pressure. This does not necessarily precede drying and can also follow it or occur concurrently. In the cases of density reduction occurring due to dryness, locusts should have already migrated, because they depend on green vegetation for feeding. Thus, in areas where the vegetation has already started to dry, the risk of gregarization by the concentration of populations is lower, and the survey teams may be redirected.

In terms of desert locust monitoring, the dynamic dryness map (as shown in Figure 9t) highlights the areas where vegetation has already dried out (grey), as well as where and for how long vegetation is in senescence (orange shades) or its density decreases (yellow-brown shades). This color code allows the interpreter to identify quickly, on the one hand, the areas from where the locusts are flying away (grey and yellow-brown shades) and, on the other hand, the areas still to be prospected, as density reduction might favor gregarization (orange shades). Those areas flagged in grey and yellow are less likely to shelter desert locusts, as they seek fresh and green vegetation to feed. Survey teams currently in those areas can be redirected to more favorable habitats (orange tones), translating into more efficiency monitoring in a cost-effective fashion.

5. Conclusions

Given the differential feeding behavior of desert locust and the extent and remoteness of their distribution area, remote sensing appears as a valuable tool for mapping their potential habitats. This study used 10-day MODIS mean composites over Mauritania to demonstrate the identification of the onset of dryness, as an indicator of a loss of attractiveness of habitats for locusts. The detection method relies on the temporal behavior of two indices: the Normalized Difference Vegetation Index, depending on green vegetation, and the Normalized Difference Tillage Index, sensitive to both green and dry vegetation. Two temporal metrics of NDVI and NDTI were selected by ROC curves for their ability to discriminate between three classes: “growth”, “density reduction” and “drying”. Three algorithms were benchmarked, and the decision tree was found to be both accurate (71%) and computationally economical. Smoothing the time series prior to the extraction of the metrics increased the overall accuracy by 5%. Several factors affected the classification accuracy: the north-south climatic gradient, the vegetation density, the period at which the vegetation was mapped and the relief. Finally, the method was applied over a southern region of Mauritania for the years 2010–2011. The resulting dynamic dryness maps appeared consistent with the seasonal vegetation cycle. The results obtained pave the way for the first operational implementation of the senescence dynamic maps and, consequently, further strengthen the capacity of the locust control centers, as well as supporting FAO’s preventative control strategy. As density appears as the key limiting factor, new and future satellites, such as Landsat-8 and the future Sentinel-2, will allow computing NDTI-like indices at 30 and 20 meters of spatial resolution, respectively, leading subsequently to a more accurate vegetation status recognition. The use of a longer time series of data is also expected to strengthen the results.