Desert Locust Cropland Damage Differentiated from Drought, with Multi-Source Remote Sensing in Ethiopia

Abstract

In 2020, Ethiopia had the worst desert locust outbreak in 25 years, leading to food insecurity. Locust research has typically focused on predicting the paths and breeding grounds based on ground surveys and remote sensing of outbreak factors. In this study, we hypothesized that it is possible to detect desert locust cropland damage through the analysis of fine-scale (5–10 m) resolution satellite remote sensing datasets. We performed our analysis on 121 swarm point locations on croplands derived from the Food and Agriculture Organization (FAO) of the United Nations, and 94 ‘non-affected’ random cropland sample points generated for this study that are distributed within 20–25 km from the ‘center’ of swarm affected sample locations. Integrated Drought Condition Indices (IDCIs) and Vegetation Health Indices (VHIs) calculated for the affected sample locations for 2000–2020 were strongly correlated (R² > 0.90) with that of the corresponding non-affected group of sample sites. Drought indices were strongly correlated with the evaluation Standardized Precipitation Evapotranspiration Indices (SPEIs), and showed that 2020 was the wettest year since 2000. In 2020, the NDVI and backscatter coefficient of cropland phenologies from the affected versus non-affected cropland sample sites showed a slightly wider, but significant gap in March (short growing season) and August-October (long growing season). Thus, slightly wider gaps in cropland phenologies between the affected and non-affected sites were likely induced from the locust damage, not drought, with fine scale data representing a larger gap.

1. Introduction

The Ethiopian economy is heavily dependent on smallholder rain-fed agriculture, which is vulnerable to natural and anthropogenic catastrophes that in turn leads to food insecurity [1]. The desert locust, Schistocerca gregaria (Förskal), is the most destructive herbivory migratory pest in the world [2]. It is highly mobile and feeds on large quantities of any kind of green vegetation, such as crops, pasture, and fodder. The desert locust has been a threat to agricultural production, causing starvation in Africa and western Asia for thousands of years [3,4]. The desert locust is normally a solitary insect that occurs in deserts and scrubs in northern Africa, the Sahel, the Arabian Peninsula, and parts of Asia to western India [3]. During this solitary phase, the desert locust has low occurrence and poses no threat to croplands. The 2019–2021 desert locust upsurge in East Africa and the Arabian Peninsula followed a series of heavy rainfall events beginning with Cyclone Mekunu in May 2018 that produced heavy rains in the Rub’al Khali of the Arabian Peninsula, ending in a drought and creating patches of vegetation [5], followed by Cyclone Luban in October 2018, leading to swarms being reported in early 2019. These conditions are conducive for the egg_laying, hatching, development and gregarious formation (swarms) of desert locusts [3,4,5]. In the gregarious phase, they can induce extensive agricultural production damage. Desert locusts may have a life cycle of two to three months, resulting in three to five generations in a year [3,4,5], and swarms can fly up to 100 km per day [3]. These factors of locust swarm development and mobility have made identifying their impact on smallholder croplands difficult with remote sensing.

Ethiopian vegetation has been disturbed by locust swarms since the end of 2019 (first swarm appearance on 25 June 2019) [6] and this round of outbreak has been the worst in the country for the past 25 years [7,8]. Locusts have invaded more than 200 Woredas (districts) which occupy over 200,000 hectares of land along the Rift Valley Escarpments, eastern and southern parts of the country [2]. The FAO has estimated the 2020 desert locust agricultural production damage in Ethiopia led to acute food insecurity for 8.6 million people between January and March 2021 [9]. The control effort of the outbreak was limited to merely minimal human power due to the complexities of war, bad politics, perilous terrain, and poor proactive institutional preparedness.

Desert locust research has typically focused on predicting their paths and breeding grounds based on ground surveys and satellite remote sensing of outbreak climate factors. Through a review of the literature, we observed a knowledge gap for detecting locust damage on croplands using satellite remote sensing datasets. Characterization of desert locust cropland damage could help efficient rehabilitation and mitigation efforts of food insecurity. Since higher rainfall and good vegetation growth favors locust infestation, we evaluated if fine spatial resolution satellite data (5–10 m), such as Planetscope, and Sentinel-1 and -2, would be able to detect damaged vegetation from locusts.

