Dynamics, Anomalies and Boundaries of the Forest-Savanna Transition: A Novel Remote Sensing-Based Multi-Angles Methodology Using Google Earth Engine

Abstract

This chapter proposes a remote sensing multi-angles methodology to assess the transition at the interface of the forest-savanna land cover. On Sentinel2-A median images of successive dry seasons, three referential and nine analytical spectral indices were computed. The change vector analysis (CVA) was performed, selecting further one magnitude per index. The averaged moving standard deviation index (aMSDI) was proposed to compare spatial intensity of anomalies among selected CVA, and then statistically assessed through spatial and no-spatial autoregression tests. The cross-correlation and simple linear combination (SCL) computations spotted the overall anomaly extent. Three machine learning algorithms, i.e., classification and regression trees (CART), random forest (RF), and support vector machine (SVM), helped mapping the distribution of each specie. As result, the CVA confirmed each index ability to add new information. The aMSDI gave the harmonized interval [0–0.083] among CVA, confirmed with all p−values=0, z−scores>2.5, clustering of anomaly pixel,and adjusted R2≤0.19. Three trends of vegetation distribution were distinguished with 88.7% overall accuracy and 0.86 kappa coefficient. Finally, extremely affected areas were spotted in upper latitudes towards Sahel and desert.

1. Introduction

The global structure and productivity of ecosystems are deeply impacted by joined climate conditions and human drivers, causing general vegetation degradation [1]. The phenomenon has kept increasing in the last four decades and eventually affects whole ecosystem, soil productivity, biological systems, biotic diversity, and other environmental systems’ ability to support human needs in concerned areas. Main indicators are the decline in parameters, such as low biomass, less ecological production, fragmentation, or lower canopy cover [2, 3].

Inside the tropics, vegetation is globally sensitive to seasonal and inter-annual variation in precipitations and temperatures. Extremes seasonally, i.e., longer rainy season and shorter dry season in lowest latitudes, versus the reverse phenomenon towards medium latitudes, influence the vegetation distribution with several phenological and physiological adaptations, including cover and status changes [4, 5, 6]. Typically, forest colonizes wetter areas, while savannas cover drier areas, with a gradual species distribution such as dense forest, tree savannas, grassy/herbaceous savannas, and isolated desert shrubs or clumps of dry grasses known as steppes. However, transitions are not rigidly determined by climate [7]. There is an extensive overlap between forest and savanna creating a mosaic of landscapes, and most studies on the subject remain widely hypothesized and modeled with controversial results, supported by questionable evidences. Biases include the high species turnover around 1000 mm to 2500 mm rainfall, the (un)stable states of forest and savanna maintained by feedbacks between tree cover and disturbances, and for the satellite-based approaches, the structural (in)difference between trees or grasses layer [8].

These specificities are challenging to spatialize at a point that sub-Saharan African ecosystems have played a key role in the development of remote sensing of vegetation for decades [9, 10, 11]. Nowadays, several satellite-based models provide scalable spectral information relevant to vegetation distribution and changes, physiology, and phenology, in broad terms, to monitor and combat land degradation, especially in African drylands [12, 13, 14, 15]. As such, numerous spectral indices measure the vegetation parameters [16]. The Normalized Difference Vegetation Index, NDVI, especially, has purposely been widely used [17]. However, some limitations like sensitivity to soil background effects and atmospheric influence as well as values saturation under dense and multilayered canopy, usually alter the NDVI capacity to simultaneously predict senesced vegetation and efficiently discriminate individual anomalies, i.e., growth, vigor, leaf area index, biochemical components (anthocyanin, carotenoids, cellulose, etc.), water content or pigmentation [18, 19, 20] with accuracy. Then some previous studies focused on identifying or modeling the direction of change as well as underlying drivers of drylands vegetation [21]. Those models applied to two or more spatially close and interwoven vegetation species, require to implementation of specific processings. To the best of our knowledge, the recent progress in modeling sub-Saharan vegetation transition introduced the term of “bistability” around lower and upper transition boundaries between forest and savanna [8]. This model is based on paleo-ecological evidences (soil, topography) and climatic parameters change and oscillation (rainfall and temperature), that influence (for the firsts) and predispose (for the seconds) these two species to coexistence. With the support of a floristic survey, the ambiguity of mischaracterizing savanna as a degraded forest was clarified at some point, by identifying, the dense forest, the “bistable” forest, the “bistable” savanna, and the proper savanna.

This study approaches the forest-savanna transition study, by investigating its different spectral behaviors and their statistical meaning, in a context lacking field data, and from an open-source/open-data perspective. The triple aim is to assess the dynamics, discriminate species disturbance based on their empirical spatial distribution, and predict their extent and boundaries. As such, assuming a blurred boundary and an overlapping spatial gradient between the two species, some phenological and physiological characteristics are considered as separately as possible in terms of anomalies, and further integrated beneath the same model, so as to locate the spots requiring permanent monitoring or sustainable actions, without mischaracterizing punctual changes, factual distribution, and most accurate delineation.

2. Materials and methods

2.1 Study area

The study was conducted on the sub-Saharan mixed ecoregion, highly dependent on varied annual precipitations (AP) and yearly medium to high average temperatures (T⁰), which both influence relative humidity (RH) changes. The area belongs to the medium Cameroon (central Africa), between latitudes 5⁰0’-8⁰5’N and longitudes 10⁰0′-15⁰5′E, a climatic transition between the agroecological zones of western highlands (AP = 1800-2500 mm; T⁰ = 19.5°C; RH = 75%), bi-modal rainforest (AP = 1700-2000 mm; T⁰ = 24°C; RH = 80%), and then Guinean high savanna (AP = 1500–1800 mm; T⁰ = 30°C; RH = 60%) to Sudano-Sahelian savanna (AP = 400–1200 mm; T⁰ = 28°C; RH = 50%) for the core area. Specifically, the vegetation density broadly reflects the climate gradient of dense moist broadleaf forest and highland forest to sparse extensive savannas featuring the co-dominance of woody shrubs, grassland in plains, and herbaceous steppe at the edge of Sahel (Figure 1) [22].

Figure 1.

