Sentinel-2 and AISA Airborne Hyperspectral Images for Mediterranean Shrubland Mapping in Catalonia

Abstract

The Mediterranean ecosystem exhibits a particular geology and climate, which is characterized by mild, rainy winters and long, very hot summers with low precipitation; it has led to the emergence of resilient plant species. Such habitats contain a preponderance of shrubs, and collectively harbor 10% of the Earth’s species, thus containing some of the most unique shrubby formations protecting against environmental natural degradation. Due to shrub species diversity, initial phases of forestland, heterogenous grasses, bare ground and stones, the monitoring of such areas is difficult. For this reason, the aim of this paper is to assess semi-automatic classifications of the shrubby formations based on multispectral Sentinel-2 and visible and near infrared (VINR) AISA-EAGLE II hyperspectral airborne images with a support of Canopy High Model (CHM) as a three-dimensional information and field-verified patterns, based on Match-T/DSM and aerial photos. Support Vector Machine (SVM) and Random Forest (RF) classifiers have been tested on a few scenarios featuring different combinations of spectral and Minimum Noise Fraction (MNF) transformed bands and vegetation indices. Referring to the results, the average overall accuracy for the SVM and AISA images (all tested data sets) was 78.23%, and for the RF: 79.85%. In the case of Sentinel-2, the SVM classifier obtained an average value of 83.63%, while RF: 85.32%; however, in the case of the shrubland, we would like to recommend the RF classifier, because the highest mean value of F1-score achieved was 91.86% (SVM offered few-percent-point worse results), and the required training time was quicker than SVM. Commonly available Sentinel-2 data offered higher accuracies for shrubland monitoring than did the airborne VNIR data.

1. Introduction

Shrubland (alternative terms: scrubs, scrubland, shrubland, shrubby formation, chaparral) covers all kind of woody formations dominated by shrubs (perennial woody plants not exceeding 5–8 m height) with herbs, grasses and geophytes [1]. In particular regions, it may be a stable type of vegetation or a transitional stage in the ecological succession after a fire, deforestation or agricultural wastelands. The world’s main areas of shrubland occur in the Mediterranean climate (wet, mild, long winters and hot, dry summers) [2]. Natural shrublands have the lowest water consumption and better resist climate aridification compared to other stands, e.g., trees, which cannot survive such harsh conditions [3]. Shrubs facilitate rainwater infiltration, protect soil from erosion and optimize a carbon exchange, thereby limiting its emission into the atmosphere [4,5]. Disruption of this vegetation causes an increase in soil surface temperature, plant water stress and erosion effect on runoff and sediment yield [6]. The importance of the occurrence of shrubs in the Mediterranean environment is also described by Torre et al. [7], analyzing airborne lidar scanning (ALS) data, vegetation profiles for a decade period for selected study plots in Barcelona’s Natural Parks, (Catalonia, NE Spain) where the relationship of shrubs vegetation cover and the occurrence of small mammals was studied. These play a key role in the sustainability of Mediterranean forests; therefore, monitoring scrubs plots is important from a management perspective under current conditions of global change. In Catalonia, as in many Mediterranean countries, the amount of precipitation is insufficient to maintain the dense vegetation cover that helps prevent a loss of habitat biodiversity [8]. Shrublands cover approximately 14% of Catalonia, being one of the main vegetated covers [1], and monitoring of shrub spatial distribution allows biomass monitoring, which provides some of the key input data for fire risk estimation; biomass and its parameters are frequently discussed issues in the case of remote sensing shrub analysis. Calvão and Palmeirim [9] identified Mediterranean scrubs (authors use this term) located in SE Portugal on Landsat images, which allowed for the calculation of a Normalized Difference Vegetation Index (NDVI); least-squares linear regression model verified by field-acquired data from experimental sites describing leaf, woody and litter biomass, canopy cover and Leaf Area Index (LAI) allowed for the capture of phenological spectral changes of the Cistus spp. and the evergreen sclerophyllous scrubs, and correlations between NDVI, canopy cover and fresh biomass was shown to oscillate around r=−0.53, p < 0.02 [9]. In the following paper, Calvão and Palmeirim [10] also confirmed monitoring capabilities of long-term Mediterranean vegetation changes based on 12 Landsat-derived vegetation indices (VIs) characterizing the state of plants (including soil properties) and biophysical parameters (biomass and leaf area indices, canopy height and cover), which allowed for the development of regression models; normalized difference VIs correlated well with canopy-oriented variables (r = 0.83–0.86) and worse correlated with leaf-related variables (r = 0.59). The analyses confirmed that the most useful indices were Ratio Vegetation Index = Simple Ratio (RVI) and NDVI, and Renormalized Difference Vegetation Index (RDVI), Modified Simple Ratio (MSR), Leaf Area Index (LAI) and Landsat bands scored worse, non-significant correlation values (<0.58). Similar results were achieved for High-Arctic shrubs [11,12] and for alpine meadow plants [13,14,15]. Sentinel-2 images were also successfully used to identify multi-temporal shrubland areas in Inner Mongolia, China [16] and in deforested equatorial areas of Gabon [17]. In the first case, the authors [16] described the importance of the phenological factor in identifying shrubs. For this purpose, the authors used the gaussian process regression (GPR, determining similarity between samples) and random forest regression (RFR, indicating a degree of mean square error (MSE) reduction of feature variables), which were the vegetation indices presented below. The conducted analyses confirmed that early spring (early April), summer (late June-mid-July) and late October are optimal periods for shrubs analysis (in our case, various growing periods were not analyzed, as we had access to the data acquired during a single AISA flight campaign, and the Sentinel-2 image was selected in accordance with the airborne data). Gan et al. [16] also selected relevant vegetation indices for this type of analysis, which were NDTI (Normalized Difference Tillage Index), GNDVI (Green Normalized Difference Vegetation Index), NDVI, NDSVI (Normalized Difference Senescent Vegetation Index), EVI2 (two-band Enhanced Vegetation Index), MSAVI, NDPI (Normalized Difference Phenology Index), NDGI (Normalized Difference Greenness Index) and NDMI (Normalized Difference Moisture Index). As the outcome, Caragana shrubs were identified with the best overall accuracy, achieving 91.7%, and RMSE achieved 5.0–5.5% (validation–test patterns). A significant role of increasing the results were found with the NDTI and NDGI indices [16]. In the case of Gabonese research, authors tested an increase of DEM and NDVI on the accuracy of Sentinel-2 image classification carried out with the Random Forest (RF) algorithm. The results confirmed that non-forest communities achieved 76.4–96.1% (Producer Accuracy), and the DEM increased outcomes by six percent points compared to the set which contained only the NDVI index [17].