Drought is a common recurrent phenomenon in Ethiopia that has brought agriculture production shortfalls which leads to food insecurity and famine [10,11,12]. The most badly hit countries were Kenya, Ethiopia, Somalia, and Uganda [10,11] and it forced the worst humanitarian crisis to occur in the past 10 years. About 13.7 million people were at risk of hunger in Ethiopia [11]. Due to the 2015/2016 drought, about 18 million people in Ethiopia were in need of food aid (18 M ≈ 20%, i.e., 1 in 5 people needed food aid) [11,12,13,14]. In order to separate cropland damage from desert locust and drought impacts, we have also analyzed drought conditions in the region. Remote sensing-based integrated drought indices have the potential to describe drought conditions comprehensively. Analysis of multiple drought-related factors, such as precipitation, temperature, soil moisture, and vegetation health enables the effective characterization of the occurrence, spatial extent and intensity of drought in a given area [15,16,17].

2. Data, Study Region, and Methodology

2.1. Datasets and Study Region

The datasets used in this study, their descriptions, and sources, are presented in

Table 1. Descriptions of datasets used in this study.

S/N	Datasets and Sources	Spatial Res.	Temporal Res.	Study Duration	Data Variables
1	Locust Hub [6]	Point locations with impacted area (ha)	As it occurs	January 2020–February 2021	Locust sites
2	Planet Scope Basemaps [20]	5 m (resampled to 50 m)	Monthly	September 2020–February 2021	NDVI
3	Sentinel 2 L2A Surface Reflectance [29]	10 m (resampled to 50 m)	5 days *	January 2020–February 2021	NDVI
4	Sentinel-1 σ° [24]	10 m (resampled to 50 m)	6 days *	January 2020–February 2021	VV, VH
5	GFSAD30 [18]	30 m	Annual	January 2020–February 2021	Cropland points
6	MOD13Q1 [30]	250 m	16 days *	January 2000–December 2020	NDVI
7	CHIRPS † [31]	5 km	Daily *	January 2000–December 2020	rf
8	FLDAS †† [28]	10 km	Daily *	January 2000–December 2020	Ta, SM
9	GLEAM ††† [32]	25 km	Monthly	January 2000–December 2020	PET

^† Climate Hazards Group InfraRed Precipitation with Station data. ^†† Famine Early Warning Systems Network (FEWS NET) Land Data Assimilation System. ^††† Global Land Evaporation Amsterdam Model. * Aggregated to monthly datasets.

. We have used in situ and reanalysis data, as well as optical, active and passive microwave satellite datasets to distinguish locust and/or climate impacts on croplands.

Desert Locust swarm site datasets were downloaded from the Food and Agriculture Organization (FAO) of the United Nations website [6]. These data have observation date, geographic coordinates, areal extent affected by each recorded swarm, and other similar information. The FAO website has desert locust data in their various developmental stages and group formation, namely, bands, hoppers, adult, and swarms. We used the swarm stage dataset, since Desert Locust crop damage impact is the highest in their gregarious (swarm) stage. The graduated orange to red-colored sites in Figure 1, are 121 sample desert locust swarm sites used in this study, while the light green-colored sites correspond to cropland sample sites from the USGS Global Food Security-Support Analysis Data at 30 m (GFSAD30 [18]; totals 94; assumed non-affected by locust) used for comparative analysis. The graduated swarm sample sites are presented by month of locust occurrence from light orange (1 January 2020) to darker red (13 January 2021). Note that the graduated colors may not necessarily show swarm occurrence pattern in time. The first locust appearance in Ethiopia was the 25th of June 2019. The January 2020–January 2021 swarm locations presented here are selected samples, not all occurrences are presented here; hence, this sample may not show the distribution pattern through time. We only show when the swarms that were analyzed with remote sensing data, and most of the swarms in this sample occurred in October 2020.

We assume the temporal and spatial detail on 5–6 day intervals at 10–50 m resolution provided by Sentinel-1, and Sentinel-2, and monthly Planet Scope base maps at 5 m resolution will enable the detection of locust swarm damage.