Location of the study: (A) Cameroon in the context of African ecoregions based on Olson et al. (2001); (B) distribution of ecoregions in Cameroon; (C) preview of the subset ecoregions; (D) Sentinel2-a median image of the subset for mid-November 2020 to end march 2021.

2.2 Data and working environment description

The study was conducted using European Spatial Agency, ESA, Sentinel2-A multispectral instrument (MSI) data, which represents a very valuable opportunity for the fine characterization and monitoring of vegetation types on large scales, but is poorly investigated for the tropical biome study [11, 23]. This sensor provides 13 varied spectral bands from 0.443 to 2.190 micrometers, a 10-day repeat cycle, and a spatial resolution up to 10 meters (Additional material 1).

The phenological dry season, globally from November to March, was selected because of its high temperatures and less rainfalls, assuming they are ideal conditions to observe the vegetation adaptations to extremes, for years 2015 to 2021. In the GEE cloud coding environment and using the JavaScript opensource simplified coding, the median reducer function was appended as the pixel-wise computation of all bi-annual collection images, based on a band per band processing [24]. Then, applying a date filter from November 15 to end of March 31, a boundary filter for the study area, and a cloud cover acceptance filter below 10%, one image of 13 bands was outputted per bi-annual periods 2015–2016 to 2020–2021, displayed a,nd converted from 16 bits to 8 bits before further processings. Offline tools, i.e., desktop software, were used to extract statistics and for some complementary processings using less memory, including final layouts. Namely, Erdas Imagine 2020, ArcGIS Pro 2020, and Microsoft Excel with extension Xlstat were specifically used.

2.3 Inter-seasonal dynamics, spectral assumptions, and casting of indices

In arid and semi-arid regions, common changes in density, spatial distribution, chlorophyll, pigmentation, water stress, anthocyanin, nitrogen, carotenoids, leaf structure, and browning or senescence differently impact the biomass [25]. Previous spectral indices-based applications investigated that, the visible (0.4–0.7 μm) wavelengths respond to photosynthetic and non-photosynthetic pigments, the NIR (0.7–1.4 μm) wavelengths respond to the cellular structure and exhibit solar-induced fluorescence (SIF) and the SWIR (1.4–3 μm) wavelengths respond to senescent non-photosynthetic vegetation. As such, the anisotropic behavior of vegetation at visible-SWIR wavelengths has been parameterized to describe vegetation structure [26].

We computed twelve spectral indices, whose three were selected as reference data according to their ability to better highlight the land cover targeted, i.e., second modified soil adjusted vegetation index, MSAVI2 [27], to assess the vegetation cover, normalized difference soil drought index NDSoDI [14], to highlight the dry bare soils, and new water index, NWI [28], to map the surface water. After testing dozens of other indices, two assumptions were emitted for an efficient last casting, such as i) two or more indices targeting the same anomaly, with different or improved (increase/decrease) spatial patterns, express new information to be consider; ii) two or more indices targeting a different anomaly should show different spatial patterns to be considered enough informative. Accordingly, nine vegetation indices were selected and computed, i.e., photosynthetic vigor index (PVR), global vegetation moisture index (GVMI), plant biochemical index (PBI), first red-edge inflection point (REIP1), modified chlorophyll reflectance index green (mCRIG), leaf water content index (LWCI), modified anthocyanin reflectance index (MARI), simple ratio pigment index (SRPI), and plant senescence reflectance index (PSRI). Table 1 gives details of their formulations with references.

Name	Equation	Primary goal	Reference
MSAVI2	$\frac{2\mathit{NIR}+1-\sqrt{{\left(2\mathit{NIR}+1\right)}^2-8\left(\mathit{NIR}-\mathit{Red}\right)}}{2}$	Map vegetation cover while reducing soil background	[27]
NDSoDI	$\frac{\left(\mathit{Red}+\mathit{\text{SWIR}}1\right)-\mathit{\text{Blue}}}{\left(\mathit{Red}+\mathit{\text{SWIR}}1\right)+\left(\mathit{\text{Blue}}+L\right)}\left(1+L\right),$; $\mathit{\text{with}}L:\left[0.4\right]$	Map dry soils while reducing surroundings source of moisture	[14]
NWI	$\frac{\mathit{\text{Blue}}-\left(\mathit{NIR}+\mathit{\text{SWIR}}1+\mathit{\text{SWIR}}2\right)}{\mathit{\text{Blue}}+\left(\mathit{NIR}+\mathit{\text{SWIR}}1+\mathit{\text{SWIR}}2\right)}$	Extract surface water	[28]
PVR	$\frac{\mathit{\text{Green}}-\mathit{Red}}{\mathit{\text{Green}}+\mathit{Red}}$	Photosynthetic vigor for crop monitoring	[29]
GVMI	$\frac{\left(\mathit{NIR}+0.1\right)-\left(\mathit{\text{SWIR}}1+0.02\right)}{\left(\mathit{NIR}+0.1\right)+\left(\mathit{\text{SWIR}}1+0.02\right)}$	Vegetation water content and evapotranspiration	[30]
PBI	$\frac{\mathit{NIR}}{\mathit{\text{Green}}}$	General biochemical reflectance	[31]
REIP1*	$700+40\left(\frac{\left(\frac{\mathit{Red}+\mathit{\text{Rededge}}3}{2}\right)-\mathit{\text{Rededge}}1}{\mathit{\text{Rededge}}2-\mathit{\text{Rededge}}1}\right)$	Canopy chlorophyll, nitrogen content and polluted soil dynamics	[32]
mCRIG	$\left({\mathit{\text{Blue}}}^{-1}-{\mathit{\text{Green}}}^{-1}\right)\mathit{NIR}$	Chlorophyll, carotenoids and anthocyanin	[33]
LWCI	$\mathit{\text{LWCI}}=\frac{log1-\left(\mathit{NIR}-\mathit{\text{SWIR}}1\right)}{-log1-\mathit{NIR}-\mathit{\text{SWIR}}1}$	Water content, change and stress	[34]
MARI	$\left({\mathit{\text{Green}}}^{-1}-{\mathit{Red}}^{-1}\right)\mathit{NIR}$	Chlorophyll, carotenoids and anthocyanin	[35]
SRPI**	$\frac{\mathit{\text{DeepBlue}}}{\mathit{\text{Rededge}}2}$	Nitrogen content, water and chlorophyll	[36]
PSRI	$\frac{\mathit{Red}-\mathit{\text{Blue}}}{\mathit{NIR}}$	Pigmentation and vigor changes dues to vegetation senescence	[37]