Estimated shrub biomass in formerly burned areas was estimated in many countries, e.g., in Portugal [18]; the authors used 10-m-Sentinel-2 NDVI to create a raster mask for identifying burnt areas, and using Minimum Distance Classifier and B2, B3, B4 and B8 bands (VIS–NIR), identified forested and shrubland areas with an overall accuracy (OA) oscillating around 80% (Cohen’s K coefficient k = 0.75) [18]. Bulk density for shrubs has also been studied to estimate carbon sequestration, which could assist in predicting the potential fire risk for the Iberian Peninsula [19]. Vega et al. [20], based on the field surveys and biomass measurements, raised the shrubs communities biomass modelling in Northern Spain and developed models for standing fuel loads by size and condition. Biomass estimation plays a particularly important role in determining the fuel source for shrublands fires and greenhouse gas emission. Sentinel-2, TanDEM-X and airborne LiDAR data and Random Forest with an option of automatic feature selection allowed for the mapping of shrublands with 83.1% overall accuracy (balanced accuracy scored 0.82) [21], and the achieved results allowed for the prediction of landscape biomass (Spearman rank correlation achieved 0.80–0.84, and normalized RMSE <16%). The authors [21] also highlighted an integration of vegetation type mapping towards biomass, water and carbon management modelling based on remote sensing data. More complicated spatial systems, in the form of 24 classes of homogeneous scrub communities of Western Australia were identified on multitemporal Sentinel-2 images [22]. The authors tested a few scenarios oriented on growing seasons, data compression methods and used input data, and the best results were achieved for Sentinel-2 spring and autumn spectral bands and tasseled cap transformed data classifying by the SVM (average value oscillated between 72–74%), worse results were noted for Nearest Neighbour (NN)–72%, RF–65% and Classification Trees (CT)–50% [22]. Similar research was conducted in Portugal in frame of a M.Sc. thesis [23]; Sequeira and his supervisors used multitemporal Sentinel-2 data, which classified 31 classes: 14 crop and 17 land cover types by RF with an overall accuracy of 76%, and the shrubland scored 70% (user accuracy, UA), 63% (producer accuracy, PA) and achieved an F1-score of 67%. The researcher also confirmed that shrubland, grasslands, bare soil and sparse vegetation were often mixed, because spectral signals (mixels) consisted of green, dry vegetation and different fractions of soil and shadow [23].

Different classifiers were used for land cover and plant community identification; based on 144 bands of hyperspectral images (in the spectral range of 380–1050 nm), Khodadadzadeh et al. [24] compared different classifiers, e.g., SVM and CNN (Convolutional Neural Network) for vegetation mapping [25,26]. The best results were achieved for: Joint Sparse Representation (86.86% of overall accuracy), while CNN and SVM scored 86% [26]. Often SVM, RF and other machine learning methods are successfully implemented to identify plant communities, as well as species achieving comparable level of overall accuracies (85–92%) [27,28,29,30]. In the papers, the authors used data spectral bands with different remote sensing indices, e.g., Enhanced Vegetation Index (EVI) or NDVI [31,32] as inputs, which made it easier to distinguish types of communities [33,34] and species [35]. In addition, Airborne Laser Scanner (ALS) data were also used, which allowed for the characterization of the multi-story spread of land cover and for the classification of the different species [36]. Integration of PlanetScope images and LiDAR point clouds gave good results for the land cover classification [37]. The authors, based on random trees method, mapped transitional shrubs with herbaceous communities areas with an accuracy of 97.19% (UA) and 87.82% (PA) [37]. The Random Forest classifier offers faster data processing time than the SVM classifier while offering comparable average accuracies to the SVMs [38]. Application of Sentinel-2 data, specifically the red-edge band for detection and mapping of mountain shrublands, was also highlighted by Bayle et al. [39], who compared the effectiveness of Normalized Anthocyanin Reflectance Index (NARI, which is a normalization of the ARI = 1/R550 – 1/R700 index [40]) to NDVI for shrubby formation mapping. NARI, indicating anthocyanin content in the canopy, achieved 0.92, and NDVI 0.58 (measuring the Area Under the Curve, AUC), which, after verification of the model, allowed for the evaluation of the final shrub cover map with a Kappa coefficient equal to 0.67, confirmed that the Sentinel-2 red-edge band delivers extended opportunities to detect shrubby formations.