Planet Scope Base maps of red, green, blue and near-infrared wavelengths at 5 m spatial resolution were provided free of charge by the Norway’s International Climate and Forests Initiative (NICFI) [20]. These data are available at monthly temporal resolution since September 2020, while it stretches back to December 2015 with a quarterly temporal resolution. The PlanetScope constellation of 3U CubeSats is comprised of Dove, Dove-R and SuperDove satellites in a 98° sun synchronous orbit with crossing times between 9:30–11:30 am [21,22].

The full Sentinel-2 mission comprises twin polar-orbiting satellites with a combined 5-day revisit time [23]. These satellites provide a global dataset since December 2018 and products are systematically generated at the ground segment. We used these level-2A product assets from the Google Earth Engine (GEE) Data Catalog to calculate NDVI. The assets contain 12 spectral bands representing surface reflectance (SR).

The Sentinel-1 mission provides a continuous all-weather, day-and-night imagery from a dual-polarization C-band Synthetic Aperture Radar (SAR) instrument at 5.405 GHz (C band) [24]. We have used data from the GEE collection that includes the S1 Ground Range Detected (GRD) scenes, processed using the Sentinel-1 Toolbox to generate a calibrated, ortho-corrected product. We processed the VV (single co-polarization, vertical transmit/vertical receive) and the HH (single co-polarization, horizontal transmit/horizontal receive) combination bands scenes.

The Moderate Resolution Imaging Spectroradiometer (MODIS) aboard the Terra and Aqua satellites was launched by NASA in 1999 and 2002, respectively, has 36 spectral bands from which different groups of data products are processed at different spatial and temporal resolutions tuned for specific applications. For this study, we used MODIS collection 6 level 3 Normalized Difference Vegetation Index (MODIS/006/MOD13Q1) data product from Google Earth Engine Data Catalog (Table 1). The MODIS NDVI product is computed from atmospherically corrected bi-directional surface reflectance that has been masked for water, clouds, heavy aerosols, and cloud shadows. Nadir bidirectional reflectance distributions functions (BRDF)-Adjusted Reflectance (NBAR) provides improved retrievals of surface reflectance through consistent normalization of multiple views of the surface to a nadir view using BRDF to model surface anisotropies [25,26].

The Climate Hazards Group InfraRed Precipitation with Stations (CHIRPS) dataset combines satellite precipitation estimates with in situ station rainfall data and has a spatial resolution of 0.05° [27]. An evaluation of CHIRPS, FLDAS, GPM/TRMM rainfall datasets with station rainfall dataset in the Ethiopian highlands showed that CHIRPS rainfall had the least bias and error (ME ≈ −0.2–0.2 mm, MAE ≈ 0.5–2 mm), and the best agreement (COR ≈ 0.8) [27].

Data from the Famine Early Warning Systems Network (FEWS-Net) Land Data Assimilation System (FLDAS) included gridded meteorological forcing fields and modeled hydrological variables [28]. The FLDAS air temperature field is derived from the National Oceanic and Atmospheric Administration (NOAA) Global Data Assimilation System (GDAS) and NASA Modern Era Reanalysis for Research and Applications version 2 (MERRA-2) datasets. The FLDAS soil moisture distribution data was derived from a combination of the Modern Era Retrospective-analysis for Research and Applications Version 2 (MERRA-2) and daily Climate Hazards Center InfraRed Precipitation with Stations (CHIRPS) data. We used the National Centers for Environmental Prediction/Oregon State University/Air Force/Hydrologic Research Lab (NOAH) model derived 0.1° × 0.1° daily Eastern Africa Region data product.

2.2. Methods

2.2.1. Desert Locust Swarm Data

For this study, generally, data variables were extracted for 2020 locust and cropland sample points using Google Earth Engine (GEE), and further processed using R. R is a programming language originally designed for statistical computing and graphics. Currently, R programming is widely used for processing remote sensing and GIS datasets. In addition, ArcGIS was used for buffer analysis around points, for zone preparation and sample selection. We generated corresponding sample points on neighboring croplands from the United States Geological Survey (USGS) Global Food Security-Support Analysis Data at 30 m (GFSAD30) using GEE. The locust datasets have swarm areal coverages between about 1 ha–35,500 ha. The following assumptions were made for this analysis: (1) locust swarm location points were taken at the geographic center of the respective swarms; and (2) cropland area is uniform within these study swarm points. Thus, (1) swarm-cover linear radius in meters:

$$\left(\mathrm{r}\right)=\sqrt{Area\left(ha\right)}\ast 50,\mathrm{i}.\mathrm{e}.,\mathrm{r}~50\mathrm{m}\text{–}10,000\mathrm{m}\left(10\mathrm{km}\right)$$

(1)

A double-ring buffer of 20 km (inner-ring) and 20–30 km (outer-ring) was created for all swarm points (Figure 1). Cropland sample points were clipped by the swarm outer-ring buffer (94 samples). To get a comparable number of swarm location sample points, we selected back swarm point samples that were within 25 km from selected cropland pixels. In such a process, we have selected 121 swarm sample points that are located between 20 km and 25 km from selected cropland sample points. Swarm points that were more than 25 km away from cropland sample points were discarded, which highly reduced the available swarm sample points. This was performed to be in line with the 25 km similar climatic influence zone assumption.

Assuming the swarm sites as center is the easier assumption that we propose, so that we are able to do a buffer analysis and select corresponding “non-affected” sites in reference to the swarm points. Fortunately, assuming the swarm sites as locust swarm edges will not affect our analysis, since we generate the corresponding “non-affected” sites 20 km away from swarm sites (between 20–25 km). Note that the reported larger locust swarm diameter is 20 km.

2.2.2. Approach for Analyzing Satellite Land Surface Remote Sensing Products

We applied Normalized Difference Vegetation Index (NDVI) [33] on both the aggregate swarm “affected” and aggregate “undisturbed” (non-affected) sample points for 2020 to estimate locust swarm damage to croplands at observational scales from 5 m to 250 m and spatially resampled these data to 50 m to 500 m to reduce atmospheric contamination. Cloud and cloud shadow masking, and speckle removal were applied to the optical and microwave products used in this study. The 250 m native spatial resolution 16-day MODIS surface reflectance data were resampled using a 2 × 2 square kernels mean value of the sample. We calculated Sentinel-2 L2A NDVI weekly data as per-pixel maximum value monthly composites. Our intent was to determine if each satellite time-series at moderate and coarse resolution had a significant departure from “normal” in swarm locations as compared to non-affected sites. We also calculated NDVI from Planet Scope Monthly Base maps for September 2020 through February 2021 and resampled to 50 m for both aggregated sites.

Sentinel-1 backscattering coefficient data from the Ground Range Detected (GRD) scenes Interferometric wide swath mode (IW) for the single co-polarization vertical transmit/vertical receive (VV) and dual-band cross-polarization vertical transmit/horizontal receive (VH) descending path was used in this study. Sentinel-1 may resolve differences in radar time series backscatter coefficients (σ°) from locust damage, as others have found utility in these data for assessing crop damage from other types of disturbance [34,35,36]. Both Sentinel-2 surface reflectance and Sentinel-1 SAR backscattering coefficients were resampled using 5 × 5 square kernels resulting in 50 m spatial resolution datasets.

2.2.3. Approach for Analyzing Climate Data

For drought analysis, we calculated soil moisture condition index (SMCI), vegetation condition index (VCI), temperature condition index (TCI), and precipitation condition index (PCI). Using these multi-indices will help us to detect the effects of all types of droughts (meteorological, hydrological, agricultural, and ecological). These indices were combined using optimal weights with respect to the corresponding Standardized Precipitation Evapotranspiration Index (SPEI) to get an Integrated Drought Condition Index (IDCI) and Vegetation Health Index (VHI) as in

Table 2. Descriptions of Remote Sensing Drought Indices.