Table 1.

Spectral indices used.

*

Was inverted after first preview, to better highlight vegetation patterns instead of soil dynamics.

**

Was adapted to SENTINEL2-A covered wavelengths by using band 6.

These indices were stacked and previewed with, MSAVI2, NDSoDI and NWI seeking the following: i) agreement or disagreement with MSAVI2; ii) absolute disagreement with NDSoDI; iii) disagreement when possible NWI (Additional material 2). On the transversal transect of 1544 pixels covering all types of vegetation and land cover, 500 pixels were sampled on different sections (100/300 pixels) for each index, to compare trends and relationships with MSAVI2 and NDSoDI. The two periods moving average was used to better visualize trends, while the simple linear regression was performed to extract the coefficient of correlation, ${\mathit{R}}^{\mathbf{2}}$, and the root mean square error, RMSE (Figure 2). Based on patterns distribution and these statistics, a threshold was defined for each index, to separate vegetated and non-vegetated areas (Additional material 3). The binarized images were used as entries for the change assessment.

Figure 2.

Relationships between each analytical spectral index and both referential ones. The bright green curve represents the dependent variable, while the two others are explanatory variables (MSAVI2 = dark green; NDSoDI = light brown). N = 500 for all variables and regression parameters. The following is noticed from the trends of curves: PVR, PBI and GVMI have a neat positive correlation with MSAVI2, but a sharp negative one with NDSoDI; conversely, MARI and PSRI curves describe the exact opposite trends (positive with NDSoDI and negative with MSAVI2); mCRIG and LWCI curves evolve in another direction cutting the MSAVI2 and NDSoDI curves, while IREIP1 curve shows no real relationship them.

2.4 Change vector analysis and spatial dynamics assessment

This process creates a difference image between two or more bands in a multi-temporal image analysis, so to detect changes in the type or the conditions of surface features. Depending on the study, change vector analysis (CVA) can use calculation principles of the Mahalanobis or the Euclidean distance as in this study, according to the following expression [38]:

$$\mathit{R}=\sqrt{{\left({\mathit{\beta \sigma }}_{\mathbf{1}}-{\mathit{\beta \sigma }}_{\mathbf{2}}\right)}^{\mathbf{2}}+{\left({\mathit{\beta \rho }}_{\mathbf{1}}-{\mathit{\beta \rho }}_{\mathbf{2}}\right)}^{\mathbf{2}}}$$E1

Where $\mathit{R}$ is magnitude of the vector change, ${\mathit{\beta \sigma }}_{\mathbf{1}}$ and ${\mathit{\beta \sigma }}_{\mathbf{2}}$ are fraction in date 1 and date 2 respectively, ${\mathit{\beta \rho }}_{\mathbf{1}}$ and ${\mathit{\beta \sigma }}_{\mathbf{2}}$ are fraction cover in date 1 date 2.

Four classes of magnitude are represented for either degradation or re-growing [39] (Additional material 4). For each binarized image, a CVA process was performed, the magnitudes of increase and decrease patterns were used, while the stability magnitudes were ignored. A total of twelve difference images per index resulted, i.e., six per selected magnitude. Considering the need of complementarity among indices, only one CVA magnitude, i.e., increase or decrease, was selected per index for further processing. Criteria used to optimize the selection were the original goal of each index, as well as its spatial and statistical relationship with MSAVI2, NDSoDI and NWI (Table 2).

Index	Targeted asset	Compared trends to MSAVI2 / NDSoDI / NWI	CVA
PVR	Age & Cover	Same with a plus / Inverted / Inverted	Decrease
GVMI	Water content & Cover	Same with a plus / Inverted / Inverted	Decrease
PBI	Inner composition & Cover	Same with a plus / Inverted / Inverted	Decrease
IREIP1	Height, Cover & Health	Different / Inverted / Inverted	Decrease
mCRIG	Health & Composition	Different / Inverted / Inverted	Increase
LWCI	Water stress	Different / Inverted / Close	Decrease
MARI	Health & Inner composition	Inverted / Close / Inverted	Increase
SRPI	Health & Inner composition	Inverted / Close / Inverted	Increase
PSRI	Age & Cover	Inverted / Close / Inverted	Increase

Table 2.

CVA magnitude selected per spectral index.

2.5 Averaged moving standard deviation index method and anomalies mapping

The moving standard deviation index, MSDI, is a filter applied to satellite images multispectral or derivative channels using the moving standard deviation calculation, generally to assess degradation [40]. One common application is the vegetation and soil of semi-arid systems, where the variability of the MSDI is used to indicate levels of habitat degradation [41, 42, 43]. MSTDI has been proven efficient to operate well in complex regions [40, 42, 44].

Here, the five derivative images per selected CVA were used as entries. A standard deviation was computed for each entry. Difference between consecutive standard deviations were calculated, giving four new images. Then, the averaged MSDI (aMSDI) was conceived as follows:

$$\mathit{\text{aMSDI}}=\frac{{\sum }_{\mathit{i}=\mathbf{1}}^{\mathit{n}}{\mathit{\sigma }}_{\mathit{i}}{\mathit{b}}_{\mathit{i}}}{\mathit{n}+\mathit{\tau }}$$E2

For all $\mathbf{1}\mathit{\le }\mathit{i}\mathit{\le }\mathit{n}$, ${\mathit{\sigma }}_{\mathit{i}}{\mathit{b}}_{\mathit{i}}$ represent the ${\mathit{i}}^{\mathit{th}}$ CVA-band contribution to the information distribution in the outputted image, whereas $\mathit{\tau }$ is the difference between the number of times each index was originally calculated (7 times, starting from the multispectral median images) and the number of resulting standard deviations computed among selected CVA-bands (5). Here, $\mathit{\tau }\mathit{=}\mathbf{2}$. Computations were performed at $\mathbf{3}\mathbf{\times }\mathbf{3}$-, $\mathbf{5}\mathbf{\times }\mathbf{5}$- and $\mathbf{7}\mathbf{\times }\mathbf{7}$-pixel moving window. The same steps were applied to both of decrease and increase CVA magnitudes of MSAVI2. The goal at this point was to assess spatial convergence and divergence trends with the nine analytical vegetation spectral indices’ aMSDI, as a first-hand validation process of anomalies distribution. Two autoregression-based tools were used as metrics at this step (ESRI, 2014):