Due to climate changes manifested by an increase in the frequency of high air temperatures and a decrease in the amount of rainfall, the Mediterranean ecosystem is exposed to water stress, which leads to a risk of vegetation, soil and rhizosphere destruction. Hence, it is necessary to carry out a detailed inventory of the current state and to monitor the intensity and directions of changes taking place. The research problem remains the optimization of tools and data that allow for constant monitoring of the environment covering the entire Mediterranean basin. Due to the amount of the available satellite data, it is necessary to assess individual input data to determine whether it is better to use all spectral bands or to reduce them using well-known MNF or PCA methods, additionally adding dedicated remote sensing vegetation indices that strengthen the information potential of analysed shrubs, which create heterogenous patches with surrounding grasses. Of course, there are also digital surface models, e.g., generated from LiDAR data or high-resolution aerial photos, but these data are not widely available for all countries. In addition, there is uncertainty about whether the vertical plant structure needs to be analysed, or if the shrubs are spectrally different enough to identify them without additional input data. Another motivation for doing this type of research is that regeneration processes in such unfavorable conditions are also difficult because the disturbed plant canopy raises soil temperature, thereby increasing evapotranspiration, and the surrounding vegetation is exposed to water stress. This leads to changes in the structure of habitats, opening up a possibility for the emergence of alien invasive or native expansive species, which could potentially destroy biodiversity. Multispectral satellite images are a natural data source for the monitoring of such vast areas, they also ensure regular repeatability according to the same conditions. However, due to the very high heterogeneity of the land cover (bushes, grasses, exposed soil and stones), acquired images represent mixels that depend mostly on a content of green pigments and water. Therefore, airborne high-resolution data, particularly hyperspectral data, offer an optimal set of information, but they are limited to local studies (also because of acquisition costs). Nevertheless, they provide valuable reference data for analyses based on satellite images and the verification of developed algorithms. For this reason, the papers aim to assess the information potential of free Sentinel-2 satellite and airborne AISA-EAGLE II hyperspectral images with the support of the Canopy Height Model (CHM) to identify heterogenous shrubland areas. The aim was to select the most accurate method and data set receiving the highest shrubland identification results for both Sentinel-2 and hyperspectral AISA data. It is extremely important to monitor these areas, especially with such an important phenomenon associated with scrubland areas, to determine the fuel loads of the scrub extensions. Many government institutions around the world, including the co-authors of the study, work on optimizing Earth Observation data and field transects in order to adjust models estimating the biovolume, biomass and biophysical parameters of shrublands, in particular, as an input for ignition and propagation models, and for the calculation of fire risk.

2. Materials and Methods

2.1. Study Area

The research areas are located in two types of biogeographical places in Catalonia (southern Spain), labelled as A (an Atlantic vegetation type) and B (a Mediterranean vegetation; Figure 1; Figure 2). The A area (Valarties valley) is located in the Aran Valley (mountain Pyrenean region), which creates a special topo-climate due to its orography and Atlantic dominancy, and has rich and varied vegetation (forests, shrubby formations, sub-, and alpine meadows). The B area covers the mountain Montsec region (the pre-Pyrenees), which has a Mediterranean climate. The calcareous mountains cover 186.96 km² (40 km long) and are catalogued as a natural park of the Natura 2000 Network: Special Protection Area (SPA) allowing coexistence with diversified flora (Serra del Montsec Protected Natural Area) [41]. The shrubby formations are formed mainly by shrubs (perennial woody plants not exceeding 5–8 m in height), and also include grasses, herbs, and geophytes, which can be natural, or sometimes anthropogenic (recovered vegetation after a fire). Natural shrubby formations occur in areas where factors, such as climate or soil, create conditions in which the growth of trees is prevented, but the growth of shrubs is possible.

Shrubby formations are spread across the entire Catalan territory, from littoral to alpine areas, some of them labelled as habitats of communitary interest in the EU. In Corine Land Cover level 2, this formation is described as scrub and/or herbaceous vegetation associations (code 3.2; shrubby formation with single trees composed of dwarf forms of birch, willow and other species). However, on level 3 of the Corine land cover map–Transitional woodland-shrub was selected (code 324), and included different-size and condition trees, shrubs; herbaceous vegetation and natural bare surfaces (Transitional woodland/shrub) [42]. Under this physiognomic type are different botanical communities, and biodiversity is dynamically changing due to the spread to abandoned agricultural areas on the one hand, and the existence of fires on the other hand. Monitoring of such ecosystems is necessary in order calculate biovolumes to infer biomass, or CO₂, as shrubs are involved in combustible models (calculating fire risk) and the modelling of post-fire shrub restoration in combination with soil analysis [43]. In the case of soil modelling, the second derivate transformation allows, by post-processing manipulation, for the identification of proper spectral properties [44,45].

2.2. Input data

Sentinel-2 and AISA airborne hyperspectral images that were atmospherically and geometrically corrected were used as input data. An AISA scanner was configured to acquire two-meter resolution for 256 bands (spectral range 394.04–997.43 nm,

Table 1. Input data used in the studies.

Sentinel-2	bands	B2	B3	B4	B5	B6	B7	B8	B8A	B11	B12
	central wavelength (nm)	490	560	665	705	740	783	842	865	1610	2190
AISA	bands	B1 to B256
	spectral range (nm)	394.04–997.43
Additional data	Land Cover Map of Catalonia [46]
	Digital Mapping Cameras DMC-3 Orthophotoimage
	Topographical Map of Catalonia 1:5000 [47]
	Habitats Map of Catalonia 1:50,000 [48]
	Canopy High Model of Catalonia (2018, 2019)

). Two AISA flight lines over NW Catalonia (Atlantic Pyrenees and pre-Pyrenees) were acquired on 24 July of 2019 and Sentinel-2 data (10 bands, Table 1) was captured the following day (25 July 2019; all images were cloud-free and good quality).

