Minding Spatial Allocation Entropy: Sentinel-2 Dense Time Series Spectral Features Outperform Vegetation Indices to Map Desert Plant Assemblages

1. Introduction

Dryland vegetation communities endure remarkable ecological resilience and adaptation strategies that are increasingly challenged by uncertain climate patterns and changing ecological conditions. Globally, drylands comprising arid and hyper-arid regions has been estimated at 27% of the terrestrial surface area [1]. Local scale mapping and conservation of, often heritage and ecologically unique, arid plant communities is gaining global attention considering the upsurge in developmental actions that undermine the spatial and temporal shifts in plant habitats due to changing climates. The differentiation of environmental conditions at local scale is an important determinant of plant biodiversity and distribution [2,3]. Such differentiation, mediated by variability in topography, promotes the formation of micro-refugia for plant species [4] – stressing the importance of microtopographic variability. Variability in the spatial configuration of terrain contributes to diversifying potential plant habitats, on one hand. On another, by decreasing the spatial range of individual habitats and their likely connectivity, it mitigates inter-species competition (for limited space and resources); thus, contributing to plant micro-refugia and species’ richness [4,5,6]. The mapping of these refugia and vegetation assemblages is imperative in monitoring spatial changes and in planning the protection and conservation of their contributed and heritage biodiversity.

Land use and land cover maps, by and large, are increasing serving as support tools for both research and end-user applications as predicting impacts of climate crisis and planning development or management strategies [7,8,9,10]. Since maps are a simplified representation of ground reality, errors and uncertainties are inevitable in their production. These errors vary in magnitude depending on one of several, often compounding, factors including data quality, data processing, and analysis procedure [11]. Thus, an important consideration of such maps is their reliability or uncertainty – often evaluated by metrics of map quality or accuracy. For effective remote sensing application, different approaches have been suggested for accounting and reporting uncertainties in classification at the pixel level [12], object-based [13,14,15] and scale-based or spatial uncertainty [16,17]. In general, different analysis and target application of land use and land cover, whether pixel- or object-based, warrant corresponding scale of accuracy assessment [18]. Thus, approaches combining multi-features (several indices, texture) at high spatiotemporal resolution are promising in differentiating classes of mixed vegetation or land cover [19,20,21,22], which may inform several intended conservation and management endeavours in drylands and desert landscapes.

In remote sensing of vegetation, the use of Vegetation Indices (VIs) have been largely successful in mapping high and dense vegetation landscapes [19,22]. Such vegetation retain leaf and related phenological cycles over prolong periods, which feeds sensitivity of VIs and other proxies as feature texture. Leaf phenology and texture is captured in aerial or satellite images and across time series to inform the differentiation and mapping of different and mixed vegetation [15,21,22,23,24]. However, such proxies of vegetation may be applicable and reliable in arid and desert landscapes to a lesser extent. In arid and desert landscapes, permanent vegetation prevails in cultivated farms and active oases, in existent montane forests, sparse woodlands. A considerable proportion of the vegetation has ephemeral distribution, with very brief periods of green leaf phenology – chlorophyll retention – during periods of rainfall and water flushes into the Wadis (water channels and catchments). Since these periods of vegetation proliferation are often erratic and short-live, field observations of vegetation occurrence and distribution may not be readily reflected in remote sensing-based proxies of vegetation status. And, these proxies may be less sensitive, or the lack thereof, to observed non-green, dried, and sparse vegetation [25,26,27]. Thus, it is unrealistic to continuously apply classic proxies and assumptions, mostly reliable to permanent vegetation distribution and density, in mapping vegetation in arid and desert landscapes. Usage of spectral-based VIs, as spatial predictors, to map and monitor arid vegetation is increasingly questionable [11,26,28] and may be unreliable for different and specific landscapes [11,20,25,27,29,30,31]. Thus, for brief and mostly dry vegetation in arid landscapes, it is unreasonable to continuously apply or assume unverified reliance on Vegetation Index (VI) proxy from multi-spectral remote sensing data. Such proxy may be reliable under certain landscape specifics – such as in areas with persistent vegetation cover - and in other circumstances provide only complementary information in discriminating vegetation assemblages and mapping their distribution.

The spatial and temporal heterogeneity of dryland vegetation are key contributors to errors in mapping their distribution. Nonetheless, literature on remote sensing of such vegetation has overlooked quantifying the uncertainties in vegetation mapping from spatial features of spectral signatures and vegetation indices. Different methods have been suggested in assessing accuracy or errors of remote sensing classification of vegetation and land cover at the pixel-based level [32,33], local spatial scale [34], object-based feature extraction [35], and map change classification [36,37]. In addition to classification error and confusion, spatial consideration of unknown space, for which no reference predictor information is available, is essential in estimating area of applicability of map outputs [16,17]. The application of remote sensing in pixel-based (soft) classification produces errors that vary across space. Thus, usage of overall accuracy estimates, and related metrics as kappa, provide less intuitive visual clues to, and interpretation of, the spatial variability and patterns in classification errors [32,38].