• The global Moran’s I index was computed to assess the spatial autocorrelation of aMSDI outputs along the 1544 pixels transect. Its integration in spatial analysis is important to avoid incorrect statistical inference from inefficient or biased parameter estimates [45]. The formula is expressed such as follows:

$$\mathit{I}=\frac{\mathit{n}}{{\mathit{S}}_{\mathit{O}}}\frac{{\sum }_{\mathit{i}=\mathbf{1}}^{\mathit{n}}{\sum }_{\mathit{j}=\mathbf{1}}^{\mathit{n}}{\mathit{w}}_{\mathit{i},\mathit{j}}{\mathit{z}}_{\mathit{i}}{\mathit{z}}_{\mathit{j}}}{{\sum }_{\mathit{i}=\mathbf{1}}^{\mathit{n}}{\mathit{z}}_{\mathit{i}}^{\mathbf{2}}}\left[-\mathbf{1},\mathbf{1}\right]$$E3

Where ${\mathit{z}}_{\mathit{i}}$ is the deviation of an attribute for pixel $\mathit{i}$ from its mean ${\mathit{x}}_{\mathit{i}}\mathit{-}\overline{\mathit{X}}$, ${\mathit{w}}_{\mathit{i}\mathit{,}\mathit{j}}$ is the spatial weight between pixel $\mathit{i}$ and $\mathit{j}$, $\mathit{n}$ is equal to the total number of pixels, and ${\mathit{S}}_{\mathit{O}}$ is the aggregate of all the spatial weights, developed as:

$$S_O=\sum _{i=1}^n\sum _{j=1}^n{w}_{i,j}$$E4

Values closer to $\mathbf{0}$ indicate little or no spatial autocorrelation, regular values indicate a clustered trend, and negative values indicate a dispersed trend. In this paper and based on the proposed model, i.e., one CVA magnitude per index and average of the MSDI among years, the clustered trends were targeted to express values closeness. The null hypothesis for the patterns analysis was emitted for a less randomness and a significant clustering, comparing how good would be the model to spot and regroup spatially related pixels of anomalies, and discriminate them from “healthier” pixels. High standard deviation (all $\mathit{z}\mathit{-}\mathit{\text{scores}}\mathit{>}\mathbf{2.5}$) and low probability (all $\mathit{p}\mathit{-}\mathit{\text{values}}\mathit{=}\mathbf{0}$) were ranges targeted, meaning that the observed spatial patterns are probably too unusual to be the result of random chance.

• Ordinary Least Squares (OLS) regression was computed as the non-spatial measure of the spatial convergence/divergence of aMSDI. This method predicts or models a dependent variable relationship with a set of explanatory variables. The regression coefficients are usually estimated by using least-square techniques such as:

$$y=\alpha +\sum _{j=1}^{n}x{\beta }_j+\epsilon$$E5

Where, $\mathit{\alpha }$ is the intercept, ${\mathit{\beta }}_{\mathit{j}}$ is the regression coefficient, and $\mathit{\epsilon }$ is the residual term, given by the difference between observed and expected value. Under these assumptions, the regression coefficient can be obtained by [46]:

$$\widehat{\mathit{\beta }}={\left({\mathit{X}}^{\mathit{T}}\mathit{X}\right)}^{-\mathbf{1}}{\mathit{X}}^{\mathit{T}}\mathit{y}$$E6

Here, each aMSDI of the nine indices was considered as dependent variable, whereas MSAVI2’s aMSDI of both CVA, were simultaneous explanators. Three OLS parameters were selected for interpretation, i.e., adjusted R-squared, corrected Akaike’s information criterion (AICc) which assesses the best-fit model between spatial and non-spatial Ordinary Least Squares (OLS) models, root mean squares error (RMSE), and three graphs, i.e., scatterplots of the relationship among variables, histograms of standard deviation probability and plots of residuals versus predicted.

At this point, the spatial correlation was assessed through the cross-correlation mapping process using only the finest $\mathbf{3}\mathbf{\times }\mathbf{3}$ -pixel moving window, that consistently shows best statistic performances than the two others [43]. The parameters of the algorithm were set for the maximum gap (maximum pixel shifting) at $\mathbf{1}$, and the maximum masked fraction (maximum fraction of the pixel within the correlation window that are allowed to be masked) at $\mathbf{0.5}$.

Further, the individual maps of anomalies were used to produce a single one. The first step consisted in stacking the aMSDI for each $\mathbf{3}\mathbf{\times }\mathbf{3}$ -pixel moving window. Then, another correlation map was computed for all the nine inputs, only selecting the correlation coefficient layer. The highest values of the outputs were assumed show agreement and disagreement among aMSDI for high and low anomalies. At last, to keep the highest values as signs for anomaly spots and the lowest ones as potentially “healthier” areas, a simple linear combination, SLC, was performed among the nine aMSDI for each moving window size, by summing the layers without a weighing factor.

2.6 Machine learning for vegetation discrimination

Based on visual trends, binarized images were regrouped in three axes of three indices of vegetation each. A principal component analysis, PCA, was performed for each ax using the covariance reducer algorithm, that reduces some number of $\mathbf{1}\mathit{-}\mathit{D}$ arrays of the same length $\mathit{N}$ to a covariance matrix of shape $\mathit{N}\mathit{\times }\mathit{N}$. This reducer uses the one-pass covariance formula [24]. The first component of each output was chosen, and a threshold was set to separate vegetated and non-vegetated areas. Then the three outputs were stacked with NDSoDI and NWI to form a new five bands image. In Ref. to both this new image and to the last study period median image (2021–2022), three classes of vegetation were defined, completed by soil/Built-up and water classes. At least two spectral curves were produced per class to confirm the trends, except for water that appeared uniform. The land use land cover, LULC, classes sampling was performed using 115 points, and a process of $\left[\mathbf{1}\mathit{:}\mathbf{1}\right]$ ratio between training samples on the PCA stacked image and testing samples on the MS image (Table 3).