As additional data, the Land Cover Map of Catalonia, Topographical Map of Catalonia, Habitats Map of Catalonia and Digital Mapping Cameras DMC-3 Orthophotoimage was used as support material for correcting the training polygons that had been prepared for classification. It was also helpful for the visual verification of the classification results. The recent version of the Land Cover Map of Catalonia [46] was used, formerly El Mapa de Cobertes del Sòl de Catalunya MCSC (2009–2018), updated from 241 classes to 41. Topographical Map of Catalonia 1:5000 [47], accuracy: 1 m in planimetry and 1.5 m in altimetry. Habitats Map of Catalonia 1:50 000 [48] has extracted shrubland and grassland covers. Digital Mapping Cameras DMC-3 Orthophotoimage captured simultaneously with AISA images, 4 bands (infrared included), 45 cm of GSD. The Cartographic and Geological Institute of Catalonia generates, by correlation, a complete model for all of Catalonia, which is a high-resolution (1m) raster dataset. The DSM is created at 1-m resolution from aerial imagery using automatic image correlation methods. The DTM was derived from ICGC LiDAR (Light Detection and Ranging) coverage flights. Match-T/DSM automatically generated dense point clouds and digital terrain and surface models (DTM/DSM) from aerial image blocks. CHM was finally obtained by subtracting the ground surface (DTM) from the visible surface (DSM).

All of the above-mentioned materials allowed us to prepare training and verification (testing the classification results) polygons (45 polygons for each class), which were verified during a field reconnaissance (

Table 2. Number of field-verified pixels for the polygons of the areas A and B.

Class	Number of Pixels–Sentinel-2		Number of Pixels–AISA		Total Polygons
	A	B	A	B	A	B
Grasslands	1750	1983	44,654	49,548	80	166
Shrubby formations	2579	3581	63,807	89,807	120	140
Soil	850	668	21,141	12,618	52	168
Deciduous forests	2120	1909	52,750	48,165	41	142
Evergreen forests	3368	1720	83,912	42,743	52	110
Total	10,667	9861	266,264	242,881	345	726

, Figure 3). The proportions between classes for each area vary, indicating an unbalanced dataset; therefore, each data set was balanced with a threshold of 200 pixels. The main reason for this is to give equal priority to each class during the training process [49,50].

The surface reflectivity is defined as the fraction of incident electromagnetic power that is reflected by the ground towards an observer and it is a function of radiation wavelength, but also of both illumination and observation geometry. The information contained in the radiometrically corrected images provided by the AISA Eagle II is essentially radiance, which is an indirect descriptor of surface properties, because it is the result of multiple interactions among the sun radiation, the atmosphere and the surface itself. The atmospheric correction is the part of the processing chain that compensates for the illumination and observation geometry dependences, and finally, provides the value of surface reflectivity for each hyperspectral band at pixel level. Briefly, for each pixel of the imaged scene, it considers both geometric and atmospheric information. The geometric part is essentially related to each pixel’s absolute coordinate and the position of the sun and the sensor in terms of azimuth and elevation when it was imaged. The atmospheric process considers the properties of the medium of propagation in terms of gas concentration, aerosols, column water vapor (CWV), etc., and is usually handled using an atmosphere simulator. For the correction of the AISA Eagle II data, version 5.0 of the MODerate resolution atmospheric TRANsmission (MODTRAN) software was used. In the case of the Sentinel-2 imagery, the atmospheric correction was performed by means of the Sen2Cor v2.8 [51] processor, which is freely distributed by the European Space Agency (ESA). Sen2Cor is a Level-2A processor which corrects single-date Sentinel-2 Level-1C Top-Of-Atmosphere (TOA) products from the effects of the atmosphere to a Level-2A Bottom-Of-Atmosphere (BOA) reflectance product. During the Sentinel-2 images correction, a 10 m × 10 m Digital Terrain Model (DTM; derived from aerial photogrammetry by the ICGC) was used in order to obtain a precise terrain correction, whereas for AISA, a DTM with 1 m resolution was used.

The research target was heterogeneous shrubby formation, which diversified spectrally due to their size, age and species, therefore, Canopy High Model (CHM) and remote sensing indices of vegetation were applied. Indices characterized canopy and condition of the vegetation (EVI) [31], chlorophyll concentration in vegetation (GNDVI) [52], relative abundance of chlorophyll (MCARI1) [53], cover of the vegetation with eliminated soil-induced variations in vegetation indices, allowing for the analysis of soil variation [54] and the measurement of the light absorption efficiency for photosynthesis (SIPI3, plant health and stress analysis) [55] (

Table 3. Remote sensing vegetation indices (VIs) calculated on Sentinel-2 and AISA data.

Acronym	Name	Formula	Reference
EVI	Enhanced Vegetation Index	2.5 × ((R800 − R670)/(R800 − (6 × R670) – (7.5 × R475) + 1))	[31]
GNDVI	Green Normalized Difference Vegetation Index	(R800 − R550)/(R800 + R550)	[52]
MCARI1	Modified Chlorophyll Absorption Ratio Index 1	1.2 × (2.5 × (R790 − R670) − 1.3 × (R790 − R550))	[53]
OSAVI	Optimized Soil Adjusted Vegetation Index	(1 + 0.16) × (R800 − R670)/(R800 + R670 + 0.16)	[54]
SIPI3	Structure Intensive Pigment Index 3	(R800 – R445)/(R800 − R650)	[55]

). Indices were selected on the basis of the literature, which indicates a wide spectrum of applications of the above-mentioned indices in the forest research [56,57,58] and shrubs [10,34] with the use of classifiers used in this study.