Uncertainties in mapping vegetation have been largely overlooked in remote sensing of vegetation in arid, semi-arid and desert landscapes [39,40,41,42]. There is scant literature that either validate or benchmarks the reliability and spatial extent of VIs-based mapping of desert plant communities for applications as biodiversity conservation and development planning. This study, therefore, questions and dissociates, based on quantified classification uncertainties and spatial disagreement, the contribution of VIs proxies versus principal spectral information in discriminating and mapping vegetation assemblages in desert landscapes. Using a dense time-series of Sentinel-2 imagery and the random forest (ensemble features) classification algorithm, we quantify classification error from both per-pixel class uncertainty by Shannon entropy (information loss) and feature-based uncertainty by intensity analysis of class transition across classification models.

2. Materials and Methods

2.1. Study Area

The study landscape is Wadi Ashar in AlUla County (in Al Madinah District) of the Kingdom of Saudi Arabia. The landscape’s physical structure and geology comprise a mix of several spectacular natural features and landforms (sandstone massifs, harrat plateaux, steep gorges, sand dunes, sandy plains, open Wadis, or seasonal drainage channels) and land uses (irrigated farms, date palm orchards, animal camps, etc.). The elevation across the landscape ranges between 760m a.s.l in the Wadis to 1240m a.s.l on the Harrat plateaus, while the terrains have slope range of 0 to 75 degrees. The administrative boundary covers a surface area of about 3,046 ha (Figure 2), and the landscape has established and planned development strategies and facilities of both public, touristic, and regal interests - exquisite hotel resorts. At its centre is the iconic Maraya Concert Hall as the centre structure of a developed chain of exquisite desert resorts which are punctuated by heritage landscape features and touristic attractions.

2.2. Vegetation Classification and Mapping

In detecting and mapping the different field-surveyed vegetation communities or assemblages, we implemented supervised classifications based on the random forest ensemble (bagging and decision tree) algorithm [43]. As predictors for the classifier algorithm, we explore spatial information of spectral bands and their derivative vegetation indices from time-series Sentinel-2 multi-spectral imagery. The details of the reference ground truth data, spatial data, and classification procedure are provided in the following subsections.

2.2.1. Plant Survey (Field) Data

Field data and observations are initially a mosaic of different key vegetation types in the study area. The vegetation types were preliminarily categorised into a high order based on three major habitats: sandy plains, sand dunes and harrats. To better understand these habitats and generate vegetation descriptions and maps, extensive field data were sampled during field surveys - the surveys were conducted from November to December 2021 and in September 2022. During the surveys, two types of plant data were collected: 1. Plot data - a complete list of the plants in a delimited plot of vegetation, with information on species cover, substrate, and other abiotic features; 2. General observations – to ensure that we capture as wide a variety of plant species found on the site as possible.

To ascertain and describe vegetation community on surveyed locations, which had associated identified plants, the plant species data, recorded as species abundance (values from 0 to 5 coded according to a DAFOR scale), were then subjected to hierarchical cluster analysis (HCA). The HCA helped define the reference classes or types of vegetation communities in the study area. Further visual observation in google earth enabled extracting additional reference points identified for roads, settlements, and other farms or agricultural areas. These, in addition to the vegetation survey points, serve as ground truth data (Table 1) for remote sensing data extraction and classification model training. The reference spatial data were extracted using reference polygons for the identified vegetation classes (described in Table A1 of online Supplementary Material).

2.2.2. Spatial Data and Processing

As spatial information, we explore Sentinel-2 (S-2) Multi-spectral Image archives provided by the European Space Agency. Available S-2 images covering the study areas were accessed using the Google Earth Engine (GEE) cloud computing platform, and the time-series images were filtered for images acquired for the area of interest in the period between 1^st January 2021 to 15^th September 2022 with cloudy pixel percentage of < 5% for each image. A total of 85 images (temporal series) were accessed following this filtering criteria, and each image collection was processed to compute 5 different vegetation indices (VIs) that provide proxy for vegetation status and cover: NDVI: Normalized Difference Vegetation Index [44], EVI: Enhanced Vegetation Index [45,46,47]. SAVI: Soil Adjusted Vegetation Index [48], MSAVI2: Modified Soil-Adjusted Vegetation Index – 2 [49] and GSAVI: Green Soil Adjusted Vegetation Index [50]. For each computed VI, we estimated and retained the max values for each pixel across the time-series [51]. Thus, the pre-processing of spatial data resulted in 5 raster images for the vegetation indices.

In addition to spatial information from VIs, we explored information lurking in the spectral signatures (for both vegetation and non-vegetation) to potentially improve the classification based on VIs. For each of the computed VI, we extracted the corresponding spectral bands to create a raster stack of spectral signatures. A total of 20 spectral bands (4 from each VI corresponding to the Blue, Green, Red, and Near Infrared bands) were subjected to a Principal Component Analysis (PCA) data reduction procedure to obtain a set of uncorrelated spatial features. Based on the PCA, 5 Principal Components (PCs) explained about 99% of the variance in the spectral signatures for the study area (Table 2), of which 4 were most significant [52]. The PCs (raster layers) estimation procedure is illustrated in Figure 3. The PCs raster layers served as additional explanatory variables for vegetation and land cover classification. Thus, in total, we used 10 raster layers (5 PCs and 5 VIs) as the spatial predictors or environmental space for vegetation classification.