Drought Indices	Data Source	Formula
PCI	CHIRPS	$\frac{({\mathrm{CHIRPS}}_{\mathrm{rf}}{}_{\mathrm{i}}-{\mathrm{CHIRPS}}_{\mathrm{rf}}{}_{\min })}{\left({\mathrm{CHIRPS}}_{\mathrm{rf}}{}_{\max }-{\mathrm{CHIRPS}}_{\mathrm{rf}}{}_{\min }\right)}$
SMCI	FLDAS	$\frac{\left({\mathrm{FLDAS}}_{\mathrm{SM}}{}_{\mathrm{i}}-{\mathrm{FLDAS}}_{\mathrm{SM}}{}_{\min }\right)}{\left({\mathrm{FLDAS}}_{\mathrm{SM}}{}_{\max }-{\mathrm{FLDAS}}_{\mathrm{SM}}{}_{\min }\right)}$
TCI	FLDAS	$\frac{({\mathrm{FLDAS}}_{\mathrm{ta}}{}_{\max }-{\mathrm{FLDAS}}_{\mathrm{ta}}{}_{\mathrm{i}})}{\left({\mathrm{FLDAS}}_{\mathrm{ta}}{}_{\max }-{\mathrm{FLDAS}}_{\mathrm{ta}}{}_{\min }\right)}$
VCI	MODIS	$\frac{({\mathrm{MODIS}}_{\mathrm{NDVI}}{}_{\mathrm{i}}-{\mathrm{MODIS}}_{\mathrm{NDVI}}{}_{\min })}{({\mathrm{MODIS}}_{\mathrm{NDVI}}{}_{\max }-{\mathrm{MODIS}}_{\mathrm{NDVI}}{}_{\min })}$
IDCI	CHIPS, FLDAS, MODIS	$\alpha \ast \mathrm{PCI}+\mathsf{\beta }\ast \mathrm{TCI}+\mathsf{\gamma }\ast \mathrm{SMCI}+(1-\alpha -\mathsf{\beta }-\mathsf{\gamma })\ast \mathrm{VCI}$
VHI	MODIS, FLDAS	$VHI=\alpha \ast \mathrm{TCI}+(1-\alpha )\ast \mathrm{VCI}$

α, β and γ represent the weight of single index while constituting the integrated drought indices; PCI: Precipitation Condition Index; SMCI: Soil Moisture Condition Index; TCI: Temperature Condition Index; VCI: Vegetation Condition Index; IDCI: Integrated Drought Condition Index; VHI: Vegetation Health Index.

. Such drought analysis methods have been well developed over several years [11,37,38,39]. Initial values for the weights were determined based on multiple linear regression of the factors with the SPEI time-series datasets. Then, other possible combinations closer to the initial combinations were made and a temporal analysis was conducted on a 5 km grid for all datasets. Note that the individual datasets used for the drought analysis have different native spatial resolution (Table 1). A 1-month, 3-month, 6-month, and 12-month seasonal accumulation of the CHIRPS precipitation data resulted in a similar temporal aggregation of the PCI and the corresponding derived SPEI. The drought condition indices were correlated with the corresponding SPEI, and thus we determined the best coefficients of the individual drought condition indices that yielded the more accurate IDCI and VCI. We calculated the SPEI from CHIRPS rainfall and the Global Land Evaporation Amsterdam Model (GLEAM) potential evapotranspiration (PET) that was mainly developed from in situ datasets [32,40,41]. In a previous study conducted in the same study area as our study, CHIRPS rainfall showed the least bias and error (ME ≈ −0.2–2 mm, MAE ≈ 1.5–2 mm), and the best agreement (COR ≈ 0.8), with station rainfall data [27]. Therefore, we are comfortable with using the CHIRPS rainfall data to drive SPEI, in combination with the PET data that in turn was derived mostly from in situ datasets, for evaluating the corresponding Integrated Drought Condition Index and Vegetation Health Index.

3. Results

In this research, we first assessed droughts in the region as it is prone to drought impacts. This was performed to determine whether vegetation damage in the region (if any) came from the impacts of locust, drought, or both.

3.1. Temporally Integrated Drought Condition and Vegetation Health Indices, and Evaluation