Due to additional sources of information, a preselection of spectral bands, especially in the case of hyperspectral images, was needed, therefore, Principal Component Analysis (PCA) and Minimum Noise Fraction (MNF) transformations were applied to reduce the number of spectral bands [59,60], which allowed for the addition of additional indices. Each of these dimensional reduction techniques works on a unique principle and has its own advantages and disadvantages. Researchers use single transformation techniques (PCA, MNF, ICA) or hybrid combinations of these techniques (PCA and ICA transforms performed on the MNF components) before classification, then compare the results [60], emphasizing that the transformation e.g., MNF allowed for the extraction of the most informative data, with high signal-to-noise ratio (SNR), especially in the first MNF bands [60,61]. In the case of PCA, bands are delivered as optimal data due to its maximized total variance from the original images [62]. According to researchers, MNF and PCA bands have great impact for an improvement of classification results [61,63] because PCA techniques work on the data variance, and MNF sorts the information based on the signal-to-noise ratio [60,64,65,66]. In the case of the AISA data, the first 15 MNF bands were highly decorrelated and used for classification, and in the case of MNF, there were the first four bands from the Sentinel-2 images (

Table 4. The 28 scenarios were tested simultaneously on Sentinel-2 and AISA processed data and supporting materials like CHM, VIs or both of them; + data sets for the Sentinel-2, + data set for the AISA, with information how many bands was used for all data, MNF, PCA and the abbreviation name for the used VIs.

	S-2	+ CHM		+ Indices		+ CHM+ Indices		AISA
all bands	10 bands +		+		+		+	256 bands +
		+		+		+
MNF	4 bands +		+		+		+	15 bands +
		+		+		+
PCA	3 bands +		+		+		+	4 bands +
		+		+		+
Indices	EVI, GNDVI, MCARI1, OSAVI, SIPI3		+					EVI, GNDVI, MCARI1, OSAVI, SIPI3
		+

) and, respectively, four AISA’s PCA and three Sentinel-2’s PCA bands (Table 4).

For the classification, a combination of all possible scenarios was tested to use both all spectral bands and remote sensing vegetation indices (in these scenarios, a significant improvement in the classification results was not expected, as the information contained in the indices was somewhat duplicated from the spectral bands). On the other hand, it was expected that the use of plant indices and transformed information contained in the MNF and PCA bands should not only accelerate the classification time due to the reduced amount of the input files, but also improve the classification accuracies, as data of different qualitative data were implemented in the procedure. An interesting element for testing was the effectiveness of the CHM; would those data allow for increased identification of the shrubs (do the bands contain sufficient spectral characteristics to classify this heterogenous formation with satisfactory accuracy?) Thus, for all the hyperspectral and multispectral images: all the bands, selected MNF or PCA bands and vegetation indices (VIs), were classified separately, each as a separate dataset. Then, each of these data sets were fused with CHM and calculated vegetation indices, or both of them, thereby creating a fused data set using Layer Stacking method. The following 28 scenarios were tested simultaneously on Sentinel-2 and AISA processed data and the supporting materials, CHM and VIs (Table 4).

2.3. The Classification Algorithm and Accuracy Assessment

RF and SVM classifiers were used in the R Studio 4.0.0; for this purpose, the raster [67], rgdal [68], randomForest [69] and e1071 packages [70] were implemented. The RF classifier is popular due to its low parameterization requirements and good classification results [71], offering better results than SVM in case of a lower amount of training pixels, e.g., 50–100 [72], noisy data or outliers [73,74,75]. This can be caused by multiple decision trees and randomized subsets [76]. In our case, signature pixel values representing different classes were selected from field-verified training polygons and randomized sampled (50:50) as training/validation patterns, limiting the phenomenon of overfitting during classification [77]. Each classification iteration used 500 trees (ntree parameter), which comes from the optimization process to obtain repetitive results [78,79]. The SVM, as a non-parametric classifier, allowed for a flexibility of training data, thus limiting errors and leading to high accuracy [28,80]. A radial basis function (RBF), due to its efficiency, was chosen [80,81] and also was recommended for Corine Land Cover classification (including 324 class–transitional woodland/shrub) [79]. The learning parameters were confirmed by the tuning process (gamma value of 0.1 and cost of 1000).

In this study, an iterative classification and accuracy assessment were applied (100 times repeated on randomized selected patterns). The results were assed based on the overall accuracy (OA), the Kappa coefficient [82], user accuracy (UA) and producer accuracy (PA) [83]. Based on the UA and PA, the values were calculated as an F1-score [84]. The final classification map was selected for each data set (AISA and Sentinel-2) based on the highest average F1-score results. For the visualization of the final maps, a 4x4-pixel filtration window was applied. This was to remove a large number of very small polygons, whose presence on small-scale maps disturbs the perception of the results (published maps), while the database contains all classified pixels and can be properly visible when zooming in on a given fragment of the map.

3. Results