2.2.3. Vegetation Type Classification and Validation

To map the vegetation assemblages (and land cover) in study landscapes, the random forest classification algorithm was implemented to delineate clusters of vegetation types, hitherto identified, and grouped as clusters of plant assemblages. The reference polygons for the vegetation assemblages or land cover classes (Table 1) were used to extract training data from the stack of raster layers (10 explanatory or predictor raster in total). Considering disparities in both the number and size of reference polygons for each class, the extracted training data or pixels were unbalanced across the considered classes.

We trained Random Forest classification models to learn distinguishing features of plant assemblages based on the ten (10) raster predictors. Considering the unbalanced tally of reference clusters and extracted spatial data between the target vegetation and land cover classes, we tuned the hyper-parameters for the random forest classification algorithm by applying forward feature selection (FFS) and spatial cross-validation procedures [16]. We used a spatial 5-fold cross-validation (k = 5), and the ground truth data had the minimum of 5 reference polygons for the “Bare rock” class; thus, ensuring each fold in cross-validation contains data from each reference class. Based on the unique IDs for each reference polygon, we implemented a spatial cross-validation, using the Leave-Location-Out Cross-Validation (LLO-CV) procedure, which specify that data from the same polygon are always grouped prevent usage in both model training and testing [16] procedure by leaving location out at each iteration of model training.

2.2.4. Classification Error and Area of Applicability

Using Dissimilarity Index (DI) in assessing the area of applicability (AOA) of the final classification map. As suggested by [16], we considered areas with predictor space (VI and Spectral features) that are unknown to the classification model are termed “Unknown spaces” and identified by high Dissimilarity Index. The dissimilarity index was computed based on weighted distance of new prediction space and the average (threshold) distance within the spatial cross-validated prediction space in the trained model. Thus, for model prediction, areas, or pixels with a DI < Threshold were classified to be within the Area of Applicability (AOA) for model prediction. The predictions within the AOA were based on the DI map and threshold (Figure 4).

2.3. Quantifying Uncertainty (Information) in Class Prediction

There is an unreliable common practice in remote sensing literature to overlook details of sampling scheme and size for accuracy assessment [18]. This is especially the case for object-based image analysis (OBIA) and less likely so for pixel-level classifications. To ascertain and delineate the contribution of spectral signature versus vegetation indices to discriminate vegetation assemblages, we compare three models or scenarios of predictor (explanatory) remote sensing variables: 1) Model VIs: Using vegetation indices (VIs) – 5 raster images, 2) Model PCs: Using principal components (PCs) of time-series spectral signatures – 5 raster images, and 3) Model VIsPCs: Using both VIs and PCs as predictor features – 10 raster images.

In quantifying the uncertainty in the pixel-based classification, a measure of uncertainty, the Shannon Entropy (H), computed from predicted class probability, was applied to evaluate information gained from the use of different spatial predictors. For pixel-based classification, the Shannon entropy or cross-entropy has been applied in remote sensing for edge detection [53], in measuring of class error [21], in quantifying land use change and dynamics [54,55,56].

The entropy of a random variable, according to information theory [57], is the average level of "information", "surprise", or "uncertainty" inherent to the variable's possible outcomes. Given a discrete random variable X, which takes values in the alphabet X and is distributed according to p: X → [0,1], the Shannon Entropy (information content, surprisal, or self-information) is estimated as follows:

$$H\left(X\right)=-\sum _{x_i\in X}{\mathrm{p}\left(x_i\right)\log }_b\mathrm{p}\left(x_i\right)$$

(1)

where: H(X) is the total amount of information in an entire probability distribution, p(x_i) is probability of classifying a pixel as class i, log is logarithm, and b is the base of the log and changing its value modifies the unit of estimated entropy. The choice of base for log, the logarithm, depends on the applications: base 2 = unit of bits or "Shannon" unit, base e = "natural units" nat, and base 10 = units of "dits" or "bans".

To ascertain the level of information or uncertainty in class predictions, the Shannon entropy was employed as a measure of the reliability in the predicted probabilities (class votes) by the random forest ensemble models. The estimated Shannon entropy (uncertainty) was compared between classes and across the different models (predictor features). To evaluate the within-class variability in uncertainty, we applied the Kernel density (distribution) profiles and non-parametric Kruskal-Wallis test with mean comparison based on Dunn test pairwise comparison and Bonferroni p-value adjustment.

2.4. Dissociating the Contribution of Predictor Features: Intensity Analysis

Intensity analysis provides a procedure for analysing land use change and transition at three different levels: time interval, category, and transition level [58,59,60]. The three successive levels of analysis form a hierarchical framework to decipher increasing details in change patterns at specific time intervals (or between maps) and assuming that the spatial extent is identical at each time point and class categories between comparable maps. To assess class uncertainty and accuracy between the three classification models, we apply the categorical and transition level intensity analysis.