The highest correlation coefficient was obtained with a 6-month Integrated Drought Condition Index (IDCI6) calculated using a coefficient of 0.7 for the Precipitation Condition Index (PCI) and coefficients of 0.1 for the Condition Indices of Temperature (TCI), Soil Moisture (SMCI) and Vegetation (VCI). A high correlation value of r = 0.95 was obtained when the resulting 6-month IDCI (IDCI6) drought index was correlated with the corresponding 6-month Standardized Precipitation Evapotranspiration Index (SPEI6;

Table 3. Correlation coefficient (r) values between the drought condition indices with different weights and the different timescales SPEI for the non-affected and affected cropland sample sites (p-value < 0.01). Note the highest r with gray-shaded bold values. Seasonal PCI and IDCI were correlated with corresponding seasonal SPEI. Also note that the 1-, 3-, 6-, and 12-months PCI and IDCI were correlated with the corresponding 1-, 3-, 6-, and 12-months of SPEI respectively (similar font colors), while the only-monthly drought indices (TCI, SMCI, VCI, and VHI; yellow background) were correlated with all the 1-, 3-, 6-, and 12-months SPEI. Also note that the numbers at end of the indices refers to months.

Correlation Coefficient (r)
Drought Indices	Weight				Non-Affected Croplands				Affected Croplands
	PCI	TCI	SMCI	VCI	SPEI1	SPEI3	SPEI6	SPEI12	SPEI1	SPEI3	SPEI6	SPEI12
PCI1 PCI3 PCI6 PCI12					0.87	0.93	0.96	0.96	0.87	0.92	0.96	0.96
TCI1					0.38	0.45	0.48	0.43	0.39	0.45	0.47	0.41
SMCI1					0.73	0.75	0.61	0.39	0.74	0.72	0.56	0.36
VCI1					0.50	0.72	0.62	0.37	0.51	0.71	0.61	0.35
IDCI1 IDCI3 IDCI6 IDCI12	0.7	0.1	0.1	0.1	0.87	0.93	0.95	0.94	0.88	0.92	0.94	0.93
	0.6	0.2	0.1	0.1	0.86	0.92	0.93	0.91	0.86	0.91	0.92	0.90
	0.6	0.1	0.2	0.1	0.87	0.92	0.93	0.90	0.87	0.91	0.92	0.89
	0.5	0.1	0.3	0.1	0.86	0.90	0.89	0.84	0.86	0.89	0.88	0.83
	0.5	0.2	0.2	0.1	0.85	0.90	0.90	0.86	0.85	0.89	0.89	0.85
	0.5	0.1	0.2	0.2	0.85	0.91	0.90	0.86	0.86	0.90	0.89	0.85
	0.4	0.2	0.2	0.2	0.83	0.90	0.88	0.81	0.83	0.88	0.86	0.79
VHI1		0.2		0.8	0.54	0.76	0.68	0.45	0.56	0.74	0.66	0.42
		0.3		0.7	0.55	0.75	0.68	0.47	0.56	0.74	0.66	0.45
		0.4		0.6	0.54	0.73	0.68	0.49	0.56	0.71	0.66	0.46
		0.5		0.5	0.53	0.69	0.65	0.49	0.54	0.68	0.64	0.46
		0.6		0.4	0.50	0.64	0.62	0.49	0.52	0.64	0.61	0.46

, Figure 2-left). In contrast, the highest correlation coefficient was obtained with the 3-month Vegetation Health Index (VHI3) calculated using coefficients of 0.2 and 0.8 for the Condition Indices of Temperature (TCI) and Vegetation (VCI) respectively (Figure 2-right; Table 3). Scatterplots of these best-correlated seasonal indices (IDCI6 and VHI3) for the affected and non-affected sites is presented in Figure 3. The fitted trendline on this scatterplot yielded an R² of 0.98 and 0.94 for the IDCI6 (Figure 3-right) and VHI3 (Figure 3-left) respectively.

The two groups of study sample sites track a very similar IDCI and VHI seasonality with an r > 0.90 (Figure 2 and Figure 3). Time series of IDCIs, VHIs, and SPEIs revealed distinct drought and non-drought years in the study region. Drought years that showed growing season drought condition index values below 0.4 [37] include 2002, 2004, 2009, 2015, and 2017 (Figure 2 and Figure 4), while 2001, 2006, 2007, 2010, 2013, 2016, and 2018–2020 were non-drought years. The monthly and 3-months seasonal average SPEI displayed the shorter-time fluctuations of meteorological drought, while the 6- and 12-months seasonal average SPEIs displayed more aggregated drought conditions (Figure 4). Both the integrated drought indices (IDCI and VHI) and the evaluation drought index (SPEI) in Figure 2 and Figure 4 clearly revealed that 2020 was the wettest year since 2000 for both groups of sites.

3.2. Cropland Phenology Time Series