The highest results of the classification of the shrubland were achieved for the Sentinel-2 images and the RF classifier (88.46% of the average value of the F1-score for all classes and 91.86% of the F1-score for the scrubs area; the set of multispectral bands with CHM; Figure 4,

Table 5. Average values of: Overall Accuracy (OA), Kappa coefficient and F1-score for AISA and Sentinel-2 data sets.

Classifier	Data Set	AISA Data Set			Sentinel 2 Data Set
		OA [%]	Kappa	F1 [%]	OA [%]	Kappa	F1 [%]
RF	Image	79.25	0.73	78.62	87.60	0.83	85.61
	MNF	81.30	0.75	80.30	85.02	0.80	82.36
	PCA	76.95	0.69	75.19	83.05	0.78	81.61
	INDICES	76.10	0.68	75.42	78.34	0.71	77.18
	INDICES+CHM	78.62	0.72	79.80	81.18	0.75	82.11
	Image+INDICES	79.24	0.73	78.91	89.76	0.86	88.46
	Image+CHM	79.87	0.74	79.04	90.11	0.87	87.95
	Image+INDICES+CHM	80.37	0.74	79.80	83.49	0.78	82.11
	MNF+CHM	81.35	0.75	79.96	85.46	0.81	83.09
	MNF+INDICES	82.96	0.77	81.64	85.18	0.80	83.56
	MNF+CHM+INDICES	84.49	0.79	82.15	88.05	0.84	86.62
	PCA+CHM	78.62	0.72	76.89	85.35	0.81	83.21
	PCA+INDICES	78.45	0.71	77.37	83.94	0.79	82.66
	PCA+CHM+INDICES	80.34	0.74	79.08	87.90	0.84	87.03
SVM	Image	78.80	0.72	78.06	85.84	0.81	83.61
	MNF	76.06	0.68	75.42	82.84	0.77	79.30
	PCA	77.97	0.71	76.64	85.00	0.80	82.83
	INDICES	77.55	0.70	76.68	79.62	0.73	78.27
	INDICES+CHM	79.29	0.72	76.13	81.56	0.76	80.18
	Image+INDICES	77.62	0.71	76.35	85.72	0.81	84.01
	Image+CHM	79.93	0.74	79.03	86.22	0.82	83.99
	Image+INDICES+CHM	77.54	0.71	76.13	82.72	0.77	80.18
	MNF+CHM	77.48	0.70	76.28	82.81	0.77	81.02
	MNF+INDICES	76.77	0.69	75.41	82.47	0.77	79.96
	MNF+CHM+INDICES	79.96	0.73	77.54	84.04	0.79	80.56
	PCA+CHM	78.80	0.72	76.58	84.38	0.79	81.02
	PCA+INDICES	78.41	0.71	77.12	83.03	0.77	80.56
	PCA+CHM+INDICES	79.01	0.72	78.08	84.58	0.80	82.78

). In the case of AISA imagery, the highest achieved mean value (82.15%) of F1-score was obtained for the all classes and 82.52% of F1-score for the scrubs area for the MNF with the vegetation indices combined with the canopy height model (Random Forest classifier, Table 5).

In addition, it is relevant to note that the five best RF classification results of the shrubland apply to multispectral data, and have vegetation indices, MNF and CHM data sets. This is related to the diversity of shrubby formation, especially in combination with the 10-m pixel and the 2-m pixel for hyperspectral data. The higher classification accuracy values were obtained for the RF classifier, which means that with similar results, the computation process could be quicker than other options. Additional accuracy superiority for Sentinel-2 data over AISA is related to the range of AISA that was used for classification at the wavelength of 394.04–997.43 nm.

The aim of the research was to assess the potential of Sentinel-2 and AISA images, as well as RF and SVM classifiers, for shrubland mapping (analyzing all classes; Table 5); however, looking at the overall accuracy and Kappa coefficient results, the main focus is on the effectiveness of the analyzed classifier, which allows for the identification of all classes. Thus, it can be noticed that the Sentinel-2 data allowed for the obtainment of an average accuracy of 84.47% (regardless of whether RF or SVM was applied), and for AISA data, the accuracy was 79.04% (AISA images covered only the VNIR range). Looking at the results of classifiers, RF was slightly more effective (79.85% for AISA data, while for SVM, it was 78.23%), and Sentinel-2 data achieved 85.31% and 83.63%, respectively. In both analyzed hyperspectral and multispectral data sets, the highest results were obtained for the RF classifier. The average overall accuracy for the RF classifier for AISA data (all tested data sets) was 79.85%, with an average Kappa value of 0.73. The SVM classifier–78.23% for overall accuracy and 0.71 Kappa values, respectively. For the data sets for Sentinel-2, the RF classifier obtained an average overall accuracy of 85.32% (mean Kappa value for all sets 0.80), while the SVM classifier had an average overall accuracy of 83.63% (mean Kappa value for all sets 0.78; Table 5).

Shrubby classification results are differentiated for both areas A and B; the highest four data sets of AISA data contained MNF+CHM+INDICES; MNF+INDICES; MNF+CHM; MNF. Nevertheless, the best data set was selected – MNF+CHM+INDICES (based on F1-score and OA and Kappa values). The overall accuracy for area A was 90.44%, and for area B–78.55%, while for Kappa, the accuracies were 0.87 and 0.71, respectively. This differentiation is due to the variable surface forms of the area and their location in both areas. It is also worth emphasizing that shrubs often mix with grasslands or forests during classification, depending on their density. This caused a decrease of accuracies for the grassland, shrubby formations and deciduous forest classes (Figure 5).

In the case of the shrubby identification on AISA hyperspectral data (MNF+CHM+Indices), it was noted that the results were closer in area A–82.83% than in B–82.20% F1-score (Figure 5). The class was mistaken for grasslands in both areas, and also in A area for deciduous and coniferous forest; while in B area, it was also mistaken for coniferous forest (similar spectral properties of the coverage of both classes;

Table 6. Confusion matrix for A and B area for each class (RF and AISA_MNF+CHM+INDICES data set). Explanation: D.f.–Deciduous forest, C.f.–Coniferous forest.

A		Reference Data
		D.f.	Grasslands	C.f.	Shrubs	Soil
Classified data	D.f.	23,401	0	944	1281	8
	Grasslands	1	18,192	1	2200	235
	C.f.	1646	1	39,750	1078	15
	Shrubs	141	1719	1112	20,199	288
	Soil	0	266	0	553	7112
B		Reference Data
		D.f.	Grasslands	C.f.	Shrubs	Soil
Classified data	D.f.	14,737	27	1202	1243	0
	Grasslands	7	14,820	1	2845	213
	C.f.	2746	14	10,511	2615	8
	Shrubs	479	2908	1638	28,610	81
	Soil	0	2810	0	582	2441

). In the case of Sentinel-2 data, a set of Image + CHM offered better results. The overall accuracy for the A area oscillated around 95.21%, and for B–85.01%, while for Kappa, the accuracies were 0.94 and 0.80, respectively. (Figure 6;

Table 7. Confusion matrix for A and B area for each class–the RF classifier for Image+CHM data set for Sentinel-2. Explanation: D.f.–Deciduous forest, C.f.–Coniferous forest.

A		Reference Data
		D.f.	Grasslands	C.f.	Shrubs	Soil
Classified data	D.f.	657	0	24	34	0
	Grasslands	4	696	0	49	6
	C.f.	6	0	1591	4	0
	Shrubs	7	39	36	1021	2
	Soil	0	5	0	1	348
B		Reference Data
		D.f.	Grasslands	C.f.	Shrubs	Soil
Classified data	D.f.	581	0	83	5	0
	Grasslands	0	639	0	45	45
	C.f.	91	1	580	59	0
	Shrubs	2	99	32	1372	7
	Soil	1	109	0	6	145

).

