Improved Mapping of Mountain Shrublands Using the Sentinel-2 Red-Edge Band

Shrub encroachment into grassland and rocky habitats is a noticeable land cover change currently underway in temperate mountains and is a matter of concern for the sustainable management of mountain biodiversity. Current land cover products tend to underestimate the extent of mountain shrublands dominated by Ericaceae (Vaccinium spp. (species) and Rhododendron ferrugineum). In addition, mountain shrubs are often confounded with grasslands. Here, we examined the potential of anthocyanin-responsive vegetation indices to provide more accurate maps of mountain shrublands in a mountain range located in the French Alps. We relied on the multi-spectral instrument onboard the Sentinel-2A and 2B satellites and the availability of red-edge bands to calculate a Normalized Anthocyanin Reflectance Index (NARI). We used this index to quantify the autumn accumulation of anthocyanin in canopies dominated by Vaccinium spp. and Rhododendron ferrugineum and compared the effectiveness of NARI to Normalized Difference Vegetation Index (NDVI) as a basis for shrubland mapping. Photointerpretation of high-resolution aerial imagery, intensive field campaigns, and floristic surveys provided complementary data to calibrate and evaluate model performance. The proposed NARI-based model performed better than the NDVI-based model with an area under the curve (AUC) of 0.92 against 0.58. Validation of shrub cover maps based on NARI resulted in a Kappa coefficient of 0.67, which outperformed existing land cover products and resulted in a ten-fold increase in estimated area occupied by Ericaceae-dominated shrublands. We conclude that the Sentinel-2 red-edge band provides novel opportunities to detect seasonal anthocyanin accumulation in plant canopies and discuss the potential of our method to quantify long-term dynamics of shrublands in alpine and arctic contexts.


Introduction
In the context of recently observed shifts in climate [1] and vegetation [2,3] in the European Alps, efficient mapping of mountain land cover (LC) over time and space using satellite imagery has emerged as a high priority for environmental monitoring and management [4][5][6][7].Several studies and projects have proposed LC maps covering the Alps based on supervised classification of ground data and available satellite imagery {CORINE Land Cover [8]; Occupation des Sols (OSO) Land Cover [9]}.Automated mapping in mountainous areas is notoriously difficult due to the effects of topography, which leads to fine-scale variation in LC types [9] as well as changes in illumination due to slope aspect and angle [10].Mapping low evergreen and deciduous shrubs (< 1 m in height) in above treeline habitats represents a challenge, given that it is often difficult to distinguish shrubland from nearby grassland using classic multispectral vegetation indices [e.g., Normalized Difference Vegetation Index (NDVI)] or aerial photos.
Shrub expansion into grassland and tundra during recent decades is a widespread phenomenon that has been reported in both arctic and alpine environments in response to warming temperatures [11], changes in snow cover duration, and shifts in land-use [12] or grazing pressure [13].In addition to climate change, the European Alps have experienced widespread land abandonment during recent decades, contributing to increased woody plant cover in subalpine pastures [14][15][16].Cannone et al. [17] reported an increase in shrub cover of 5.6% per decade between 2400 m and 2500 m above sea level (a.s.l) starting from 1953 to 2003 in the Italian Alps.More recently, in the French Alps, Carlson et al. [2] reported an increase in plant photosynthetic activity between 1984 and 2015 in cliff, scree, and grassland habitats and hypothesized a contribution of shrub dynamics to observed greening trends between 2200 m and 2500 m a.s.l.Replacement of species-rich grasslands by shrubs can drastically reduce local biodiversity [18] and strongly impact the quality of available pasture resource [19].While there is growing evidence that shrub expansion is underway in the Alps, landscape scale analysis of shrub distribution and dynamics is currently impeded by our inability to systematically distinguish shrub cover from other forms of vegetation in above treeline habitats using widely available satellite imagery.From both scientific and land management perspectives, there is a clear need to improve remote sensing-based methods for mapping shrub cover.
The accumulation of red pigment in certain shrub species at the start of the growing season and especially during the fall represents a promising opportunity to distinguish shrub cover from surrounding grasslands.This process is common in many regions, from tropical [20] to arctic [21] biomes.In temperate alpine zones, the accumulation of red pigment anthocyanins is responsible for visible reddening of heath vegetation [22].A high relative content of anthocyanins in leaves is partly caused by a combination of low temperatures and excessive irradiance [23,24], which occurs in deciduous plants at the beginning of the growing season during late spring (May-June) in juvenile leaves and at the end of the season in senescing leaves during the fall (late September-November) for temperate alpine zones.Due to the same mechanism, evergreen species can also exhibit leaf reddening during the same periods or throughout the winter [25,26].In deciduous Ericaceae shrubs, including Vaccinium myrtillus and Vaccinium uliginosum, leaves remain red from fall to foliar abscission or spring due to decreases in chlorophyll and an increase in anthocyanins in leaf tissue.For evergreen Ericaceae shrubs such as Rhododendron ferrugineum, Vaccinium vitis-idaea, and Calluna vulgaris, the reddening of leaves during the fall and the winter is less visible due to maintained chlorophyll levels.Nonetheless, we hypothesize that the accumulation of anthocyanins in the leaf tissues of these species will modify reflectance values during fall months.In the case of non-Ericaceae evergreen shrubs such as Juniperus communis or Juniperus sibirica, we do not expect to be able to distinguish shrub cover from grassland due to hypothesized reduced levels of anthocyanin accumulation.In the Alps, Ericaceae species that redden in late fall season represent dominant species for many shrublands ecosystems above the treeline [27], making this region well-suited for testing and applying novel methods aimed at mapping Ericaceae-dominated shrublands.
Multiple studies have demonstrated the use of imaging spectroscopy to estimate anthocyanin content in leaf tissues [28][29][30][31][32][33][34].Based on the spectral signal of leaves containing a high concentration of anthocyanins, Gitelson et al. [28] proposed the Anthocyanin Reflectance Index (ARI) based on the green and red-edge spectral responses of leaves.Gitelson et al. [28] found that reciprocal reflectance near the green band 550 nm (R550) −1 accurately measured chlorophyll content in leaves with minute amounts of anthocyanins but further demonstrated that green reflectance values were disturbed in leaves containing high concentrations of anthocyanins.In addition, Gitelson et al. [28] found that reciprocal reflectance near 700 nm (R700) −1 was closely related to chlorophyll content in leaves regardless of anthocyanin concentration.Accordingly, anthocyanin content was found to be proportional to the difference (R550) −1 -(R700) −1 .A modified version of ARI was proposed to minimize the dependency on leaf thickness and density; however, preliminary analyses in our study system showed no significant differences between the two models.Accordingly, we chose to focus on ARI as the basis for our analysis.Our work extends the results of Gitelson et al. [28], derived from leaf spectroscopy, to the canopy level, with the aim of quantifying spatial variation in anthocyanin concentrations in the context of heterogeneous mountain plant cover.
Although ARI appears to be effective for non-destructive estimation of leaf anthocyanin content [34], due to the necessity of using a red-edge band, ARI cannot be calculated using widely available multispectral optical satellite sensors (e.g., Landsat, SPOT).The Medium Resolution Imaging Spectrometer (MERIS), one of the payloads on the European Space Agency's Envisat, or Sentinel-3 [35], a medium resolution satellite that is the continuation of the Envisat mission, both have a red-edge band, and thus ARI could be calculated.However, the spatial heterogeneity of mountain land cover necessitates a high spatial resolution.The five day revisit time and the relatively high spatial resolution of Sentinel-2 (10-20m) combined with thirteen available spectral bands ranging from visible to short-wave red provide novel opportunities for observation of vegetation dynamics in temperate mountain ecosystems [36].
In this paper, we examined the efficiency of a modified ARI index for detecting and mapping Ericaceae-dominated shrublands in the Belledonne massif located in the French Alps.Our primary goal was to test and validate a novel method aimed at distinguishing Ericaceae-dominated shrublands from adjacent mountain grasslands.Our objectives were to (i) identify the most appropriate period of the year for quantifying leaf reddening utilizing three years of Sentinel-2 imagery, (ii) evaluate the performance of the Normalized Anthocyanin Reflectance Index (NARI)-based method using photo-interpretation and floristic data, and (iii) compare our output maps to existing land cover products.