The Sentinel-2 site-average NDVI for 2020 for both the affected and non-affected sites tracked a similar seasonal pattern (Figure 5 top-center). Both NDVI time series showed a dual growing season (short Belg, and long Kiremt growing seasons). The affected sites NDVI time-series is lower (– 0.04 NDVI October 2020) than the corresponding non-affected sites NDVI in July, August, and September (Figure 5 center). The site-to-site standard error (SE) was slightly higher during the main growing season.

The MODIS NDVI also displayed the dual growing season (Figure 5 top-right). However, even though the affected sites NDVI was slightly lower than that of the non-affected sites during September and October, the difference was not significant (Figure 5 right).

Planetscope monthly mosaics did not display seasonality in affected sites, and only displayed a very slight seasonality in non-affected sites relative to MODIS and Sentinel-2 NDVI. Nonetheless, the gap is present between the affected and non-affected sites in the long growing (Kiremt) season and comparable to Sentniel-2.

The Sentinel-1 vertical-transmit and vertical-receive (VV; Figure 6 left) and vertical-transmit and horizontal-receive (VH; Figure 6 right) radar time series backscatter coefficients (σ°) also tracked similar cropland phenology for the two groups of sample sites. The backscatter coefficient timeseries showed a consistent distinct difference between the affected and non-affected sample sites throughout the year, in contrast to MODIS NDVI observations. The difference in σ° between the two groups of sites was higher and significant in April–May, July, and September–October (Figure 6). However, the non-affected sites σ° was consistently higher than that of the affected sites for the entire 2015–2021 data record (Supplementary Figure S1).

4. Discussion

4.1. Drought Conditions in Ethiopia

Throughout recorded history, East African agricultural production has suffered from recurrent drought conditions. The recent 2015/2016 and 2009 droughts are among many in the region [10,11,12,13,14]. In contrast, other years experienced floods. The integrated drought indices analyzed in this research, namely, the Integrated Drought Condition Index (IDCI) and Vegetation Health Index (VHI), and the evaluation drought index (Standardized Precipitation-Evapotranspiration Index–SPEI) clearly revealed that the year 2020 was the wettest year since 2000 for both the locust non-affected and affected groups of sites (Figure 2 and Figure 4). Therefore, the slight inferiority in NDVI and backscatter coefficient for the desert locust-affected sites in 2020 was not due to drought.

In a drought analysis drought indices timeseries, monthly and shorter season drought indices showed the shorter-time fluctuations of meteorological drought, while the longer-time indices showed agricultural drought. At the end of 2020 and beginning of 2021 in this study, some sites experienced meteorological drought conditions that might continue well into 2021 to bring agricultural drought.

4.2. Cropland Dynamics

A site-to-site NDVI standard error (SE) across Ethiopia was slightly higher during the main growing season. This difference might be attributed to the slight differences in growing season and crop types between the northern semi-humid and southeastern semi-arid parts of Ethiopia. Even though, for the locust-affected sites, MODIS NDVI was slightly lower than that of the non-affected sites during September and October in 2020 in this research, the difference was not significant. Thus, the finer composited spatial resolution Sentinel-2 NDVI (50 vs. 500 m) performed better than the corresponding MODIS NDVI to discern the locust damage on croplands in the study area. The MODIS NDVI in the study area displayed a higher phenology amplitude throughout the year in 2020 (above 3.5) compared to the corresponding Sentinel-2 NDVI. This can be attributed to the fact that the 2 × 2 kernel aggregated 250 m MODIS NDVI (yielding 500 m spatial resolution) may have natural vegetation contamination compared to that of the 5 × 5 kernel aggregated 10 m Sentinel-2 NDVI (yielding 50 m spatial resolution). In contrast, Sentinel-2 NDVI revealed a higher dynamic range potentially due to improved sensitivity of the finer spatial resolution to vegetation dynamics throughout the growing season.

The limited duration of Planet Scope monthly Basemaps product (September 2020 onwards) reduced our ability to assess the NDVI gap throughout the year 2020 (unlike the other two NDVI products analyzed in this research–MODIS and Sentiniel-2), and ability to detect both growing seasons. However, the gap is present in the long growing (Kiremt) season and comparable to Sentinel-2. These results capture the current potential for Planet products to be used for locust swarm damage assessment. Future products from Planet with refined atmospheric correction models applied to derive surface reflectance could enhance the crop monitoring aspect of these data, and potentially enhance the locust signal we are trying to resolve.