The category level intensity analysis examines how the loss intensity, L_ti, from category i and the gain intensity, G_tj, to category j compares to a uniform intensity, S_t, during each map interval [M_t, M_t+1]. If L_ti < S_t, then L_ti is dormant i.e. category i experiences loss less intensively than if the change during map interval [M_t, M_t+1] were uniformly distributed across the spatial extent. If L_ti > S_t, then L_ti is active, implying that category i experiences loss more intensively than if the change during map interval [M_t, M_t+1] had uniform distribution across the analysis domain or spatial extent of map. Likewise for the gain intensity, if G_tj < S_t, then G_tj is dormant; and G_tj is active when if G_tj > S_t. The computation of S_t, L_ti and G_tj is illustrated in Eq. (3), Eq. (4) and Eq. (5), respectively.

$$U={\displaystyle \frac{\left(sizeofchangeduringintervals\right)100\mathrm{\%}}{\left(sizeofspatialextent\right)\left(durationofallintervals\right)}}={\displaystyle \frac{\sum _{t=1}^{T-1}\left\{\sum _{i=1}^J\left[\left(\sum _{j=1}^J{C}_{tij}\right)-C_{tii}\right]\right\}100\mathrm{\%}}{\left[\sum _{i=1}^J\left(\sum _{j=1}^J{C}_{tij}\right)\right]\left(M_T-M_1\right)}}$$

(2)

$$S_t={\displaystyle \frac{\left(sizeofchangeduring\left[M_t,M_{t+1}\right]\right)100\mathrm{\%}}{\left(sizeofspatialextent\right)\left(durationof\left[M_t,M_{t+1}\right]\right)}}={\displaystyle \frac{\left\{\sum _{i=1}^J\left[\left(\sum _{j=1}^J{C}_{tij}\right)-C_{tii}\right]\right\}100\mathrm{\%}}{\left[\sum _{i=1}^J\left(\sum _{j=1}^J{C}_{tij}\right)\right]\left(M_{t+1}-M_t\right)}}$$

(3)

$$L_{ti}={\displaystyle \frac{\left(sizeoflossofiduring\left[M_t,M_{t+1}\right]\right)100\mathrm{\%}}{\left(sizeofiattime{M}_t\right)\left(durationof\left[M_t,M_{t+1}\right]\right)}}={\displaystyle \frac{\left[\left(\sum _{j=1}^J{C}_{tij}\right)-C_{tii}\right]100\mathrm{\%}}{\left(\sum _{j=1}^J{C}_{tij}\right)\left(M_{t+1}-M_t\right)}}$$

(4)

$$G_{tj}={\displaystyle \frac{\left(sizeofgainofjduring\left[M_t,M_{t+1}\right]\right)100\mathrm{\%}}{\left(sizeofjat{M}_{t+1}\right)\left(durationof\left[M_t,M_{t+1}\right]\right)}}={\displaystyle \frac{\left[\left(\sum _{i=1}^J{C}_{tij}\right)-C_{tij}\right]100\mathrm{\%}}{\left(\sum _{i=1}^J{C}_{tij}\right)\left(M_{t+1}-M_t\right)}}$$

(5)

At the transition level, the analysis examines how the transition intensity, Rtij, from category i to category j compares to a uniform transition intensity Wtj given the gain of category j during interval [Mt, Mt+1]. If Rtij < Wtj, then the gain of j transitions from i less intensively during time interval [Mt, Mt+1] than if the gain of j were to have transitioned uniformly from the space that is not j at time Mt - i.e. the gain of j avoids i. Conversely, If Rtij > Wtj, this implies the gain of j transitions from i more intensively during time interval [Mt, Mt+1] than if the gain of j were to have transitioned uniformly from the space that is not j at time Mt i.e. the gain of j targets i. Equations (6) and (7) illustrate, respectively, estimates of Rtij and Wtj. The order of subscripts j and i in Ctji of the denominator of Eq. (7) allows that the summation over i subtracts category j at the initial time Mt. The relationship among the different equations and intensity estimates is illustrated in Table 3 and Table 4.

$$R_{tij}={\displaystyle \frac{\left(sizeoftransitionfromitojduring\left[M_t,M_{t+1}\right]\right)100\mathrm{\%}}{\left(sizeofiattime{Y}_t\right)\left(durationof\left[M_t,M_{t+1}\right]\right)}}={\displaystyle \frac{\left(C_{tij}\right)100\mathrm{\%}}{\left(\sum _{j=1}^j{C}_{tij}\right)\left(M_{t+1}-M_t\right)}}$$

(6)