Study Area
The study area is located in the Belledonne and the Grandes Rousses massifs in the French Alps (Figure 1).The elevation gradient ranges from 250 m above sea level (a.s.l.) at the valley floor to 3442 m a.s.l.The surface area is 1377 km 2 , of which 536 km 2 is composed of above treeline vegetation (above 1600 m and under 2500 m a.s.l.).The Belledonne and the Grandes Rousses massifs are predominantly made of acidic siliceous rocks and exhibit large areas of Ericaceae-dominated shrublands above the treeline.We identified three test zones along a north-south gradient (Figure 1).

Remote Sensing Data
Sentinel-2 is an Earth observation mission composed of two identical satellites, Sentinel-2A and Sentinel-2B, respectively launched on 23 June 2015 and 7 March 2017.The 5 day revisit time available since March 2017 was combined with 10, 20, or 60 m spatial resolution of spectral bands.We acquired 84 Sentinel-2A and 39 Sentinel-2B (collectively referred to hereafter as S2) optical images for the T31TGL tile from January 2016 to December 2018 in order to compute NARI and NDVI time series for grasslands and shrublands during a 3 year period (N = 123 scenes).From this collection, 11 images from mid-summer to late fall of 2018 were retained to compute further analyses.The multispectral S2 images are composed of 4 spectral bands (B2: 490 nm, B3: 560 nm, B4: 665 nm, B8: 842 nm) at 10 m spatial resolution and 6 spectral bands at 20 m from the Red Edge (B5: 705 nm, B6: 740 nm, and B7: 783 nm) to the NIR and the SWIR domains (B8a: 865 nm, B11: 1610 nm, B12: 2190 nm).Data were downloaded from the French national THEIA platform (http://www.theialand.fr) at level 2A (i.e., orthorectified product in surface reflectance) with clouds and clouds shadows masks provided at 10 and 20 m spatial resolutions.The clouds, the cloud shadows masks, and the atmospheric correction were performed using MAJA (Multi-sensor Atmospheric Correction and Cloud Screening-Atmospheric and Topographic CORrection, MACCS-ATCOR Joint Algorithm).This algorithm includes a correction of atmospheric effects by combining multi-temporal and multi-spectral criteria to estimate the aerosol optical thickness [37] and also applies a topographic correction to account for changes in illumination caused by slope

Remote Sensing Data
Sentinel-2 is an Earth observation mission composed of two identical satellites, Sentinel-2A and Sentinel-2B, respectively launched on 23 June 2015 and 7 March 2017.The 5 day revisit time available since March 2017 was combined with 10, 20, or 60 m spatial resolution of spectral bands.We acquired 84 Sentinel-2A and 39 Sentinel-2B (collectively referred to hereafter as S2) optical images for the T31TGL tile from January 2016 to December 2018 in order to compute NARI and NDVI time series for grasslands and shrublands during a 3 year period (N = 123 scenes).From this collection, 11 images from mid-summer to late fall of 2018 were retained to compute further analyses.The multispectral S2 images are composed of 4 spectral bands (B2: 490 nm, B3: 560 nm, B4: 665 nm, B8: 842 nm) at 10 m spatial resolution and 6 spectral bands at 20 m from the Red Edge (B5: 705 nm, B6: 740 nm, and B7: 783 nm) to the NIR and the SWIR domains (B8a: 865 nm, B11: 1610 nm, B12: 2190 nm).Data were downloaded from the French national THEIA platform (http://www.theialand.fr) at level 2A (i.e., orthorectified product in surface reflectance) with clouds and clouds shadows masks provided at 10 and 20 m spatial resolutions.The clouds, the cloud shadows masks, and the atmospheric correction were performed using MAJA (Multi-sensor Atmospheric Correction and Cloud Screening-Atmospheric and Topographic CORrection, MACCS-ATCOR Joint Algorithm).This algorithm includes a correction of atmospheric effects by combining multi-temporal and multi-spectral criteria to estimate the aerosol optical thickness [37] and also applies a topographic correction to account for changes in illumination caused by slope angle and aspect [38] in order to provide "flat reflectance".A multi-temporal approach (MAJA) has been preferred to multi-spectral (e.g., Sen2cor), as it outperforms it when highly repetitive observations with constant viewing angles are available [39].
To calculate vegetation indices, five bands were used: three at 10 m resolution (B3, B4, and B8) and two at 20 m resolution (B5 and B11).The 20 m resolution bands were resampled to 10 meters using bilinear interpolation.Bands values were converted from 64 bit signed integer format to reflectance varying from 0 to 1.We also added a constant (0.05) to visible bands B3, B4, and B5 to correct for noise due to surface reflectance correction (comm.pers.O. Hagolle, https://github.com/olivierhagolle).For each atmospheric and terrain corrected image from 2016 to 2018, we calculated (i) the Normalized Anthocyanin Reflectance Index (NARI), which is a normalization of the ARI index presented in Gitelson et al. (2004; ARI = 1/R550 -1/R700) and is sensitive to plant canopy anthocyanin content; (ii) the Normalized Difference Vegetation Index (NDVI) [40], as it has been widely used for land cover classification (e.g., [41]),.In addition, we estimated the (iii) Normalized Difference Snow Index (NDSI) to identify the presence or the absence of snow cover in a pixel based on a threshold of 0.4 [42].While this method is known to be inefficient in forested areas [43], it was appropriate to use in our case, as our studied ecosystems are located above treeline.NARI, NDVI, and NDSI were calculated as follows: where R is the reflectance of the green (B3), the red (B4), the red-edge (B5), the near-infrared (B8), and the shortwave-infrared (B11) channels.Spectral responses of vegetation communities identified in the field as pure end-member pixels of Ericaceae-dominated shrubland and grassland communities were assessed using a hyperspectral image acquired on 28 September 2015 at the Col de Balme (1950 m a.s.l.) in the nearby Mont-Blanc massif (Project ISOTHERM).The image was acquired using a Fenix hyperspectral camera and consisted of 622 bands ranging from 400 to 2500 nm with an average bandwidth of 30 nm.Raw imagery was processed using the Airborne Processing Library (APL) software and delivered at level 1b (georeferenced and radiometrically corrected).

Snow, Cloud, and Topographic Masks
For each scene, we applied masks to exclude snow-covered, shadowed, and cloud-covered pixels from our analysis.Cloud masks computed for each date were provided by THEIA.We computed a hillshade mask using a 25 meter digital elevation model [DEM, European Digital Elevation Model (EU-DEM), version 1.0] resampled to 10 m using bilinear interpolation to align with S2 bands.Slope and aspect were computed using the Horn algorithm [44] implemented in R using the terrain function in the raster package [45].Hillshade was calculated for each scene by combining solar azimuth and elevation (computed as 90 -zenith), which were used as inputs in the hillShade function in the raster package (based on Horn [44]).Mask values ranged from 0 to 1 and represented illumination, with a value of 0 representing complete darkness.Given that poorly illuminated pixels can lead to erroneous reflectance values, we applied a manually calibrated threshold to generate a shadow mask (pixels with an illumination value < 0.25 were excluded).Snow cover maps derived from NDSI (as described above) were used for snow masking.
We calculated seasonal variations in the percentage of snow cover extent, cloud, and shadow masks for all scenes over the three years.A cubic smooth spline (smoothing parameter of 0.4) was fitted to the data in order to estimate the daily change in snow cover extent and cloud cover and to quantify surface loss due to snow and cloud cover over the three years.Pixels of the resulting classification below 1600 m a.s.l. and with NDVI < 0.2 (calculated using 08-28-2018 as a reference) were excluded from the analysis, as they corresponded to below treeline vegetation or sparsely vegetated areas.Finally, water surfaces were also masked using a 10 meter buffer mask.

Anthocyanin Mapping Protocol
After applying masks, we implemented three approaches to discriminate shrublands from grasslands using NARI.First, we assessed trends in NARI from mid-summer, when anthocyanins levels are expected to be low (hereafter NARI-low date; 19 July 2018), to fall, when anthocyanins levels are expected to be higher (hereafter NARI-high date; 27 September 2018).We utilized a total of six dates and applied a median-based linear model (mblm) using the Theil-Sen single median algorithm [46] fitted between NARI and date for unmasked pixels.A pixel was excluded from the analysis if there were fewer than three available NARI values in the time series.The Theil-Sen procedure is a rank-based test that calculates the non-parametric slope and the intercept of time-series [46].It is known to be robust for small, heterogeneous datasets containing outliers.We used pixel-level slope coefficients to calibrate a classification threshold for identifying pixels dominated by Ericaceae shrubs, hereafter referred as the "mblm" method.Second, we calculated a difference between NARI-low and NARI-high dates to quantify change in anthocyanin accumulation, hereafter referred as the "delta" method.Third, we directly used the index value of the NARI-high date, hereafter referred as the "Monodate" method.Finally, we applied the same three approaches to time series of NDVI in order to compare the effectiveness of the two indices for detecting shrub cover.Threshold values used for classification were obtained using receiver operating characteristic (ROC) curves presented in the Model Validation section.

Photo-Interpretation and Ground Reference Data
The six models were evaluated using two independent sources of data.First, ten polygons of each cover type were extracted from the CarHAB mapping program [47].These polygons were generated using a segmentation approach applied to a high-resolution satellite image and were subsequently validated by field observations of vegetation [48].Second, we digitized a dataset of 294 polygons using photointerpretation of a high-resolution image (Google, DigitalGlobe) from 10 April 2016 in Google Earth Pro during the anthocyanin peak.This image enabled reliable interpretation of Ericaceae-dominated shrublands due to the leaf reddening caused by anthocyanin production (Figure S1).In total, we analyzed 52 grass and 242 shrub polygons, respectively, representing 269 hectares of grassland with an average patch size of 5 hectares (±9.5 ha) and a mean distance between polygon centroids of 550 m (±516 m) and 234 hectares of Ericaceae-dominated shrubland with an average patch size of 1 hectare (±2.4 ha) and a mean distance between polygon centroids of 137 m (±100 m).As average distance between ground samples was superior to average patch size, we assumed spatial independence between samples.We used 70% of the grassland and the shrubland polygons for model calibration and used the remaining 30% of the dataset for model validation, which represented, on average, 57 hectares of grasslands and 52 hectares of shrubland available for validation.
The second source of independent observations consisted of floristic surveys carried out by the National Alpine Botanical Conservatory (CBNA).The dataset included 347 10 x 10 m community plots containing 451 vascular plant species sampled between 2008 and 2016, half of which were sampled in 2016.Each community plot contained a visual estimate of species cover using Braun-Blanquet [49] classes.(1: less than 1%; 2:1 to 5%; 3:5 to 25%; 4:25 to 50%; 5:50 to 75%; 6: up to 75%).We converted these values to relative abundance using mean cover class percentages.Within sampling plots, an expert botanist identified vascular plants and the relative abundance of each species.We identified 174 plots as shrublands, which were characterized by a relative abundance above 50% of target Ericaceae species Vaccinium myrtillus, Vaccinium vitis-idea, Vaccinium uliginosum microphyllum, and Rhododendron ferrurigeum.We identified 173 plots as grasslands, which contained a relative abundance above 25% of target perennial grass species Nardus stricta and Patzkea paniculata and below 1% of Ericaceae shrubs.
Due to potential signal distortion caused by rock or tree and uncertainties in plot location (Figure S2 shows Google Earth images of two typical cases of distorted signal due to surrounding environmental conditions), the 347 plots were visually verified using very high-resolution images from Google Earth Pro in order to select above treeline plots with homogenous land cover.Furthermore, plots containing more than 5% cover of coniferous tree species Pinus cembra, Pinus mugo, and Picea abies were identified as open forest and were excluded from the analysis.Thus, from the 174 plots initially assigned as shrublands, 80 were in open or closed forest, 36 in screes, and 28 masked by clouds, resulting in 30 plots of homogenous above treeline shrub cover.For the 173 plots assigned as grasslands, 21 were in forests, 23 masked, and 59 in screes, leaving a remaining 70 plots for analysis.For each plot, model values were extracted by averaging values extracted using a 10 m buffer around the plot coordinates.Shrubland plots have a mean elevation of 1861 ± 148 m and are mainly composed of Ericaceae species, which are expected to exhibit late fall reddening and grasses that are common beside shrubs (e.g., Nardus stricta).Grassland plots have a mean elevation of 1921 ± 226 m and mostly include above tree-line perennial graminoids that are commonly associated with shrublands and thus are often confused with shrublands in widely used land cover products (e.g., Nardus stricta or Patzkea paniculata).Results for both groups were compared using the Wilcoxon Signed-Rank test.

Evaluation of Model Performance
We used ROC curves to test the performance of all six approaches (three for NDVI and three for NARI) using the pROC package in R [50].ROC curves allowed us to identify the optimal approach with the highest success rate for distinguishing between shrublands and grasslands.The three approaches, described as "mblm", "Delta", and "Monodate", were applied to NARI and NDVI time series and were used as explanatory variables with respect to shrubland and grassland field plots identified above.To identify the best method, we first computed the area under the curve (AUC) and tested for significant differences from 0.5, i.e., randomness, using the Mann-Whitney test and following the process outlined in Manson and Graham [51,52].Second, we assessed whether or not the sensitivity of ROC curves was statistically different at a given level of specificity corresponding to a 5% chance of committing a commission error (0.95, α = 0.05) using a bootstrap method with 2000 replicates as used in Pepe et al. [53].Third, as the number of masked pixels varied based on the method used, we compared the surface area loss due to cloud, snow, and shadow masks over the study area and for the six mapping approaches.We then retained the model minimizing surfaces loss and maximizing AUC and sensitivity and applied the optimized threshold to discriminate between shrublands and grasslands.
In order to assess the quality of the produced shrub cover maps, we compared the classification to the remaining 30% of ground reference data.Model performance metrics were calculated using 100 repetitions, during which one pixel was randomly selected within each ground-reference polygon to avoid spatial dependencies.Furthermore, to avoid edge effects, only pixels that completely overlapped polygons were retained.Because polygons with small surfaces (inferior to 0.1 ha) or a Polsby-popper measure [54] too low (inferior to 0.25) did not allow for multiple non-autocorrelated or completely overlapped pixels, 18 polygons were discarded from the quality assessment, leaving a total of 71 pixels per run.For each confusion matrix, we calculated several indices, including (i) Recall, which is the fraction of correctly classified pixels with regard to photo-interpreted polygons, (ii) Precision, which represents the fraction of correctly classified pixels with regard to all pixels classified as such in the image, (iii) Omission error rate, which represents the fraction of pixels incorrectly classified as negative, (iv) Commission error rate, which represents the fraction of pixels wrongly classified as positive, (v) Fscore (also known as F-1 score or F-measure), which is the harmonic mean of Recall and Precision, and (vi) the Kappa coefficient (K), which expresses a relative difference between the total accuracy compared to random accuracy based on an arbitrary threshold value.
(ii) Precision = TP / (TP+FP), where TP is True Positive, FP is False Positive, TN is True Negative, TP is True Positive, TA is Total Accuracy, and RA is Random Accuracy.According to Landis and Koch [55], a Kappa coefficient of 0.6 or more can be considered as acceptable.

Land Cover Comparisons
Three existing LC products were compared to our classification: (i) two classes of the 100 meter resolution 2018 CORINE Land Cover (CLC) available for Europe [8], transitional woodland-shrub (324) and moors and heathland (322), which we considered to be the equivalent of shrublands in our classification, (ii) the mixed tree-shrub-herbaceous class (100) of the 2015 Global European Space Agency Climate Change Initiative (CCI) Land Cover map at a spatial resolution of 300 m [56], and (iii) the woody moorland (36) class of the 10 meter resolution 2017 Occupation des Sols (OSO) LC map, which is only available for France [9].In order to determine the relevance of comparison with our mapping results, we evaluated the ability of coarse resolution LC products (CLC and CCI) to detect shrub patches.We first converted our raster classification to polygons, then selected shrub patches with surfaces superior to 1 hectare from our classification, leading to a total of 760 polygons.For each pixel of 100 and 300 m of spatial resolution, we measured the fraction covered by shrub polygons in order to select pixels that were dominantly covered by shrubs at the desired spatial resolution.We defined shrub-dominated pixels as those covered by 75% or more of shrub patches and used these pixels as a basis of comparison with land cover maps (Figure S3 provides further details regarding this analysis).This approach resulted in 905, 4046, and 439,073 shrubland pixels as a basis of comparison with CCI, CLC, and OSO products.Our classification was resampled to coarser spatial resolution for CLC and CCI comparisons (respectively, 100 and 300 m).For OSO, CLC, and CCI, to avoid spatial dependencies, approximately 10% of total pixels were randomly sampled with a distance criterion based on nonaligned systematic sampling (package sp in R, [57,58]).We repeated the process 100 times and produced a confusion matrix and resulting statistical analyses.Because calculating flat pixel-based surface area ignores the effects of topography and surface roughness and leads to major underestimation of features in mountainous areas, we used a triangulated irregular network (TIN) to calculate surface areas accounting for topography.

Spectral Responses
Comparison between end-member pixels corresponding to leaves with and without anthocyanins (respectively dominated by Vaccinium uligonosum and Nardus stricta) demonstrated coherent results with respect to Gittelson et al. [27].Reflectance in the green band (550 nm) was lower in the case of the shrub canopy, while no substantial difference was observed in the red-edge spectrum (700 nm), suggesting similar chlorophyll content in the two canopies (Figure 2).NARI and NDVI values for the end-members pixels during the late fall yielded similar values to the results presented in Figure 3 (NDVI = 0.46 and 0.44; NARI = 0.18 and 0.25 for N. stricta and V. uliginosum endmembers, respectively).

Time Series
NDVI dynamics extracted for photo-interpreted polygons of shrublands and grasslands were consistent across the three study years (Figure 3, top panels).Over the course of the growing season, shrublands exhibited slightly lower NDVI values but with similar dynamics compared to grasslands.In contrast, NARI dynamics varied throughout year between grasslands and shrublands (Figure 3, mid panels).In 2016 and 2018, NARI rapidly increased in shrublands starting in early September, peaking in late September and declining between mid to late October.The beginning of the growing season was characterized by high values of NARI for shrublands and slightly lower values for grasslands for the year 2018.Similar trends were not observed for 2016 and 2017 due to high cloud cover leading to pixels obstruction at key dates and an insufficient number of scenes.In contrast to 2016, 2017 and 2018 were characterized by frequent cloud cover during late fall leading to a lack of information (respectively 55% and 60% of surfaces loss at supposed peak due to combined masks for 2017 and 2018, and only 27% for 2016).Nonetheless, for 2018, the supplementary images obtained from Sentinel-2B enabled capturing the late fall reddening of leaves.During the early growing season, high obstruction was observed due to combination of snow and clouds, while obstruction toward the end of the growing season was mostly due to increasing shadows on northern slopes.Four dates with low obstruction from clouds and hillshade were available for 2016 (

Time Series
NDVI dynamics extracted for photo-interpreted polygons of shrublands and grasslands were consistent across the three study years (Figure 3, top panels).Over the course of the growing season, shrublands exhibited slightly lower NDVI values but with similar dynamics compared to grasslands.In contrast, NARI dynamics varied throughout year between grasslands and shrublands (Figure 3, mid panels).In 2016 and 2018, NARI rapidly increased in shrublands starting in early September, peaking in late September and declining between mid to late October.The beginning of the growing season was characterized by high values of NARI for shrublands and slightly lower values for grasslands for the year 2018.Similar trends were not observed for 2016 and 2017 due to high cloud cover leading to pixels obstruction at key dates and an insufficient number of scenes.In contrast to 2016, 2017 and 2018 were characterized by frequent cloud cover during late fall leading to a lack of information (respectively 55% and 60% of surfaces loss at supposed peak due to combined masks for 2017 and 2018, and only 27% for 2016).Nonetheless, for 2018, the supplementary images obtained from Sentinel-2B enabled capturing the late fall reddening of leaves.During the early growing season, high obstruction was observed due to combination of snow and clouds, while obstruction toward the end of the growing season was mostly due to increasing shadows on northern slopes.Four dates with low obstruction from clouds and hillshade were available for 2016 (

Best Model and Classification Validation
NDVI-based models showed low AUC (0.61 ± 0.01 on average) and were not significantly different from a random model (Probability-value > 0.05; Table 1).Also, sensitivity was close to zero for Delta and mblm NDVI models and was 0.23 for the Monodate approach.By contrast, NARI-based models show high AUC (0.92 ± 0.01 on average) with significant difference from a random model (P-val < 1e16), meaning that NARI-based models are more efficient to discriminate shrublands and grasslands (Figure 4).The NARI-based Delta model showed the best results with an AUC of 0.98 and a sensitivity of 0.93, but because this method used only two dates out of the six available, it showed 19.44% of surface loss due to cloud and hillshade.Conversely, the NARI-based mblm model showed low surface loss of 8.35% because the model was calibrated using up to six dates, depending on local cloud cover (Table 1).The remaining surface loss was mainly due to shadows on northern slopes.No significant differences were found between Delta and mblm NARI-based models.Accordingly, the NARI-based mblm model was retained with an optimal classification threshold of 0.0004, i.e., pixels with slope values superior to the threshold were classified as shrublands.

Best Model and Classification Validation
NDVI-based models showed low AUC (0.61 ± 0.01 on average) and were not significantly different from a random model (Probability-value > 0.05; Table 1).Also, sensitivity was close to zero for Delta and mblm NDVI models and was 0.23 for the Monodate approach.By contrast, NARI-based models show high AUC (0.92 ± 0.01 on average) with significant difference from a random model (P-val < 1e16), meaning that NARI-based models are more efficient to discriminate shrublands and grasslands (Figure 4).The NARI-based Delta model showed the best results with an AUC of 0.98 and a sensitivity of 0.93, but because this method used only two dates out of the six available, it showed 19.44% of surface loss due to cloud and hillshade.Conversely, the NARI-based mblm model showed low surface loss of 8.35% because the model was calibrated using up to six dates, depending on local cloud cover (Table 1).The remaining surface loss was mainly due to shadows on northern slopes.No significant differences were found between Delta and mblm NARI-based models.Accordingly, the NARI-based mblm model was retained with an optimal classification threshold of 0.0004, i.e., pixels with slope values superior to the threshold were classified as shrublands.Table 1.Statistical results from the six models comparisons, including: area under the curve (AUC) presented with 95% confidence intervals, P-value correspondence (*** : <0.001; NS: not significant), the identified statistically optimal threshold, the sensitivity (with specificity fixed as 0.95, α = 0.05) for the selected threshold, and the percentage of masked pixels in the final land cover (LC) map.For the NARI mblm model, pixels with values above the optimal threshold were classified as shrublands, and the rest were classified as grasslands.Based on confusion matrices calculated using the 30% remaining polygons of the ground reference dataset (Table 2), we obtained excellent classification results (Figure 5).The Fscore was 0.90 ± 0.02 with a recall of 0.84 ± 0.03 and a precision of 0.97 ± 0.02, which denotes very low misclassification.Sensitivity and specificity were both high (0.93 ± 0.05 and 0.84 ± 0.03, respectively), indicating a low number of commission and omission errors (respectively, 0.15 ± 0.03 and 0.06 ± 0.05), which means that our classifier was missing approximatively 15% of the shrublands in our study site and was misclassifying approximatively 6% of the grasslands.The Kappa coefficient was 0.67 ± 0.05, which denotes a good classification rate according to Landis and Koch (1977).Table 1.Statistical results from the six models comparisons, including: area under the curve (AUC) presented with 95% confidence intervals, P-value correspondence (***: <0.001; NS: not significant), the identified statistically optimal threshold, the sensitivity (with specificity fixed as 0.95, α = 0.05) for the selected threshold, and the percentage of masked pixels in the final land cover (LC) map.For the NARI mblm model, pixels with values above the optimal threshold were classified as shrublands, and the rest were classified as grasslands.Based on confusion matrices calculated using the 30% remaining polygons of the ground reference dataset (Table 2), we obtained excellent classification results (Figure 5).The Fscore was 0.90 ± 0.02 with a recall of 0.84 ± 0.03 and a precision of 0.97 ± 0.02, which denotes very low misclassification.Sensitivity and specificity were both high (0.93 ± 0.05 and 0.84 ± 0.03, respectively), indicating a low number of commission and omission errors (respectively, 0.15 ± 0.03 and 0.06 ± 0.05), which means that our classifier was missing approximatively 15% of the shrublands in our study site and was misclassifying approximatively 6% of the grasslands.The Kappa coefficient was 0.67 ± 0.05, which denotes a good classification rate according to Landis and Koch (1977).

Comparison with Existing Land Cover Products
Statistical metrics calculated between NARI-based classification and available LC products showed very high omission errors rate of shrublands (Table 3).For OSO, CCI, and CLC, the precision was 0.10 ± 0.0006, 0.08 ± 0.01, and 0.20 ± 0.15, respectively, which shows that, on average, 88% of the pixels assigned as shrubland (of adequate size) were misclassified.The omission error rates were 0.78 ± 0.001, 0.96 ± 0.006, and 0.98 ± 0.12, which means that, on average, available LC products misclassified 90% of the shrublands that they were expected to detect in mountainous areas.The Kappa coefficients were 0.09 ± 0.0009, 0.03 ± 0.009, and 0.02 ± 0.02 for OSO, CLC, and CCI, respectively, which denotes a poor classification rate.Given that our omission error rate of our classification was 15%, the estimated surface area of shrublands in our study site was approximatively 102.35 ± 2.67 km², which was approximately 10 times more than shrubland area estimated by existing LC products.Shrublands that were misclassified in available LC products were primarily confused with grasslands (76, 54, and 55% for OSO, CCI, and CLC, respectively) or sparsely vegetated area (22%) for CLC and open coniferous forest (23%) for CCI.Examples of shrubland mapping from NARI-based classification and LC products used for comparison for three test zones show clear underestimation of shrublands surfaces in widely use LC products (Figure 6).

Comparison with Existing Land Cover Products
Statistical metrics calculated between NARI-based classification and available LC products showed very high omission errors rate of shrublands (Table 3).For OSO, CCI, and CLC, the precision was 0.10 ± 0.0006, 0.08 ± 0.01, and 0.20 ± 0.15, respectively, which shows that, on average, 88% of the pixels assigned as shrubland (of adequate size) were misclassified.The omission error rates were 0.78 ± 0.001, 0.96 ± 0.006, and 0.98 ± 0.12, which means that, on average, available LC products misclassified 90% of the shrublands that they were expected to detect in mountainous areas.The Kappa coefficients were 0.09 ± 0.0009, 0.03 ± 0.009, and 0.02 ± 0.02 for OSO, CLC, and CCI, respectively, which denotes a poor classification rate.Given that our omission error rate of our classification was 15%, the estimated surface area of shrublands in our study site was approximatively 102.35 ± 2.67 km 2 , which was approximately 10 times more than shrubland area estimated by existing LC products.Shrublands that were misclassified in available LC products were primarily confused with grasslands (76, 54, and 55% for OSO, CCI, and CLC, respectively

Validation Relative to Floristic Plots
NARI slope values extracted for 30 shrublands and 70 grasslands of the floristic survey were significantly different (Figure 7; P-val < 2.2e-16).The upper panels of Figure 7 show clear examples of discrimination between shrubland (red) and grassland (green) using NARI slopes values.Omission error 0.98 ± 0.12 Kappa coefficient 0.02 ± 0.02

Validation Relative to Floristic Plots
NARI slope values extracted for 30 shrublands and 70 grasslands of the floristic survey were significantly different (Figure 7; P-val < 2.2e-16).The upper panels of Figure 7 show clear examples of discrimination between shrubland (red) and grassland (green) using NARI slopes values.

Discussion
Mapping above treeline shrublands in the Alps using multi-spectral satellite imagery represents a long standing and important challenge [9], given that shrublands represent a substantial portion of above treeline land cover, and shrub encroachment dynamics are of concern for land and conservation managers.Our study used an innovative method to differentiate Ericaceae-dominated shrublands from grasslands, thus enabling improved mapping of shrub cover in mountainous areas.We validated our approach using two independent sources of data, including photo-interpreted polygons produced using very high resolution late fall images and ground-level floristic survey data.In addition to highlighting the utility of the Sentinel-2 red-edge band as well as the importance of high spatial and temporal resolutions, we largely enhanced the detection of shrublands in our study area by demonstrating a 10-fold increase in estimated shrub cover compared to existing LC products.

Method Caveats and Perspectives
Our method relies on the sensors' capacity to detect the sharp change in leaf reflectance around 700 nm (red-edge) [59].It is a key wavelength range sensitive to vegetation conditions and has proven to be useful for quantitative assessment of vegetation properties [60], such as Leaf Area Index (LAI) [61][62][63], chlorophyll or nitrogen content [61,64], chlorophyll-a absorption in turbid and productive waters [65], mapping crop types [66], burn severity [67], and LC mapping [9,[68][69][70][71].Our work represents a novel application aimed at quantifying anthocyanin content in above treeline plant canopies.Given that reciprocal reflectance in the red-edge is only sensitive to chlorophyll and reciprocal reflectance in green is sensitive to anthocyanin, a combination of both bands allows for assessment of pigment concentration [28].Nonetheless, the formula is only viable because reflectance at 550 nm is almost equal to that at 700 nm in leaves with minute quantities of

Discussion
Mapping above treeline shrublands in the Alps using multi-spectral satellite imagery represents a long standing and important challenge [9], given that shrublands represent a substantial portion of above treeline land cover, and shrub encroachment dynamics are of concern for land and conservation managers.Our study used an innovative method to differentiate Ericaceae-dominated shrublands from grasslands, thus enabling improved mapping of shrub cover in mountainous areas.We validated our approach using two independent sources of data, including photo-interpreted polygons produced using very high resolution late fall images and ground-level floristic survey data.In addition to highlighting the utility of the Sentinel-2 red-edge band as well as the importance of high spatial and temporal resolutions, we largely enhanced the detection of shrublands in our study area by demonstrating a 10-fold increase in estimated shrub cover compared to existing LC products.

Method Caveats and Perspectives
Our method relies on the sensors' capacity to detect the sharp change in leaf reflectance around 700 nm (red-edge) [59].It is a key wavelength range sensitive to vegetation conditions and has proven to be useful for quantitative assessment of vegetation properties [60], such as Leaf Area Index (LAI) [61][62][63], chlorophyll or nitrogen content [61,64], chlorophyll-a absorption in turbid and productive waters [65], mapping crop types [66], burn severity [67], and LC mapping [9,[68][69][70][71].Our work represents a novel application aimed at quantifying anthocyanin content in above treeline plant canopies.Given that reciprocal reflectance in the red-edge is only sensitive to chlorophyll and reciprocal reflectance in green is sensitive to anthocyanin, a combination of both bands allows for assessment of pigment concentration [28].Nonetheless, the formula is only viable because reflectance at 550 nm is almost equal to that at 700 nm in leaves with minute quantities of anthocyanins [72].Thus, the capacity to detect high relative content of anthocyanin is highly dependent on the central wavelength of the band and te hbandwidth.In Gitelson et al. [28], ARI is developed using bands with red-edge central wavelength at 700 ± 7.5 nm.Because the red-edge Sentinel-2 band is wider (approximately 705 ± 15 nm), noise in the response can be expected.Band 4 (665 ± 15 nm) and band 6 (740 ± 7.5 nm) of S2 have shown no capacity to detect anthocyanin accumulation and thus replace band 5. Thus, the sensitivity of mapping results to bandwidth should be further investigated using hyperspectral data to refine shrubland detection.
Classification of shrublands as one entity suggests that each contributing species exhibits similar responses.Because various levels of reddening are observed between most contributing species (e.g., Vaccinium spp.and Rhododendron ferrugineum), and red-edge position (REP) [73] varies between species due to differences in chlorophyll concentration or biomass [74,75], different phenomena are measured by the Sentinel-2 red-edge band, and various canopy contributions can be expected.In our work, species-level contributions were not investigated due to the constraints of working at the canopy level with multi-spectral satellite imagery.Further analysis of species-level contributions to spectral signatures at the canopy level while accounting for variability over the course of the growing season will require complementary measurements using a field spectrometer [76].Furthermore, identifying target end members using a field spectrometer would enable sub-pixel fractional mapping of shrub, grasslands, and rocky habitat mosaics, which are ubiquitous features of alpine landscapes.
Among the various tools available for land cover mapping [77], multi-temporal and high-resolution satellite platforms such as Sentinel-2 have been shown to be highly effective [9] due to added information on vegetation phenology [78][79][80].While seasonal changes in NDVI have already been integrated into land cover mapping protocols [9,81], accounting for seasonal dynamics of plants pigments, such as anthocyanins, represents a novel approach for enhancing high-resolution land cover maps in mountainous areas.Here, we took advantage of the unique dynamics of anthocyanin accumulation in the late fall leaves of shrublands using a multi-temporal approach, demonstrating the potential of classification algorithms using pixel-level trends in seasonal plant pigments (chlorophylls, carotenoids, and anthocyanins) [30].Alternatively, this unique red pigment activity can be utilized for mapping using a single, very high spatial resolution image acquired at the appropriate season; however, this approach is more difficult to automate, as the timing of peak pigment accumulation varies from year to year.

Quantifying Shrub Dynamics
The encroachment of shrublands into grasslands in the Alps [17] appears to result from the interaction of multiple factors, including land-use [12,13] and climate change [82,83].At the local scale, shrub encroachment into species-rich subalpine grassland leads to reduced biodiversity [84-86] and diminished pastoral resources [19].Accordingly, the ability to track shrub dynamics over space and time and understand drivers of shrub expansion is essential to land managers.Compared to currently used LC products, our method allows for improved estimation of the spatial distribution of shrublands at a fine scale and thus moving forward provides decision-makers with a novel tool for land management [87].Using Sentinel-2, NARI-based shrubland classification can be produced on a yearly basis beginning in 2016, thus allowing for monitoring of shrub dynamics and potential encroachment into grasslands during the coming years.
Understanding the local drivers of shrub expansion, however, necessitates quantifying vegetation dynamics during recent decades.At coarser spatial resolution, MERIS and Sentinel-3 [35] can give continuous NARI time series from 2002 to present, but they are not adapted to complex topographic areas.The Landsat Time Series (LTS), which covers the period from 1984 to present, offers the opportunity to study decadal-scale changes in mountain vegetation [2].Detecting anthocyanin accumulation within larger 30 m Landsat pixels and with a reduced number of spectral bands represents a significant methodological challenge.A possibility would consist of mapping fractional shrub cover within Landsat pixels, which would require careful calibration with hyperspectral data of end member canopies and as well as calibration of spectral unmixing algorithms across different Landsat satellites.Finally, leaves that concentrate anthocyanins often display characteristics typical of leaves growing under shaded conditions, such as lower chlorophyll a:b ratios [88,89] and higher total chlorophyll content [90,91], which represents a further opportunity.
Quantifying land cover changes requires that observed changes in radiance are resulting from changes in LC and not from variations in atmospheric conditions, soil moisture, sun-surface-sensor geometry, leaf pigmentation [92], or sensor performance [93].Our method used to classify pixels as shrublands or non-shrublands is highly dependent on the spectral intensity registered at the sensor, and several non-LC related variables can influence mapping results.We are confident that using high performance algorithms [37] prevents errors related to atmospheric conditions.Nonetheless, we recommend applying the bidirectional reflectance distribution function (BRDF) correction mostly to account for "day of year" effects [94], as "view angle" and "mean local time drift" effects are not expected for Sentinel-2 at this scale.Several methods are available to correct for such effects using Sentinel-2 images [95,96].Due to the progressive accumulation of anthocyanins during an extended period (approximately two months) and the limited number of available scenes, we expect that date availability will vary from year to year.Nonetheless, we are confident that trend estimation using robust models such as mblm can partly override this issue.Finally, because anthocyanin concentration is related to climate variables [97,98], the timing and the intensity of peak accumulation of anthocyanins are expected to vary between years, leading to changes in resulting land cover maps that would be non-representative of shifts in vegetation.Although quantifying this variability remains to be carried out, our initial recommendation for long-term monitoring is to produce averaged shrub cover maps over five year periods in order to smooth inter-annual variability in climate and pigment dynamics.Finally, because temperature and irradiation vary along gradients of elevation, slope aspect, and latitude [99], substantial spatial variations in reddening intensity are expected for a given species across broad study areas.An efficient methodology for long-term monitoring of shrub dynamics at biogeographic scales will thus need to account for these considerations.

Perspectives for Shrub Mapping in the Arctic
Circum-Arctic measures of vegetation dynamics by satellites indicate widespread patterns of recent greening [100][101][102][103] associated with warming trends [104].Greening trends have been partly attributed to increases in shrub cover [13] with various responses depending on shrub species and investigated region [105,106].In the Eastern Canadian Arctic, greening trends are mostly due to increases in Betula glandulosa cover [13,107], which affect competition between species [108] and soil thermal regimes [109,110].Climate-vegetation feedbacks in the Arctic have important implications for the global climate system, given that shrub expansion into arctic tundra is hypothesized to have direct effects on rates of permafrost thaw and associated release of methane [111].In the Scandinavian Arctic, multiple artico-alpine Ericaceae shrubs are present (e.g., Vaccinium uliginosum, Vaccinium vitis-idea, and Empetrum nigrum) and could contribute to observed shrub expansion into tundra.Further testing and calibration would be required to evaluate the sensitivity of NARI-based methods for mapping additional dominant arctic shrub species, such as willows (Salix spp.), green alder (Alnus viridis), and dwarf birch (Betula nana).As in the Alps, systematic mapping of shrub cover in the Arctic is both lacking and necessary in order to upscale field-based observations of shrub dynamics to the landscape scale.We encourage further testing and validation of NARI-based approaches in an arctic context in order to determine if this method could be integrated into existing monitoring schemes carried out in the Arctic [112].

Conclusions
This work presents an innovative method for differentiating above-treeline shrublands from grasslands in a temperate alpine context based on the detection of anthocyanins in Ericaceae-dominated plant canopies during the late fall.We propose that our method could form the basis for improved monitoring of mountain shrubland dynamics in the coming years during the ongoing Sentinel-2 mission.We demonstrate the poor capacity of widely used land cover maps to identify shrublands as well as the gain in performance associated with our approach.Further work is necessary to establish a relation between anthocyanin content in leaves and NARI using field measurements of leaf spectroscopy and considering local climate conditions, topographic variables, and species contributions to canopy-level reflectance.An additional challenge consists of tracking shrub dynamics during recent decades using cross-calibration between Sentinel-2 and the Landsat archive.

Supplementary Materials:
The following are available online at http://www.mdpi.com/2072-4292/11/23/2807/s1, Figure S1: Google Earth images (Map: Google, DigitalGlobe 2016-10-04) showing the reddening phenomenon due to anthocyanins accumulation used to produces our photo-interpretated polygons of shrublands and grasslands.Figure S2: Google Earth images (Map: Google, DigitalGlobe 2016-10-04) of two typical cases of distorted signal due to surroundings environmental conditions such as in scree (A) or forest edge (B).The floristic relevé dominated by shrublands that needed to be masked are presented in red. Figure S3: Fraction of pixels covered by shrublands for CLC (A) and CCI (B).The threshold of 75% is identified as a vertical dashed line.

Figure 1 .
Figure 1.Location of the study area within the French Alps.Photo-interpreted polygons and floristic surveys plots are shown in red and green for shrublands and grasslands.

Figure 1 .
Figure 1.Location of the study area within the French Alps.Photo-interpreted polygons and floristic surveys plots are shown in red and green for shrublands and grasslands.

Figure 2 .
Figure 2. Reflectance distribution of end-members pixels representing Ericaceae-dominated shrublands (dominated by Vaccinium uligonosum) and grasslands (dominated by Nardus stricta) extracted from a hyperspectral image acquired on 28 September 2015 at the Col de Balme [1950 m at sea level (a.s.l.)] in the nearby Mont-Blanc massif.Sentinel-2A bands B3 and B5 used to calculate Normalized Anthocyanin Reflectance Index (NARI) are indicated by grey rectangles in panel (A).Photo of late fall grasslands dominated by Nardus stricta (B) and late fall Vaccinium uligonosum (C; photo credit: Cédric Dentant).

Figure 2 .
Figure 2. Reflectance distribution of end-members pixels representing Ericaceae-dominated shrublands (dominated by Vaccinium uligonosum) and grasslands (dominated by Nardus stricta) extracted from a hyperspectral image acquired on 28 September 2015 at the Col de Balme [1950 m at sea level (a.s.l.)] in the nearby Mont-Blanc massif.Sentinel-2A bands B3 and B5 used to calculate Normalized Anthocyanin Reflectance Index (NARI) are indicated by grey rectangles in panel (A).Photo of late fall grasslands dominated by Nardus stricta (B) and late fall Vaccinium uligonosum (C; photo credit: Cédric Dentant).

Figure 3 .
Figure 3. Normalized Difference Vegetation Index (NDVI) (top), NARI (middle), and masked surfaces (bottom) distribution over the three years for shrublands (red) and grasslands (green).The error bars represent the 99% confidence interval.Percentage of masked surfaces are based on the surfaces of the studied area presented in Figure 1, representing 1377 km².

Figure 3 .
Figure 3. Normalized Difference Vegetation Index (NDVI) (top), NARI (middle), and masked surfaces (bottom) distribution over the three years for shrublands (red) and grasslands (green).The error bars represent the 99% confidence interval.Percentage of masked surfaces are based on the surfaces of the studied area presented in Figure 1, representing 1377 km 2 .

Figure 4 .
Figure 4. Kernel density plots for shrublands (red) and grasslands (green) for the three NARI (left) and the three NDVI (right) models.Vertical lines represent the statistically optimal classification threshold.

Figure 4 .
Figure 4. Kernel density plots for shrublands (red) and grasslands (green) for the three NARI (left) and the three NDVI (right) models.Vertical lines represent the statistically optimal classification threshold.

Figure 5 .
Figure 5. Statistical metrics used to evaluate the performance of the NARI mblm-based model of shrubland classification.Metrics were calculated 100 times using a randomized 50% selection of the photo-interpreted data retained for validation.

Table 3 .
Confusion matrices between NARI-based classification and shrublands classes of three widely-use land cover [Occupation des Sols (OSO), Corinne Land Cover (CLC), and CCI in row] in number of pixels.Omission errors and Kappa coefficients are indicated below each corresponding table.

Figure 5 .
Figure 5. Statistical metrics used to evaluate the performance of the NARI mblm-based model of shrubland classification.Metrics were calculated 100 times using a randomized 50% selection of the photo-interpreted data retained for validation.
) or sparsely vegetated area (22%) for CLC and open coniferous forest (23%) for CCI.Examples of shrubland mapping from NARI-based classification and LC products used for comparison for three test zones show clear underestimation of shrublands surfaces in widely use LC products (Figure 6).

Figure 6 .
Figure 6.Example LC maps of shrublands for three test zones (columns), representing from top to bottom rows: continuous slope values of the NARI model, Corinne Land Cover (CLC), Global ESA CCI Land Cover (CCI), Occupation des Sols (OSO), and NARI-based classification.Shrublands are represented in red.

Figure 6 .
Figure 6.Example LC maps of shrublands for three test zones (columns), representing from top to bottom rows: continuous slope values of the NARI model, Corinne Land Cover (CLC), Global ESA CCI Land Cover (CCI), Occupation des Sols (OSO), and NARI-based classification.Shrublands are represented in red.

Figure 7 .
Figure 7. Examples of vegetation plots identified as shrubland (red) and grassland (green) juxtaposed (left) or not (right) with NARI slopes values and an aerial photograph.Pixels classified as shrubland are indicated by the red outline.Distribution of NARI slope values for shrubland (red, N=30) and grassland (green, N = 70) vegetation plots.

Figure 7 .
Figure 7. Examples of vegetation plots identified as shrubland (red) and grassland (green) juxtaposed (left) or not (right) with NARI slopes values and an aerial photograph.Pixels classified as shrubland are indicated by the red outline.Distribution of NARI slope values for shrubland (red, N=30) and grassland (green, N = 70) vegetation plots.

Table 2 .
Confusion matrix between NARI-based classification using mblm model and the remaining ground reference data of shrublands and grasslands used for model validation.

Table 2 .
Confusion matrix between NARI-based classification using mblm model and the remaining ground reference data of shrublands and grasslands used for model validation.

Table 3 .
Confusion matrices between NARI-based classification and shrublands classes of three widely-use land cover [Occupation des Sols (OSO), Corinne Land Cover (CLC), and CCI in row] in number of pixels.Omission errors and Kappa coefficients are indicated below each corresponding table.