Sentinel-2 images created with additional features (Image+CHM; Image+INDICES; MNF+CHM+INDICES; PCA+CHM+INDICES) optimal sets for shrubby formation mapping; but the highest scores were achieved for Sentinel-2 Image and CHM. Some shrub areas were misclassified in area A as deciduous forest and grasslands, and in area B, more as coniferous forest than grasslands (Figure 6, Table 7). This also confirms the same relationship in the matrices using the example of the classification results on AISA data discussed earlier.

The classification of area A and B on Sentinel-2 data and AISA hyperspectral data gave a satisfactory result, and it was presented in the form of a map with approximated terrain segments for better visualization of the result (Figure 7 and Figure 8). Both areas were clipped to the same range for AISA data. Figure 7 shows the result for area A for both the Sentinel-2 data (Image + CHM data set) on the left side of the figure and the AISA data on the right side of the figure (MNF+CHM+INDICES data set). In addition, the approximate sections of the classified area for these two data cases are shown in the center. Figure 8 shows the same set of results for area B.

4. Discussion

Thus far, there have been not many studies on the identification of shrubby formation for the Mediterranean zone focusing on shrubby species. Due to the Corine Land Cover program, class 324 of transitional woodland/shrub is well-analyzed (Table 6). As an example, the outcomes of our team [79] highlighted the superiority of SVM RBF over the RF classifier, up to five percentage points, and Sentinel-2 images allowed for the identification of Catalonian transitional woodland/shrub (class 324) and sclerophyllous vegetation (323) with an average F1-score oscillating between 73% and 88% (F1-score achieved 80% for the 60 m resolution). In the case of the SVM RBF algorithm and multi-temporal Sentinel-2 images, the class 324 achieved OA: 76%, Kappa coefficient 0.74 [79]. In the present work, the average for both A and B areas’ F1-scores was ten percentage points higher because we focused on smaller, homogenous patterns of scrubs with original Sentinel-2 pixel sizes.

Our result of the classification of the Catalan shrubs on the Sentinel-2 images reached the highest values on all spectral bands–85.61% (F1-score, RF). Adding the vegetation indices, the accuracy increased to 88.46%, and for the image and CHM: 87.95%. Similar values were achieved based on MNF and PCA combinations. The SVM classifier achieved a few percentage points lower, as the highest was 84.01%. In the case of Portuguese shrubland, Aranha et al. [18] used the Minimum Distance classifier and Sentinel-2 data with 10 m resolution, and obtained very comparable results: 92% user accuracy and producer accuracy–84%, which corresponds to the F1-score value of 87.82%. The authors additionally recommended Sentinel-2’s B2, B3, B4 and B8 for the classification. In the conclusion, it was also emphasized that the studied biomass increased the fire hazard after the following five years [18]. As presented above, our results are better than those based on the similar vegetation indices, PCA and Tasseled Cap bands [22]; the authors mapped 24 Australian classes of scrubs vegetation based on Sentinel-2 data, and the SVM with the best accuracy was 72–74% [22].