$$W_{tj}={\displaystyle \frac{\left(sizeofgainofjduring\left[M_t,M_{t+1}\right]\right)100\mathrm{\%}}{\left(sizeofnotjattime{M}_t\right)\left(durationof\left[M_t,M_{t+1}\right]\right)}}={\displaystyle \frac{\left[\left(\sum _{i=1}^J{C}_{tij}\right)-C_{tjj}\right]100\mathrm{\%}}{\left\{\sum _{i=1}^J\left[\left(\sum _{j=1}^J{C}_{tij}\right)-C_{tji}\right]\right\}\left(M_{t+1}-M_t\right)}}$$

(7)

The conceptual organization of the various intensities and land use dynamic degrees during time intervals [M_t, M_t+1] is presented in Table 4. At the transition level, the intensity analysis determines if the gain of a category j either avoids or targets its transition from another category I; thus, the intensity analysis compares a uniform transition intensity (W_ti) to its corresponding transition intensities within each column j, R_tij. For the category level analysis, a comparison is made between change percentage in the time interval S_t and both the loss intensities (L_ti) and the gain intensities (G_tj) to determine whether each loss and gain is dormant or active.

The image classification, uncertainty analysis and intensity analyses were conducted within the R Studio Environment using several available libraries of the R programming language [62,63,64].

3. Results

3.1. Vegetation and Land Cover Classification

Plant assemblages (or communities) were discriminated, based on time-series of both spectral information and vegetation indices, with reliable accuracy. Overall, the different vegetation types were predicted with an overall accuracy of 0.834 (95% CI: 0.8178, 0.8493), a kappa of 0.699, and a global F1-Score of 0.535 as detailed on Table 5. The classification error was within reported reasonable range in remote sensing of vegetation in arid landscapes [41,65,66]. However, there was high variability in individual class accuracy across different classification models. All the models showed a high specificity and sensitivity in discriminating the monoculture farmlands and orchards (Table 6). For mixed plant assemblages, the classification uncertainty and accuracy were variable between the different models. The thematic map of 11 vegetation and land cover categories, including the unclassified areas – outside AOA, for the three classification models are depicted in Figure 5.

3.2. Classification Uncertainty

The Kernel density distribution, within the quantile range 2.5 to 97.5%, of the estimated Shannon entropy (classification uncertainty) for the vegetation categories is reported in Figure 6. A low classification uncertainty is observed for class representing mono-specific plant assemblages and/or their mosaic structure in farms or orchards (FO). For all the other considered categories, the distribution of the uncertainty is either skewed or multi-model across the different classification models (top row of Figure 6). Classification model based on vegetation indices (top middle panel in Figure 6) showed a normal distribution in the range of uncertainty (information complexity) for both uniform (FO) and high (VTW) vegetation categories. The range of statistics representing vegetation proxy and spectral information of the FO class was wide enough to reliably predict the spatial and structural variability in the class. In general, the kernel distribution of Shannon entropy (H) is indicative of high variability in classification uncertainty, for some vegetation assemblages, both in the pixel domain and local spatial levels. For the FO category, using VIs as predictor showed lowest uncertainty estimate (4.76%) in comparison to using the spectral components - PCs (6.13%) - and the combined predictors - VIsPCs (11.16%).

3.3. Change Intensity by Spectral information Versus Vegetation Indices

Overall differences in class quantity and allocation (exchange and shift) - between the classification models (maps) - are presenting in Figure 7. The exchange component represents allocation difference due to pairwise confusions while shift component estimates allocation difference from non-pairwise confusions. The class allocation components, exchange and shift, accounted most of the difference between all three thematic maps. Highest class disagreement was observed between the thematic maps from spectral components (PCs) and vegetation indices (VIs); compared to the latter, the former had lower allocation difference to map predictions based on combined predictors (VIsPCs). Although the overall accuracy of the VIs-based map was higher than that based on PCs predictors, the latter showed more consistency in discriminating other vegetation assemblages excepting farms and orchards. With reference to the combined (VIsPCs) model, the overall quantity disagreement in the PCs-based classification (7.60%) was less than that observed in the VIs -based classification (9.68%)., for the PCs to VIs map transition, the quantity component derives largely from a net gain of 2 non-vegetation categories - R and BR - and net loss in 2 vegetation categories - OVRW and OMSSS - (Figure 8). The contribution of such classes to the quantity component is backed by the heterogeneous nature and distribution of vegetation in the desert landscapes as modulated by major landforms and topography. Such spatial heterogeneity in vegetation and landforms is well captured in time-series image composite and the derived spatial features from principal (spectral) components.

The Land cover and use level change intensity between the map predictions is presented in Figure 8A. In general, we observed an active change intensity (gain and loss) for most classes across the different models. However, for some classes, large change intensity is attributable to their comparably small change area (size) in comparison (initial) map than in the reference map. Such as observed for the following classes: HSS, DSSDC, BGU, OutsideAOA. Figure 8B presents the three components of the quantified disagreement between reference and comparison maps: quantity, exchange, and shift. The components are stacked to highlight their respective contribution to overall class disagreement between the maps reference and comparison maps. In other words, their joint contribution conveys, for the purpose of this analysis, the proportion of correctly classified pixels (often calculated as 1 – total proportion of disagreement). Thus, uncertainty in class predictions based on spectral components versus vegetation indices centred around class allocation disagreements and less so from proportion correction (quantity disagreement).