Training:Testing					Total
Vegetation1	Vegetation2	Vegetation3	Soil/Built-up	Water
100:100	100:100	100:100	100:100	100:100	500:500

Table 3.

Samples distribution.

Three machine learning algorithms were performed, such as the Classification and Regression Trees (CART), Random Forest (RF) and Support vector machine, (SVM) [47, 48, 49]. The metrics used to assess general performance of each classifier were overall accuracy, OA, and the kappa coefficient, KC [50, 51]. Whereas, the thorough assessment of their efficiency on individual LULC classes was done using the error matrix to extract the producer accuracy (PA) and the user accuracy (UA). Eqs. (7) and (8) formulate the main quality measures of the learning:

$$\mathrm{OA}=\frac{\mathit{\text{Total number of correct samples}}}{\mathit{\text{Total number of samples}}\left(\%\right)}$$E7

$$\mathrm{KC}=\frac{{\mathit{\epsilon }}_1-{\mathit{\epsilon }}_{\mathbf{2}}}{1-{\mathit{\epsilon }}_{\mathbf{2}}}$$E8

with,

$${\mathit{\epsilon }}_{\mathbf{1}}=\frac{{\sum }_{\mathit{i}=\mathbf{1}}^{\mathit{n}}{\mathit{D}}_{\mathit{ii}}}{\mathit{N}},\text{and},{\mathit{\epsilon }}_{\mathbf{2}}=\frac{{\sum }_{\mathit{i}=\mathbf{1}}^{\mathit{n}}{\mathit{D}}_{\mathit{i}+}{\mathit{D}}_{+\mathit{i}}}{{\mathit{N}}^{\mathbf{2}}}$$

Where, ${\mathit{D}}_{\mathit{ii}}$ is the number of observations in row $\mathit{i}$ and column $\mathit{i}$ of the confusion matrix, $\mathit{n}$ is the number of rows in the error matrix, $\mathit{N}$ is total number of counts in the confusion matrix, ${\mathit{x}}_{\mathit{i}\mathit{+}}$ is the marginal total of row $\mathit{i}$, and ${\mathit{x}}_{\mathit{+}\mathit{i}}$ is the marginal total of column $\mathit{i}$. Main steps of the processing are summarized in Figure 3.

Figure 3.

Overview of the study design.

3. Results

3.1 Magnitudes of spatial dynamics

Globally, each index showed a different interseason sensitivity to dynamics for the two main magnitudes of change in the vegetation distribution. The indices detecting vigor, biochemical composition, and water content, i.e., PVR, PBI and GVMI, that are highly correlated with MSAVI2, invaded the upper medium area for the increasing magnitude, while their increasing patterns were more concentrated on the south and lower medium area (Figure 4a and b). Besides, the indices detecting canopy chlorophyll, carotenoids, anthocyanin, nitrogen, and water stress, with a positive (mCRIG) or non-significant (IREIP1, LWCI) correlation to MSAVI2, gave a general decreasing or increasing patterns on the whole subset (Figure 4a and b). Whereas, the indices targeting the lack of chlorophyll, carotenoids, anthocyanin, pigment or the plant senescence, and that were negatively correlated to MSAVI2, i.e., MARI, SRPI and PSRI, cover the South area for the decreasing patterns, and the North part for the increasing patterns (Figure 4a and b).

Figure 4.

a. Decrease trends of CVA per index. The background layers are from 2021 to 2022. b. Increase trends of CVA per index. The background layers are from 2021 to 2022.

Moreover, all magnitudes agreed to the same total areas for all indices, when adding the unchanged magnitudes. For the extreme cases, LWCI indicates a brutal decrease in between 2015 and 2016 and 2016–2017 (Change1, Figure 5a), from 61,609 km² to 22,030 km², i.e., 58% of the total change for the study period (Figure 5b). While, GVMI indicates an important increase between 2019 and 2020 and 2020–2021 (Change 5, Figure 5a), from 1439 km² to 41,215 km², i.e., 69.4% of the total change for the study period (Figure 5b). Then all changes were complementary between the two magnitudes, except for IREIP1 in the last period, 2020–2021/2021–2022, which the decrease was unsignificant (4km²) compared to the increase (Figure 5a&b). Above all, the magnitudes of change were compared in terms of hierarchy tree, to measure indices that better assess the spatial progression or regression. For the decrease magnitudes, LWCI assessed regression for 4/5 intervals, except during the period 2018–2019 to 2019–2020 (Change 4, Figure 5a), which recorded a recover followed by another regression. While for the increase magnitudes, mCRIG gave the most accurate assessment of a spatial progression for 5/5 intervals, which is the overall best score magnitudes. These statements justify their dominant position in the tree-map (Figure 5c).

Figure 5.

Areas(a), percentages(b) and hierarchy of changes (c) based on the CVA decrease (row1) and increase (row2) magnitudes.

3.2 Trends and intensity of anomalies

3.2.1 Contribution of aMSDI and spatial autoregression significance

The outputs of the aMSDI were all obtained in the same interval, [0–0.083] for each CVA-band, with a maximum value confirmed above ${\mathbf{10}}^{\mathit{-}\mathbf{6}}$, in every pixel window size. It should be noticed that the binarized inputs do not interrupt the spatial gradient of the phenomena described and the values showed a continuum of intensities. Therefore, although the results are unitless, the lowest value, [0], predict no potential issues, and the highest value, [0.083], predict a critical intensity of the targeted degradation or anomaly. Globaly, the spatial distribution of seasonal degradation/anomaly evolves in the south-north direction (Additional material 5 a-d).

From the spatial autocorrelation model synthesized in Table 4, the expected Moran’s I index and p-values are identical for all, −0.000648 and 0, respectively, while variance is almost the same for all, 0.000457 or 0.000458. Besides, final Moran’s I index varies between 0.81 for MSAVI2 aMSDI of increase CVA, and 0.95 for PVR; whereas z-scores are all above 2.5, i.e., between 37.9 for MSAVI2 aMSDI of increase CVA, and 44.3 for PVR. Therefore, the spatial patterns of every highlighted vegetation’s distribution issue by each aMSDI, described an important clustering of pixels with same and closest values at a confidence level of 99%, stating that the phenomenon is far from randomness as wished (ESRI, 2014) (Additional material 6). Added to this, the proposed model of aMSDI has been significantly stable around all parameters, while its values scale is easy to compare and interpret among outputs of different analytical spectral indices.