The proposed solutions confirm that the obtained accuracy values are competitive in relation to the references (

Table 8. Comparison of the obtained results with achievements of other researchers.

Author	Data Used (Data Set)	Object of Classification	Classifier	Measures of Accuracy
Our results	AISA; Sentinel-2	Shrubby formation	RF	F1: 73.70–84.80%; UA: 81.10–89.90%; PA: 75.20–79.90%
			SVM	F1: 73.10–78.20%; UA: 77.20–81.60%; PA: 75.80–79.30%
			RF	F1: 84.14–91.67%; UA: 86.01–97.20%; PA: 88.20–93.07%
			SVM	F1: 81.50–88.90%; UA: 81.40–89.20%; PA: 87.40–91.20%
[85]	UAV data, ALS,	Mediterranean Shrub Species	Decision Tree Extra Trees Gradient Boosting Random Forest MultiLayer Perceptron	Mean Cross-Validation Score: 78–80% 81–86% 81–87% 81–83 % 78–81%
[86]	Sentinel-2 Sentinel-1	Shrubland	RF	UA: 75.01%; PA: 66.8%
[87]	UAV RGB (original bands, original bands + spectral indices, and original bands + spectral indices + texture indices) + PCA	Shrubby species	RF SVM	OA: 78%; Kappa: 0.65 OA: 76%; Kappa: 0.63
[79]	Sentinel-2	Corine Land Cover; 324: transitional woodland/shrub	SVM RF	Romania: F1: 69–87%; UA: 80–91%; PA: 56–94% Catalonia (Tarragona): F1: 89%; UA: 97%; PA: 82% Poland: F1: 87%; UA: 80%; PA-94% Romania: F1: 32–70% Catalonia (Tarragona): F1: 81% Poland: F1: 41%
[88]	Sentinel-2	Mediterranean shrub	RF SVM QDA	OA: 81.9%; Kappa: 0.80 OA: 46%; Kappa: 0.40 OA: 68%; Kappa: 0.60
[89]	AISA-FENIX	Shrubland	SVM	OA: 86.1%
[90]	S-2, RGB	Shrub land (Shrub)	DeepLabV3+ SVM (GLCM) RF (GLCM)	OA: 71.57% OA: 76.70% OA: 73.64%
[91]	Sentinel-2 Landsat 8	Corine Land Cover; 324: Transitional woodland scrub	RF	UA: 88.7%; PA: 77.60 % UA: 62%; PA: 46.50 %
[23]	Sentinel-2	Shrubland	RF	F1: 67%; UA: 70%; PA: 63%
[92]	Sentinel-2	Shrubland	RF	F1: 67%; UA: 70%; PA: 63%
[18]	Sentinel-2	Shrubs	Minimum Distance Classifier	UA: 92%; PA: 84%
[21]	Sentinel-2 + TanDEM-X + Lidar	Shrublands	RF	OA: 83% (+/−15%); Kappa: 0.81; Shrublands Balanced Accuracy: 81.9%
[22]	Sentinel-2	Shrublands	SVM NN RF CT	OA: 68 % OA: 62 % OA: 63 % OA: 43 %
[93]	Sentinel-2	Shrubland	RF	UA: 79.70%; PA: 79.70%; OA: 89.20%; Kappa: 0.87
[94]	Bi-temporal ZiYuan-3 multispectral and stereo images	Scrubs	SVM	UA: 84.4%; PA: 71.7%
[95]	Landsat OLI Landsat OLI + Lidar	Shrublands	RF	UA: 51%; PA: 86%; OA: 47% UA: 90%; PA: 85%; OA: 77%
[96]	Sentinel-2, Sentinel-2 + LIDAR	Shrub	RF	UA: 59%; PA: 81% UA: 68%; PA: 84%
[97]	Sentinel-2 Sentinel-1	Shrub land	RF	UA: 79.04%; PA: 80.22% UA: 87.21%; PA: 86.58%
[98]	Sentinel-2 + DEM	Sparsely vegetated areas (includingorchards and shrubs)	ML SVM	UA: 93.96%; PA: 21.24% UA: 57.88%; PA: 45.05%

).

The results of the classification also depend on the vegetation period in which the airborne data were acquired. It was confirmed that the optimal image acquisition period for plant communities mapping is from the second half of August to the first half of September in Central European vegetation periods, due to the complex physiognomy [13,14,15,78]. This was also confirmed by Calvão and Palmeirim [9], who observed a summer stress (an intense radiation and a water deficit), which changed the registered spectra of different types of vegetation included in the shrubby formations. It was also noted, as by previous researchers Gan et al. [16], that shrub land cover classification is not an easy task using remote sensing tools. In general, the shrub canopy is heterogeneous due to combinations with other targets (grasses, soil, single trees and shades), and an additional difficulty is the problematic separation of shrubs from grasslands, as they show similar spectral characteristics.

The optimization of tools and data allows for the use of Sentinel-2 data for constant monitoring of the shrubland cover, and also ensures regular repeatability according to the same conditions. This is extremely important for determining the fuel loads of the scrub extensions, models estimating biovolume, biomass and biophysical parameters of shrublands, and for the calculation of fire risk.

The proposed method, which is based on PCA/MNF bands, remote sensing indices, Canopy High Model data and Random Forest classifier, offered 5–20% higher accuracy than those obtained so far in the literature. It is important to highlight Sentinel-2’s dense time series, which allows it to compensate for its lower spectral and spatial resolution compared to the expensive, but much more accurate, airborne hyperspectral AISA data. Due to the high solar radiation and air temperatures, Mediterranean plants have a number of protective mechanisms depending on species, e.g., they are covered by waxes, or their leaves are smaller and change their position due to the sun’s radiation, which make their identification difficult because the woody parts and the soil generate mixel signals depend on shrub species. On the other hand, the shrubby formation provides a biodiversity and protects soil and living organisms from excessive overheating.

5. Conclusions

The present study confirmed a sufficient potential of Sentinel-2 and airborne AISA data to monitor the sensitive and heterogeneous ecosystem of the Mediterranean shrubland. Sentinel-2 images offer better results than AISA, but it should be remembered that the scope of hyperspectral data in this case concerned only visible and near infrared (VNIR) ranges. Outcomes of the Random Forest classifier indicated that differences between average values of the overall accuracies for all data sets amounted to a few percentage points (6–7%). Our research confirmed that remote sensing is the right tool to conduct permanent and objective monitoring, because:

• The shrubby formation was identified with an average overall accuracy above 91.86% based on the Sentinel-2 data; the best data set contained all bands + CHM and Random Forest as the classifier.
• AISA images offered 84.49% overall accuracy, the best data set: MNF + the vegetation indices (EVI, GNDVI, MCARI1, OSAVI, SIPI3) combined with the canopy height model and Random Forest classifier.
• Random Forest offered slightly better results than Support Vector Machine classifier; average overall accuracy for the SVM and AISA images (all tested data sets) was 78.23%, and for the RF: 79.85%. In the case of Sentinel-2, the SVM classifier obtained an average overall accuracy of 83.63%, while the RF classifier: 85.32%.
• We recommend the use of indices for characterizing the canopy and condition of the vegetation (EVI), chlorophyll concentration in vegetation (GNDVI), relative abundance of chlorophyll (MCARI1), cover of the vegetation with eliminated soil-induced variations in vegetation indices (OSAVI) and measure of the efficiency of which vegetation can use incident light for photosynthesis (SIPI3).
• Based on the results of the study and the literature review, the act of monitoring shrubby formation is ecologically important for monitoring water and for carbon management, as well as for preventing fires and soil erosion.