4. Discussion

4.1. Vegetation Discrimination Accuracy

The classification models had accuracies in the range widely reported in remote sensing literature [18]. However, reliance on class accuracy to validate land use maps is a subject of increasing debate. Controversies in the acceptable range of classification accuracy hinge on complexities in validation related to (1) variability in mapping goals, (2) the errors of mapping as a function of the target classes, (3) heterogeneity within and between landscapes, and (4) the highly variable quality of remotely sensed data between classification experiments [68]. Thus, as a guideline to aim classification results, a generic accuracy threshold is not feasible across studies [18]. A widely used metric for accuracy assessment, the kappa statistics, which compares the observed accuracy to a baseline of random accuracy – to be properly estimated as a proportion correct for each class equalling one over the number of considered categories. This comparison is usually irrelevant and misleading as random classification is not usually the alternative scenario in mapping experiments [38]. Based on accuracy and Kappa statistics, the VI-based classification had a higher overall accuracy that would imply better classification and reliability. Similarly, an accuracy of 66% was reported using higher resolution multi-spectral imagery to map 28 vegetation classes in arid landscapes [41]. However, for classification based on VIs, the results in this study emphasize class uncertainty (or entropy) distribution and class allocation (or configuration) disagreement when compared to the classification by combination of feature selection and extraction (Figure 6 and Figure 9). Besides confirming high variability of pixel-based (soft) classification error across space, the Shannon Entropy distribution also communicates spatial uncertainty in individual class predictions. Such variability and uncertainty in class allocation is largely overlooked by mere usage and reliance on accuracy values of pixel-based land use classification.

The application of vegetation indices to spatially discriminate and map vegetation types, from landscape level to global scales, abounds in literature. The widely used vegetation indices, as those in our analysis, rely on the NIR band to distinguish specific vegetation types and to assess their status and structural changes in space and time. However, arid and hyper-arid landscapes exhibit a predominantly sparse vegetation distribution that characterised by a short duration of photosynthetic activity - largely greyer to brown than green foliage, relating to either the phases of plant senescence or adaptation strategies through variable duration of non-vegetative states. These features present challenges to the vegetation index approach of vegetation identification and discrimination in such landscapes. Thus, sporadic spatial and temporal congruence between VIs and vegetation phenological status may explain the high uncertainty and variability of the VI-based classification model (Figure 5). This is especially observable of categories with low vegetation, short-live photosynthetic activity and extensively sparse or non-vegetative state – MHCTS, OMSSS, DSSDC. As these vegetation assemblages comprise herbaceous cover, shrubs mosaics, and woodlands with comparably short-lived leafy or green states, they may convey little to no sensitivity of legacy proxies developed and tailored to capture more permanent and dense vegetation. The disconnect or link between conventional vegetation indices and different vegetation types and structure warrant discussions. For instance, research reports inconsistencies between RGB-only vegetation indices and actual crop health or NDVI [69]. Their observations were, however, based on comparison of vegetation in agricultural fields (corn versus rice). There is need for alternative, non-NIR reliant approaches of discriminating arid or rangeland vegetation for wide range applications in at landscape and other jurisdictional levels. Our results support and compound the need to reconsider vegetation mapping by tailoring and handcrafting spectral information to landscape specifics, especially in the standpoint of discriminating and monitoring the sparse and ecologically fragile vegetation in hyper-arid landscapes.

Contrary to the other vegetation and land cover categories, the class-specific accuracy for farms and orchards (FO) did not vary as much with choice of spatial predictors. Our results highlight that the VIs-based classification, with comparably high overall accuracy (Kappa of 65.91%), had high variability in class-specific prediction probability. Vegetation assemblages with low height and density showed inconsistency in reliable (low uncertainty) prediction across the classification models. This uncertain accuracy reflects temporal and spatial variability in combination of vegetation structure and phenological cycles In desert landscapes, most vegetations have distinctive spatial patterns - their type and mechanism of formation are modulated by geomorphic processes amongst other landscape processes [67]. Most of the vegetation communities are a mixture of different species, with often different growth forms, providing different functions. Thus, distinctive plant assemblages often inhabit different niches distributed across micro-habitats in a landscape. The high uncertainty and quantity disagreement observed for vegetation index-based classification model is indicative of typical vegetation structural and forms that superpose variability in vegetation index. In an arid landscape in Africa, [30] observed that spatial and temporal variations in species dominance is likely a noise contributor in modelling relationship between NDVI and biomass.

Several factors may explain the variation in classification accuracy or classification confusion (Table 6) between the considered vegetation types. Class confusion may be linked to either the status of the vegetation either at field sampling, spatial variability in class spectral and VI signature, or the resultant of information from temporal image composite - observably a mishmash within a spectrum from green to dried-out vegetation. Class error or uncertainty is attributable to the low sensitivity (or the lack thereof) of some vegetation indices on the one side. For instance, as is remarkable with vegetation classes with short-live green leaf phenology. On another note, confusion accounted by the status of the considered vegetation assemblages across temporal and spatial scales - a characteristic mixed spectral signature; these may derive spectral information or indices which are transitional between one or more classes. Uniform or transitional spectral information of vegetation or land cover may result from one or a combination of following: spatial location, the vegetation structure at field sampling and across seasons, and changes resulting from previous and on-going developmental and maintenance landscaping.