aMSDI	Moran’s I	Expected Moran’s I	Variance	Z-Scores > +2.58	P-Values
MSAVI2_Decr	0.9	−0.000648	0.000457	42.2	0.00
MSAVI2_Incr	0.81	−0.000648	0.000457	37.9	0.00
PVR	0.95	−0.000648	0.000458	44.3	0.00
PBI	0.85	−0.000648	0.000458	39.6	0.00
GVMI	0.84	−0.000648	0.000457	39.5	0.00
IREIP1	0.76	−0.000648	0.000457	35.4	0.00
mCRIG	0.86	−0.000648	0.000457	40.3	0.00
LWCI	0.91	−0.000648	0.000457	42.5	0.00
MARI	0.91	−0.000648	0.000457	42.5	0.00
SRPI	0.86	−0.000648	0.000457	40.4	0.00
PSRI	0.84	−0.000648	0.000458	39.3	0.00

Table 4.

Spatial autocorrelation model report. N = 1544 for all variables and parameters.

3.2.2 Spatial convergences of anomalies and non-spatial autoregression analysis

First outputs of the cross-correlation mapping algorithm informed on the spatial relationship between each MSAVI2’s aMSDI and individual aMSDI of analytical indices. About the MSAVI2 decrease CVA, the highest positive correlation of aMSDI was with PVR (R = 0.4; R² = 0.16; RMSE = 0.026), while PBI (R = -0.23; R² = 0.05; RMSE = 0.027) and MARI (R = -0.22; R² = 0.05; RMSE = 0.027) were found totally decorrelated (Figure 6). As such as, the most noticeable anomaly impacting the forest-savannah vegetation cover is the decrease of vigor and other subsequent characteristics highlighted by PVR and close algorithms. Less issues might be caused by biochemical composition including anthocyanin. Concerning the MSAVI2 increase CVA, the highest positive correlation of aMSDI was established with LWCI (R = 0.53; R² = 0.28; RMSE = 0.016), at the opposite of decorrelation to MARI (R = -0.56; R² = 0.32; RMSE = 0.015) and IREIP1 (R = -0.55; R² = 0.31; RMSE = 0.015) (Figure 6). With its particular distribution on the East area, LWCI express that water content influence to the increase of vegetation, while anthocyanin and canopy chlorophyll have no influence.

Figure 6.

Cross-correlation maps between aMSDIs of MSAVI2 and the selected aMSDI for each index.

The multi-regression performed between each analytical aMSDI and simultaneously both MSAVI2’s aMSDI on $\mathbf{3}\mathit{\times }\mathbf{3}$-pixel moving window, were consistent but specifics by case. Generally, the relationships were all highly significant (all $\mathit{p}\mathit{=}\mathbf{0}$). The corrected Akaike Information Coefficient (AICc) vary between −7795 and − 6642, for IREIP1 and PVR respectively; the intercept values are in between [0.005–0.03]; the regression coefficient of each MSAVI2 is more significant with indices potentially spotting high and low anomalies at the same areas, but both negatives for LWCI and MARI (Figure 7). The adjusted correlation coefficients were low (all ${\mathit{R}}^{\mathbf{2}}\mathit{\le }\mathbf{0.19}$) and those of aMSDI sharing identical patterns with MSAVI2 (PVR, PBI and GVMI) recorded a normal probability distribution of standardized residuals around $\mathbf{0}$, whereas other are skewed to the right. In both cases, there is a dominant peak, marking the separation among affected and non-affected areas, or between covered and uncovered areas. All the graphs of residuals show less randomness and more clustering per class of standard deviation (Additional material 6).

Figure 7.

OLS regression results synthesis (Eqs. (5) and (6)). All estimates have, p = 0 and n = 1544. Values in parentheses indicate standard error estimates for regression parameters. MSAVI_Decr/MSAVI_Incr = values of aMSDI for decreasing or increasing patterns of the second modified soil-adjusted vegetation index; $R^2=$ adjusted correlation coefficient; RMSE = root mean square error; ∆AICc = indicates the corrected Akaike information criterion, thus the spatial model improvement over the non-spatial form.

Spatial synthesis of the cross-correlation mapping algorithm confirms the distribution of patterns and global trends. The combined correlation among the nine indices’ aMSDI were found strong at different degrees, for both highest and lowest values, then indicating convergences of significant or non-significant anomalies at the South-west and North-East of the area (Figure 8). For the three windows of calculation, they were clearly separated from predicted vegetation statuses, i.e., highly and lowly affected, soil, built-up and water features (Figure 8). From the cross-correlation map synthesis, the simple linear combination (SLC) outputs for each window discriminated positive from negative spatial trends. Then, the largest spots of extremely severe anomalies are located in the south-west, south-east and north-east areas, well separated on $\mathbf{5}\mathit{\times }\mathbf{5}$ and $\mathbf{7}\mathit{\times }\mathbf{7}$-pixel moving windows. Subsequently at the finest $\mathbf{3}\mathit{\times }\mathbf{3}$-pixel moving window, the spots of anomalies are continuously concentrated in the south-west to north-east direction (Figure 8).

Figure 8.

Anomaly trends and spots. The synthesized correlation maps on the first row shows patchwork of high correlation in the medium area for the $\mathbf{3}\mathit{\times }\mathbf{3}$-pixel moving window, but two large continuous spots of convergence in the north-east and south-west according to $\mathbf{5}\mathit{\times }\mathbf{5}$ - and $\mathbf{7}\mathit{\times }\mathbf{7}$-pixel moving windows. The SLC outputs on the second row definitely separate high values from low values correlations.

3.3 Spatial extent and predicted delineation

The repeated sampling of vegetation gave spectral curves with at least three pairs of same trends for vegetation (Figure 9a). As spatial evidence, the patterns of ML were identical for the three algorithms, with three distinguished classes of vegetation, well separated from soil/built-up and water (Figure 9b). The three classifiers performed with a high and identical OA of 88.7%, for a KC of 0.86 (Figure 9b). These measures confirm the spectral agreement between classified and multispectral image reflectance of the five land features. When confronting PA and UA, that also recorded identical values for all classifiers per class (Figure 9b), it can be assumed that regrouping indices thresholded images is a good option to identify the forest-savanna transition and interpenetrated vegetation. Presumably and based on the world ecoregions map as well as the empirical vegetation density, the three classes of vegetation can be identified as follow in south-north latitudes direction (Olson et al., 2001): i) first ax matches tropical and sub-tropical moist broadleaf forests; ii) second ax is dominated by grasslands, mixed with savannas; iii) third ax mixes savannas and shrublands.