The backscatter coefficient (σ°) for affected sites was consistently lower both during 2015–2019 when no locusts were present and during 2020 when locusts were present, suggesting the differences seen were not mainly due to locust damage but because the affected sites were intrinsically different, indicating locust habitat preferences. However, these gaps in 2020 were wider during April–May, July, and September–October, which is consistent with the other NDVI indices above that occurred during the desert locust active duration. While the precise characteristics of these preferences are uncertain, the lower coefficients are consistent with a lower vegetation biomass/structure and lower aggregated soil moisture. Lower vegetation biomass may aid such behaviors as basking to increase body temperatures while well-drained soils in more open areas can be preferred sites for laying [42,43].

The vegetation damage gaps we detected in this study during the small and main growing season in Ethiopia are in line with other studies. In a World Bank-Supported High-Frequency Phone Surveys report, there were two waves of locust invasion in Ethiopia [2]. The first invasion was from January to May, while the second round started in late September and peaked in October to November. The first round of locust invaded ~180–240 Woredas (Districts) and damaged over 200,000 hectares of mainly croplands, but also grazing lands and trees in eastern and southern Ethiopia. This round of locust vegetation damage corresponds to the March small growing season (Belg) phenology gap between the locust-affected and non-affected groups of sample sites in this study. The second round of invasion was reported to be 20 times more severe than the first round, and this was exacerbated by COVID-19 restrictions, war, and flooding that hindered locust management [2]. According to the World Bank report, due to the locust crop damage and related factors in Ethiopia, poverty levels were expected to rise and make over 2 million people to fall into poverty. The scale of Productive Safety Net Programme (PSNP) was expected to increase from 9 million to as high as 15 million [2]. Established in 2005, the PSNP is a multi-donor trust fund Ethiopia’s Government multi-billion dollar food security, public works, and social safety net program for millions in need across Ethiopia [44,45,46]. This study also found wider NDVI time-series gap in the main growing season (particularly October–November) between the locust-affected and non-affected sample sites.

5. Conclusions

NDVI and backscatter coefficient phenology from locust-affected sites (121) and adjacent non-affected sites (94) in Ethiopia for the year 2020 showed closely similar climate patterns. The gaps between affected and non-affected sites are slightly wider during March (shorter growing season) and August–October (late long growing season) with Sentinel-2 10 m resolution NDVI data. We conclude this can be attributed to the active locust swarms during these times in our study domain. Our meteorological and agricultural drought indices in the study area showed that the 2019/2020 growing season was in a non-drought phase. NDVI difference from Planet Scope Basemaps at 5 m resolution in the major locust-affected and important crop producing parts of the country also revealed that the 2020 NDVI in the late long growing season (September–October) was inferior compared to that of the 2021, while corresponding climatic conditions were comparable.

Through our analysis we found that the largest challenge in large area locust damage assessment with remote sensing includes: (1) locust infestations occur following strong rainfall that can support good crop vegetation growth; (2) locusts move toward good crop growing areas for feed; and (3) damage is variable, crops can regenerate after invasion, and the sensitivity of damage from 10 m (Sentinel-2) to 250 m (MODIS) is not consistent. Locust outbreaks continue to be difficult to detect with remote sensing, yet it appears Sentinel-2 had the ability to detect early and late season NDVI differences from the 121 sites we investigated. These combined factors make it very difficult to discriminate locust-affected and non-affected croplands with moderate and fine scale multi-temporal satellite remote sensing. We found phenology gaps between the two groups of samples are wider for finer resolution Planet, Sentinel-1 and Sentinel-2 datasets compared to that of the moderate resolution MODIS dataset, but the signal of locust outbreak remains quite small. Analysis of optical and/or microwave meter to sub-meter level very-high spatial resolution datasets collected at dense temporal intervals acquired from satellites (e.g., from Planet Scope [3–7 m], future WorldView Legion constellations, etc.), aerial (manned and unmanned aerial vehicles), and/or ground-based technologies may provide a more robust signal to detect cropland locust damage. Analyzing locust band cropland damage might also be an important consideration for future study, as locust hoppers are less mobile than adults.