Developmental and landscaping works during periods of field surveys and beyond, all captured within by the image time-series window, results in roads and buildings of different type and status. These contribute complexities in spectral information at both the individual class levels and spatial context. Such observations may help explain the high uncertainty and miss-classification for the road and build-up areas in the study area. Besides, heterogeneity in vegetation structure has been reported to be prominent during periods of high vegetation growth [70]; however, such scenario is less intuitive in arid or desert landscapes. Colour similarities between desert plant assemblages in non-green (non-photosynthesising) state and non-vegetation areas may explain the low accuracy or misclassification using vegetation indices and the ambiguity in discriminating these classes. Variabilities in spectral information and indices are reflected in the wide distribution of class uncertainties or Shannon entropy (Figure 6). Though such temporal variabilities hold true for the vegetation classes, besides the contributed disturbance of landscaping actions to both suppress and enhance growth of some vegetation categories, the high uncertainties derive from structural and phenological heterogeneity in time and space. The period cover by time-series image composite provided adequate spatial information and feature to capturing temporal variations in both intra- and inter-annual vegetation distribution. Using vegetation proxies in assessing vegetation sensitivity to climate change, research suggests parabolic variations in along aridity index gradient and differences between vegetation types in arid and semi-arid regions [71]. We observed a low but variable class uncertainty in FO class - representing monospecific vegetation assemblages and/or their mosaic (Figures 5). For this class, the prediction uncertainty (entropy) and sensitivity of vegetation proxies is attributable to differing degrees of anthropogenic activities, predominantly the changes in structure and pattern crops and trees, irrigation and fallow regimes, and land use management.

4.2. Class Uncertainty Distribution

The classification uncertainty was estimated as the Shannon entropy of the distribution of class probability votes for each pixel. This was extrapolated to the estimated (average) uncertainty in the class predictions for different classification models (or predictor variables). The distribution of uncertainty (entropy) varied for the different land cover/use categories across the classification models. Based on the dense time series spectral information and vegetation indices at each pixel (in image composites), the estimate of uncertainty is interpreted as representative of the probable distribution error relating to the sampling for each considered vegetation assemblage or land cover class. In making inference about population statistics, sampling may provide different estimates of target population parameters. Thus, the kernel distribution may reflect, and be interpreted as, variation in the Shannon entropy or class uncertainty due to the sampling procedure - for the distribution of reference vegetation assemblages. The VI-based classification shows a narrow distribution of class uncertainty for the high vegetation categories. i.e. FO and VTW. For the other considered classes, high variability in uncertainty was observed for low vegetation and non-vegetation categories. The model, thus, reasonably captures the spatial variability in high vegetation categories (Figure 6) and less so for the low and mixed vegetation classes. The latter classes may be potential targets for sampling tailored aimed at improving classification error i.e. sampling mainly areas with feature (spectral and indices values) that are unknown or not learned by the classification algorithms. Nonetheless, the observed uncertainty distribution is not interpreted to be exhaustive for each considered classes and the study landscape - as is considered the accuracy of the classification algorithm. Though the Random forest algorithm is known to outperform other algorithms used in vegetation classification [66], there is scant or lack of literature confirming such performance for desert vegetation. Our results, while reducing this information void, highlight considerations in optimal feature selection relating to the context of analysis, vegetation specifics and expected model performance. A high variability in class-specific accuracy was observable across all considered classification models, that reflects an improvement in reported classification accuracies for similar landscapes [41,66].

The map prediction from spectral components, with comparable lower accuracy and Kappa (80.32% and 65.42% respectively) showed lower prediction uncertainty for individual classes (Figure 6). The performance of spectral components, from multispectral bands, has been reliable in land cover/use mapping based on image from different sensors and feature selection [72]; however, the measured performance has been limited to accuracy measures. While revealing the performance of spectral components over spectral indices, for the study landscape, our results compliments concerns of misleading reliance on accuracy and Kappa value in remote sensing application [38]. The classifier model with spectral components, deriving from principal component analysis of time-series spectral signatures, are representative of data dimension reduction approaches or feature selections. In the scenarios of limited sample size and spectral indices, feature selection, through Recursive Feature Elimination (RFE) improves LC classification accuracy in arid areas by removing redundant variables [65]. In our analysis, redundant spectral information was eliminated through principal component dimensionality reduction – the spectral components (PCs). A merit of multispectral data reduction, by feature selection (principal components), is the non-reliance on specific spectral bands and the reliability in differentiating both natural (including vegetation) and artificial features [72]. By contrast, application of vegetation indices, as a feature extraction approach is constrained by a select of specific spectral bands to match analysis objectives and target features. The wavelength range of spectral bands may vary when merging images captured by different sensors. Thus, a select of VI proxies will often provide a constrained or limited sensitivity range in relation to specific landscape features and their variability in space and time.