Figure 9.

Spectral curves (A), spatial patterns (B) and statistics (C) of the three axes of vegetation on LULC. In details: (A) CART; (B) RF; (C) SVM: (D) stacked axes of vegetation with soil (NDSoDI) and water (NWI) features for 2021–2022; E) Sentinel2-a image of 2021–2022.

Concerning the vegetation types, discrimination, and extent, RF and SVM gave the exact areas for the three classes, CART agreed with them for the second and third axes, whereas the stacked derivative only agreed with all the algorithms concerning the third ax, i.e., 3807.3 km² (Figure 9c). The discrepancies especially with the stacked image could be explained by the soil features footprint that is different between the NDSoDI displayed and classified image, as well as the urge overlapping noticed earlier between second and third axes. In fact, the interpenetration of the predicted grassland/savanna (axe2) and savanna/shrubland (axe3), which might better reflect the transiting vegetation behavior targeted along the study, raises on the other hand assumptions on the accurate assessment of areas. Their UA (grassland/Savanna = 73.3% and savanna/shrubland = 94.4%) then support this assumption (Figure 9b).

4. Discussions

The efficiency of the whole proposed methodology was assessed and discussed on selected aspects and the comparison with existing methods was basically empirical. At first, depending on their goals, previous MSDI-based studies analyzed only the standard deviation of the red and near infrared wavelengths, while those integrating vegetation indices were limited to three of them [40, 43, 52]. Because the goal of proposing aMSDI in this study was to assess consecutives dry season anomalies and discriminate them from empirical statuses of the forest-savanna specificities, we integrated nine spectral indices, selected on the basis of targeted phenological or physiological weaknesses, and whose computations basically integrate several wavelengths. Interestingly, although only one CVA magnitude was chosen per index, all individual models showed the expected visual convergence of similarity or dissimilarity trends.

Moreover, previous applications stated that the common calculation of MSDI on raw spectral index, gives outputs with a minimum value of zero and a maximum value determined by those of the pixels evaluated [43]. Consequently, outputs value cannot be directly compared. Here, by applying the averaging process to binarized CVA, [0,1], this study alternatively addressed these issues for a multidate analysis. With same or divergent visual patterns, all outputs were scaled inside the identical interval, [0–0.083], although showing convergent or divergent spatial patterns. The significance of the spatial and non-spatial autoregression models, has helped to confirm the inner variability of each aMSDI although identical scales of values and apparent same trends between others, as well as convergences/divergences of trends with the reference aMSDI (MSAVI2). Consequently, the output spotting anomalies was proof of complementary among individual contributions and spatial agreements.

Besides, common attempts of mapping distribution, typology, and delineation of forest and savanna, have always been supported by fieldwork, based on climate parameters, as well as including paleo-ecological evidences and detailed floristic survey to be efficient [8]. The methodology presented in this paper has predicted three different axes of vegetation, resulting from the PCA processing and thresholding (Table 5). For each of the six study periods, the first ax in the south part is composed by a dense and potentially healthier vegetation, highly correlated with the referential data MSAVI2. Whereas the other two axes, more and more sparse towards north, are divergent with the first one and somehow each with another (Figure 10a&b). Nevertheless, their individual trends foresee some overlapping, in the center and in the northern areas. A sampling of each ax along the transect of 1544 pixels, showed how interweaved and complexes are the boundaries among forest-savanna species (Figure 10b).

Vegetation trends	Regrouping	PCA threshold
1^st axe	PVR, PBI, GVMI	$\ge 0.65$
2^nd axe	IREIP1, MARI, PRSI	$\ge 0.95$
3^rd axe	mCRIG, LWCI, SRPI	$\ge 0.75$

Table 5.

Vegetation axes, proposed groups and PCA thresholds used to binarize.

Figure 10.

a. the three main axes of vegetation’s spatial distribution. b. Pixels’ value along the transect of the vegetation’s axes.

To answer the interrogations behind these ambiguities, a simple multicollinearity test was run, showing how independent one ax is from another. When the correlation between two independent variables is considerably high, it is a problem in the modeling process. The VIF (variance inflation factor) and tolerance were used for diagnosis. VIF is the reciprocal of tolerance, knowing that, tolerance is $\mathbf{1}\mathit{-}{\mathit{R}}^{\mathbf{2}}$. We used the lowest known VIF, $\mathit{<}\mathbf{3}$, while expecting the highest tolerance, so to measure independence. Therefore, while the occurrence percentiles of value 1 on each ax highlights interweaving in between [62.5–78.3] for axes 1 and 2, and [69.5–99.97] for all, low VIF ($\mathbf{1.33}\le \mathit{VIF}\le \mathbf{1.35}$) and high tolerance ($\mathbf{0.74}\le \mathit{Tol}\le \mathbf{0.89}$) confirm the total separability, i.e., less collinearity among axes (Table 6).

Axe	VIF	Tolerance	Percentile for value 1
AXE1	1.33	0.75	[49.9–99.97]
AXE2	1.12	0.89	[62.5–99.97]
AXE3	1.35	0.74	[78.3–99.97]

Table 6.

Multicollinearity test results.

Although from this study, we cannot properly use the qualifier of “bistable” forest or savanna, because it highly depends upon climate and paleo-ecological parameters, it is important to notice how ambiguous is the distribution and blurry are the boundaries. Thoroughly, on any ML output, three zooms distributed on three different latitudes helped to notice different types of transitions (Figure 11). Between the lower latitudes 5⁰30’-6⁰30’North, the transition is from the first (moist broadleaf forest) to third ax of vegetation (shrubland savanna), although the second ax (grassland savanna) would have been “expected.” Between the middle latitudes 6⁰30’-7⁰30’North, the transition mixes in the below area, the “-‘unexpected” third ax (shrubland savanna) with the ‘expected’ first ax (moist broadleaf forest) of vegetation, before the wide expansion of the ‘expected’ second ax (grassland savanna). At this point, the only “expected” transition was inside the upper latitudes 7⁰30’-8⁰30’North, where the second ax (grassland savanna) gradually gave way to third ax (shrubland savanna) of vegetation. These elements of analysis support the qualifier of “bistable” area, while still questioning the anisotropic distribution with latitudes, and encouraging the finest scale of analysis, i.e., spatial and spectral resolution.