4.3. Class Spatial Allocation and Intensity

To resolve the contribution of vegetation indices versus spectral signature (principal) components in class predictions, intensity analysis was applied to assess the spatial configuration of vegetation classes or land cover [58,59]. This enabled quantifying the class change and dynamics between the classification models and the underlying change patterns and intensity related to the spatial predictors. The marginal difference in overall accuracies between the classification models supports direct comparison of class prediction and spatial allocation in the respective thematic maps (Figure 5). The change pattern in the predicted class supports our understanding of the contribution of the different predictor features, or their interaction, to pixel class at different spatial scales.

Overall, the change in each map interval was accounted mostly by allocation differences due to class exchange. The VIs-VIsPCs map interval had the smallest quantity differences but the largest proportion of allocation (exchange and shift) differences (Figure 7). Most of this quantity difference is contributed by FO and OVRW categories (Figure 8B); while the former category has an overall quantity gain, the latter is a net quantity loss. Each bar in Figure 8B demonstrates that less than a quarter of the overall disagreement is attributable to quantity disagreement. Quantity disagreement refers the amount of difference between the reference map and a comparison map that is due to the less than perfect match in the proportions of the categories. An objective in this study was to estimate and compare the contribution of predictor variables to land use/cover change based on the three classification models; in this regard, class allocation transition due to class exchange and shift disagreements were of much less significance. Exchange disagreement consists of a transition from category i to category j in some pixels and a transition from category j to category i in an identical number of other pixels; meanwhile, Shift refers to the difference after subtracting quantity and exchange from the overall. Thus, although the overall classification error is in the range 17.29% to 19.68% across the classification models, this compares with merely reported accuracies for land use/cover mapping in similar landscapes [41,66]. As observable in Figure 8, pairs of categories may contribute to exchange disagreement between maps. Thus, it is useful to separate overall change into its constituent components given that each component is interpretable in relation to application specifics. For instance, in change characterization, error assessment and map comparison. Insights into the effect of category aggregation on model performance can help eliminate less contributing or informative factor to the model performance. The observed performance in map comparisons and the contributions to class quantity and allocation transition (Figure 8) suggest that spectral information may provide optimal spatial features in discriminating and mapping vegetation assemblages in arid and hyper-arid landscapes. The results are reliable for both the purpose of analysis and as a reference distribution map (Figure 9) for different intended practical applications.

5. Conclusions

The need for monitoring and management of vegetation assemblages in desert and arid landscapes have gained country-specific and international attention due to spreading developmental actions that impact the ecological heritage in these landscapes, besides ubiquitous ecological challenges elicited by erratic climate patterns. In this study, we question and quantify the uncertainty associated with, of contributed by, using vegetation indices versus spectral information, from Sentinel-2 imagery, as proxies for mapping vegetation communities of drylands. The results suggest considerations in translating multi-spectral remote sensing data to context-specific vegetation mapping, with emphasis on the ephemeral and ecologically erratic vegetation distribution in desert landscapes. Our findings underscore the distribution of classification uncertainty and spatial difference (quantity and spatial allocation) related to the choice of predictor features between spectral information and vegetation indices.

The discrimination and mapping of plant species distribution on drylands, modulated by environmental variables such as soil and topography, depends on the choice of spatial predictor variables or metrics, besides other parameters relating to model selection, desired level of map details, complexities in the target landscapes, etc. The choice of predictor variables constitutes persistent burden and challenges to resource managers, who require informed guidance relating to available reference data, the scale of analysis, and to reliably match models and/or knowledge with context. In context of remote sensing of drylands, vegetation types are widely mapped according to a select metrics of vegetation and environmental differentiation. By questioning the reliability of classic metrics of vegetation across different landscape, this study suggests considerations of feature extraction versus feature selection in mapping and quantifying vegetation communities in arid landscapes. Thus, with consideration of landscape specifics, the findings highlight a rethinking in a) combining predictors for different features for vegetation mapping and b) using information loss (entropy) to decipher spatial uncertainties in class allocation. Though the resultant maps are not perceived as potential inputs for models that seek to maximize overall accuracy, the results contribute relevant insights into the characterisation and mapping of desert vegetation communities and woodlands by integrating remote-sensing data.

The analysis and quantification of class uncertainty may lend a-priori information in understanding the application and limits of remotely sensed data to landscape specifics. In quantifying the level of uncertainty in mapping vegetation assemblages, using different spatial predictors, the presented analyses fall short of considering an exhaustive number of possible proxies in identifying and mapping the considered vegetation communities. Nonetheless, this study contributes insights into using and optimising feature extraction (VIs) versus feature selection (PCs) for pixel-based classification of arid vegetation communities. These results are perceived to leverage usage of Sentinel-2 satellite imagery to inform and support the interpretation of vegetation maps intended for operational management and developmental planning in hyper-arid landscapes..