Figure 11.

Zooms on the transitions, a sign of anisotropic distribution with latitudes. Yellow square = lower latitudes (5⁰30’-6⁰30’N) transition; red square = middle latitudes (6⁰30’-7⁰30’N) transition; blue = upper latitudes (7⁰30’-8⁰30’N) transition.

Finally, the display of anomalies with the LULC classes disambiguated the confusion of savanna and degraded forest. The observation was made by overlaying the highest and the lowest values of anomalies in the most concentrate area, on the SVM output. On three spots covered by grassland, shrubland and bare soil, the modeled extremely severe anomalies concern just a part of each class. Whereas, on two spots of lower to no-anomalies, savanna as well as bare soils are partially concerned (Figure 12).

Figure 12.

Anomalies versus LULC classes. (A) SVM classification map. (B) Two extreme classes of anomalies at $3\times 3$- pixel moving window size. (C) & (D) subset of comparison among classes of LULC and anomalies. (E) Stacking results with the following details: Black, yellow, purple, and ginger pink circles = extremely severe anomaly spots extended to part of the entire vegetation, in the grassland/savanna –savanna/shrubland transition area; orange circle = low or inexistent anomaly spot in the savanna/shrubland transition, characterized as drier and more exposed vegetation to degradation; cyan circle = low or inexistent anomaly on bare soil.

5. Conclusions

This study has conducted an experimentation on the forest-savanna vegetation, with the goal of assessing dynamics, assuming anomalies and predicting boundaries. On Google Earth Engine platform and using Sentinel2-A satellite images of seven consecutive dry seasons, from 2015 to 2016 to 2021–2022, twelve spectral indices were selected according to their different phenological and physiological assessment of the vegetation, and other natural features to be discriminate. Using the processing of change vector analysis, CVA, it was successfully observed that each index brings a substantial information, to better assess increase or decrease patterns of the vegetation cover. Further, proposing the averaged moving standard deviation index, aMSDI, to face potential issues of simple MSDI, the scale of spatial trends appraisal was found identical between the same interval [0–0.083] for all pixel window sizes, while keeping spatial trends as specific as they are for each selected CVA. As confirmation, all $\mathit{p}\mathit{-}\mathit{\text{values}}\mathit{=}\mathbf{0}$, $\mathit{z}\mathit{-}\mathit{\text{scores}}\mathit{>}\mathbf{2.5}$ there is a high clustering between anomaly pixels, whereas low adjusted R² among each analytical index aMSDI and MSAVI2 ones validate the performance of the model. Besides, three main trends of vegetation emerged, i.e., moist broadleaf forest in the south, grassland mixed to savanna in the core and savanna mixed to shrubland in the north, based on CART, RF and SVM classifiers performed on thresholded, PCA regrouped and stacked bands, with 88.7% OA and 0.86 KC. Finally, taking all the nine aMSDI as entries, a paired cross-correlation mapping helped to identify same general trends for high and low values. Whereas, the application of simple linear combination, SLC, highlighted the important spots of anomalies distribution in the northern part of the subset towards Sahel and desert, but less concentrated in the southern part towards moist forest area. Because the forest-savanna anisotropy with latitudes remains questionable depending on the scale of the study, it can be inferred from the study that, although aMSDI method shows capabilities in semi-dry areas, the choice of indices to use is the responsibility of the author. Whereas, a strong correlation remains to be investigated between the seasonal anomalies and the causes or drivers.

Acknowledgments

To USGS for the free availability of satellite data. To all experts-developers of Google Earth Engine platform APIs. To StatsN’Maps Consulting Firm, for the logistic support. To our laboratories.

Conflict of interest

The authors declare no conflict of interest.

Links to the code are available upon request.

Additional material 1. Characteristics of Sentinel-2A bands.

	Sentinel-2A MSI
Name	Range (μm)	Bands (Resolution/m)
Coastal aerosol	0.421 − 0.457	B1(60)
Blue	0.439 − 0.535	B2(10)
Green	0.537 − 0.582	B3(10)
Red	0.646 − 0.685	B4(10)
Red_edge1	0.694 − 0.714	B5(10)
Red_edge2	0.731 − 0.749	B6(10)
Red_edge3	0.768 − 0.796	B7(10)
NIR wide	0.767 − 0.908	B8(10)
NIR narrow	0.848 − 0.881	B8A(20)
Water vapor	0.931 − 0.958	B9(60)
Cirrus	1.338 − 1.414	B10(60)
SWIR1	1.539 − 1.681	B11(20)
SWIR2	2.072 − 2.312	B12(20)

Additional material 2. The spectral indices used. The first row is the supporting or reference data.

Additional material 3. Binarized indices and thresholds defined.

Additional material 4. Direction and magnitude of change as proposed by Kuzera et al. (2005).

Additional material 5a. Averaged MSDI patterns of MSAVI2.

Additional material 5b. Averaged MSDI patterns for PVR, PBI and GVMI.

Additional material 5c. Averaged MSDI patterns for IREIP1, mCRIG and LWCI.

Additional material 5d. Averaged MSDI patterns for MARI, SRPI and PSRI.

Additional material 6. Spatial autoregression sample plots of aMSDI for PVR (left) and MSAVI2 increase patterns (right). red square = spatial correlation targeted, for n = 1544 pixels.

Authors contributions

Conception, A.H.N.M.; study design, A.H.N.M.; acquisition of data, A.H.N.M., I.C.N.P., F.C.L.T., L.M.B., M.T. and J.V.M.M.; execution, A.H.N.M.; analysis and interpretation, A.H.N.M.; writing-original draft preparation, A.H.N.M.; writing-review and editing A.H.N.M., I.C.N.P., F.C.L.T., L.M.B., M.T. and J.V.M.M.; All authors have read and agreed to the published version of the manuscript.