Mapping Forest Species in the Central Middle Atlas of Morocco ( Azrou Forest ) through Remote Sensing Techniques

The studies of forest ecosystems from remotely-sensed data are of great interest to researchers because of ecosystem services provided by this ecosystem, including protection of soils and vegetation, climate stabilization, and regulation of the hydrological cycle. In this context, our study focused on the use of a spectral angle mapper (SAM) classification method for mapping species in the Azrou Forest, Central Middle Atlas, Morocco. A Sentinel-2A image combined with ground reference data were used in this research. Four classes (holm oak, cedar forest, bare soil, and others-unclassified) were selected; they represent, respectively, 27, 11, 24, and 38% of the study area. The overall accuracy of classification was estimated to be around 99.72%. This work explored the potential of the SAM classification combined with Sentinel-2A data for mapping land cover in the Azrou Forest ecosystem.


Introduction
Forests play a range of services that are crucial for humans and different kinds of fauna.They can provide food and fiber; regulate the hydrological cycle; protect watersheds and their vegetation, water flows, and several types of ecological services [1][2][3][4][5].However, natural equilibria in ecosystem services provided by our forests have been dramatically reduced because of natural and human-induced causes [6].
The Food and Agriculture Organization for United Nations [7] mentions that global forest area fell by 129 million hectares (3.1%) in the period 1990-2015, to just under 4 billion hectares.
Forest cover since the 1950s in the Mediterranean area has been lost at a rate of about 30% [8].However, according to Hansen et al. [9] the dry tropical biome represented 20% of the total forest cover in 2000, and were being lost at rate of 2.9% between 2000 and 2005.In order to protect and to use those services, the National Forest Inventories have been used and applied in many countries in the world.Despite the role of National Forest Inventories for obtaining information on forest ecosystems, previous studies have shown that they are costly and time-consuming [10].
To overcome this, several authors have made serious efforts to seek and develop new approaches for the identification of degraded areas in forests [11,12], to map and classify different forest species through the use of satellite imagery [13][14][15].Several classification methods have been applied in land cover mapping, commonly distinguished in a first approach: those focused on artificial intelligence, such as artificial neural networks (ANN) [16] and those oriented towards supervised or unsupervised algorithms, including maximum likelihood (ML) classification [17], decision tree (DT) classifier [18][19][20], the k-nearest neighbors (K-NN) method [21], and support vector machine (SVM) [22][23][24].Among those algorithms, our approach was based on the spectral angle mapper (SAM) because of its quickness and ease of use, independence from illumination variations, as well as its use in several studies [25].In addition, according to Bonn [15], mapping and monitoring species in forest ecosystems of the Middle Atlas Mountains of Morocco is complicated due to considerable physical and environmental variability; both at spatial and spectral scales.This variability limits the use of the conventional methods and spectral indices for this task.To better overcome these challenges, they have used ASTER (Advanced Spaceborne Thermal Emission and Reflection Radiometer) images to compare the performances of linear spectral mixture analysis (LSMA) and SAM for mapping and discrimination of species in this forest ecosystem.Hence, Sentinel-2A with high spatial, spectral, and temporal resolution combined with SAM can be used to improve discrimination of species in such areas of study.
Several authors employed different remotely-sensed data and different algorithms to map forest types in non-tropical areas.For instance, Shen et al. [26] investigated and demonstrated the mapping potential of the forest ecosystem at the tree species level from high spatial resolution hyperspectral images (Airborne Imaging Spectrometer for Applications-AISA), in Hachioji, Japan.The mapping performance of eight conventional classification methods were tested-including SAM, which achieved the best result.Additionally, Kachmar et al. [27] used Landsat 5 TM data, in order to classify dominant forest cover types in the Naeba Mountains of Japan considering the SAM classifier.Silva [28] used a dataset composed by a several Landsat-5 TM and Rapid-Eye images from the Paraíba Valley (Brazil).They used training data from 2011 to build a spectral reference library considering the Landsat-5 TM spectral signatures for each land use/land cover LULC (eight classes were considered).The spectral library was then used to classify the time series of Landsat-5 TM images (years 1985, 1994, 1995, 2005, and 2011) with SAM.Troyer et al. [29] aimed to create a higher spatial resolution LULC dataset for the entire Little Miami River watershed in Southeast Ohio for a number of studies being conducted by the U.S. Environmental Protection Agency (U.S. EPA).The LULC classification was derived from 82 flight lines of Compact Airborne Spectrographic Imager (CASI) hyperspectral imagery.They used several classification algorithms, including maximum likelihood classification (MLC), SAM, and the classification and regression tree (CART), however, no single algorithm alone proved to be capable of successfully classifying the hyperspectral data.Nevertheless, SAM produced the best results.However, in the works previously cited, none of them used Sentinel-2A data combined with SAM.SAM treats the spectra as vectors in a space with dimensionality equal to the number of bands used.Using the SAM classifier and spectrally-suitable forest training areas, forest cover types are classified and their accuracies are related to topographic correction methods, applied and localized land use, and land cover change occurring in the study area [27].Moreover, SAM's advantage over more 'traditional' classifiers is its relative insensitivity to illumination and albedo effects inherent with remotely-sensed imagery [30].
Until recently, the availability of Sentinel-2A as an instrument with high spatial, spectral, and temporal resolution has opened up other opportunities for this aim.This is the novelty presented in this work, the combination of Sentinel-2A data with SAM in order to map forest species in a non-tropical area.
The aim of this work was to map and classify the dominant forest species in the Azrou Forest in the Central Middle Atlas region of Morocco by using high spatial resolution satellite images (Sentinel -2A), based on the spectral angle mapping (SAM) supervised classification method.In this context, our methodology represents a powerful step because the use of SAM for this focus with Sentinel-2A has not yet been fully explored and has not been tested, until now, in this study area.Furthermore, this approach has valuable implications for forest cover mapping in other areas, specifically where field surveys and field data are costly and time consuming.

Study Area
The study area is located in the Central Middle Atlas of Morocco (Figure 1), between 5 • 34 1.415" W/5 • 52 5.464" W and 37 As a result of forest management plan, the species composition of the area is highly heterogeneous, mainly covered by evergreen cedar (Cedrus atlantica); pure stands or mixed with evergreen stands of (holm) oak (Quercus ilex L, and/or with Quercus canariensis).
Cedrus atlantica, known as a noble species, representing an important ecological and economic value in Morocco, covers an area of over 130,000 ha distributed in Morocco (Rif, Middle Atlas, and northeast of the High Atlas) and Algeria [31].It grows in cool, moist environments at elevations between 1300 and 2600 m in the Rif, Middle Atlas, and High Atlas mountains [31][32][33], where the amount of annual rainfall ranges from 500 to 2000 mm and the minimum temperature of the coldest month is between −8 • C and −1 • C [31,34].Whereas evergreen oak forest constitutes the natural vegetation under semi-arid climates [35].
Still, according to the same plan, the Azrou forest is composed of all types of substrates (limestone, dolomite, and basalt) combined with a wide variety of soils ranging from deep soils to shallow rocky soils.

Study Area
The study area is located in the Central Middle Atlas of Morocco (Figure 1), between 5°34′1.415″W/5°52′5.464″W and 37°1′25.700″N/36°59′27.554″N, and it covers approximately 24,588 ha.The climate is typically Mediterranean with annual rainfall varying from 563 mm to 1122 mm, while the maximum temperature is around 30.3 °C at Ifrane Meteorological Station and 43 °C at Tagounit Meteorological Station.July and August are the driest months in the study area.
As a result of forest management plan, the species composition of the area is highly heterogeneous, mainly covered by evergreen cedar (Cedrus atlantica); pure stands or mixed with evergreen stands of (holm) oak (Quercus ilex L, and/or with Quercus canariensis).
Cedrus atlantica, known as a noble species, representing an important ecological and economic value in Morocco, covers an area of over 130,000 ha distributed in Morocco (Rif, Middle Atlas, and northeast of the High Atlas) and Algeria [31].It grows in cool, moist environments at elevations between 1300 and 2600 m in the Rif, Middle Atlas, and High Atlas mountains [31][32][33], where the amount of annual rainfall ranges from 500 to 2000 mm and the minimum temperature of the coldest month is between −8 °C and −1 °C [31,34].Whereas evergreen oak forest constitutes the natural vegetation under semi-arid climates [35].
Still, according to the same plan, the Azrou forest is composed of all types of substrates (limestone, dolomite, and basalt) combined with a wide variety of soils ranging from deep soils to shallow rocky soils.

Data Description
A Sentinel-2A image acquired on 28 August 2015 was used in this study.It consists of 13 spectral bands from the visible to the SWIR; of which 4 bands have a resolution of 10 m, 6 bands have a resolution of 20 m and 3 bands have a resolution of 60 m.

Pre-Processing Stage
The pre-processing of the Sentinel 2A image was performed under QGIS software, through the Semi-Automatic Classification Plugin (SCP) interface, developed by Luca Congedo [36].The pre-

Pre-Processing Stage
The pre-processing of the Sentinel 2A image was performed under QGIS software, through the Semi-Automatic Classification Plugin (SCP) interface, developed by Luca Congedo [36].The pre-processing stage consists of several steps, from the image download to the final image (map).The first step consists of converting the digital number (DN) to the top of atmospheric (TOA) reflectance, followed by the atmospheric correction, where the dark object subtraction (DOS) was considered.The DOS method assumes that atmospheric path radiance is the radiance value measured by the satellite for the darkest object within the image, usually clear water bodies or areas in complete shadow [37].This model corrects only the atmospheric additive scattering component.

Supervised Classification
The SAM classification algorithm was used in this study.As a supervised classification technique, it is highly dependent on the identification of the training areas obtained from the observation of a field spectrometer, or are taken directly from a remote sensing image with sufficient field data, or from spectral libraries [30,38].Developed by Kruse et al. [30], SAM determines rapidly and easily the spectral similarity between the image spectra to reference reflectance spectra (Figure 2) [38][39][40][41] and it is based on the number of bands used in the processed image [42][43][44].Small angle values indicate greater similarity between pixel and reference spectra [45].
This method is relatively insensitive to changes in illumination in the scene [46].Mathematically, SAM can be expressed by Equation ( 1) where α is the angle formed between the reference spectrum and image spectrum; X is the image spectrum; and Y is the reference spectrum.Each pixel will be assigned to the class according to the lowest spectral angle value pixels are similar.Pixels with a measurement greater than the specified maximum divergence threshold are not classified [47].
ISPRS Int.J. Geo-Inf.2017, 6, 275 4 of 10 processing stage consists of several steps, from the image download to the final image (map).The first step consists of converting the digital number (DN) to the top of atmospheric (TOA) reflectance, followed by the atmospheric correction, where the dark object subtraction (DOS) was considered.
The DOS method assumes that atmospheric path radiance is the radiance value measured by the satellite for the darkest object within the image, usually clear water bodies or areas in complete shadow [37].This model corrects only the atmospheric additive scattering component.

Supervised Classification
The SAM classification algorithm was used in this study.As a supervised classification technique, it is highly dependent on the identification of the training areas obtained from the observation of a field spectrometer, or are taken directly from a remote sensing image with sufficient field data, or from spectral libraries [30,38].Developed by Kruse et al. [30], SAM determines rapidly and easily the spectral similarity between the image spectra to reference reflectance spectra (Figure 2) [38][39][40][41] and it is based on the number of bands used in the processed image [42][43][44].Small angle values indicate greater similarity between pixel and reference spectra [45].
This method is relatively insensitive to changes in illumination in the scene [46].Mathematically, SAM can be expressed by Equation ( 1) where α is the angle formed between the reference spectrum and image spectrum; X is the image spectrum; and Y is the reference spectrum.Each pixel will be assigned to the class according to the lowest spectral angle value pixels are similar.Pixels with a measurement greater than the specified maximum divergence threshold are not classified [47].
Figure 2. Spectral angle formed between the reference spectrum and the image spectrum (adapted from [30]).

Results
The SAM classification algorithm was applied to obtain a land cover classification map (Figure 3) from the high spatial resolution satellite image Sentinel-2A considered in this work.Figure 4 shows the spectral responses of the different classes selected extracted directly from the Sentinel-2A image.
In this study, 7419 pixels were considered as bare soil, 29,147 pixels as cedar forest, and 33,668 pixels as holm oak.The separability scores of the ROIs were calculated using the Jeffries-Matusita (J-M) distance between the ROIs.Therefore, all ROIs used in this study have good separability and are reasonable for use as training samples, as shown in Table 1.

Results
The SAM classification algorithm was applied to obtain a land cover classification map (Figure 3) from the high spatial resolution satellite image Sentinel-2A considered in this work.Figure 4 shows the spectral responses of the different classes selected extracted directly from the Sentinel-2A image.
In this study, 7419 pixels were considered as bare soil, 29,147 pixels as cedar forest, and 33,668 pixels as holm oak.The separability scores of the ROIs were calculated using the Jeffries-Matusita (J-M) distance between the ROIs.Therefore, all ROIs used in this study have good separability and are reasonable for use as training samples, as shown in Table 1.The output classes considered in the classification process: holm oak, cedar forest, bare soil, and others.Of these, cedar forest makes up 11% (2747 ha) of the area; this mainly occurs in the northeast and southwest of the study area.However, holm oak makes up approximately 27% (6603 ha) of the    The output classes considered in the classification process: holm oak, cedar forest, bare soil, and others.Of these, cedar forest makes up 11% (2747 ha) of the area; this mainly occurs in the northeast and southwest of the study area.However, holm oak makes up approximately 27% (6603 ha) of the  The output classes considered in the classification process: holm oak, cedar forest, bare soil, and others.Of these, cedar forest makes up 11% (2747 ha) of the area; this mainly occurs in the northeast and southwest of the study area.However, holm oak makes up approximately 27% (6603 ha) of the land cover; this is represented throughout the region.Bare soil class accounts for approximately 24% of the Azrou forest and the unclassified class presents 38% of the land cover in the area.
Accuracy assessment is an important step in the classification process.The main objective is to determine how effectively pixels were grouped into the correct classes.A confusion matrix is a square array of numbers set out in rows and columns that express the number of sample units assigned to a particular category relative to the actual category as verified by ground truth information usually collected from ground visits, or aerial photographs and satellite images, or a reference dataset.
The kappa statistic gives a measure that indicates whether the confusion matrix is significantly different from a random result.Kappa analysis can also be used to compare different matrices from different classifiers and to determine whether one result is significantly better than another [48].
Classification accuracy statistics are summarized in Table 2.The overall accuracy (OA) of classification was 99.72% and the kappa statistics was 0.99.A confusion matrix (Table 2) reveals that the forest classes were mapped with high accuracy.8,286,532 out of 8,309,854 pixels were correctly classified.It further shows that the overall accuracy (99.72%) and individual class accuracies are encouraging.The commission and omission errors of the four classes are reported in Table 2 (in percent).The unclassified (others) class was classified with commission and omission errors of 0.07% and 0.08%, respectively.For the bare soil class, commission, and omission errors are 1.12% and 1.06%, respectively.For the holm oak class commission and omission errors are 0.93% and 0.95%, respectively, with most of the commission (2.04%) and omission (2.06%) occurring with the cedar forest class.
According to studies carried out in 2007 by the High Commission for Water and Forests and the fight against desertification [49], pure cedar stands occupy an area of 1497.41 ha, i.e., approximately 8.41% of the total forest area, and pure oak stands cover an area of 4419.77ha, i.e., approximately 24.82% of the total area of the forest.The "other" section, which includes secondary species, represents a total area of 4392.45 ha, i.e., approximately 24.67% of the total area of the forest.
Comparing the present land cover SAM classification map to the study realized by the High Commission for Water, Forests, and for Combating Desertification in 2007, they also found the same results regarding the land cover map with similarity in the number of hectares for each class.
Our results and adopted method indicate that a SAM classification is a valuable and good method to map and classify forest regions.Moreover, our methodology represents a powerful step because the use of SAM for this focus with Sentinel-2A has not yet been fully explored and has not been tested until now in this study area.
In previous research, SAM has been successfully used [50][51][52][53].Chikhaloui et al. [50] explains that using SAM sometimes mixed the slightly-and the moderately-degraded soils into homogeneous regions because of the spectral similarity approach of this method [54].Additionally, this confusion can be due to the use of only one angular threshold (0.20 rad) for the different classes [50].
The results obtained by Matthew et al. [55] showed that the spectral angle mapper had poor performance (<51% OA) in comparison to linear discriminant analysis (LDA) and maximum likelihood (ML) classifiers for tropical rain forest trees.Similar results are found by Dhaval et al. [56].
Petropoulos et al. [57] found that the object-based classification outperformed the SAM by 7.91% OA for mapping land use/cover in the Mediterranean region.However, they noticed the required level of expertise and time in its implementation, as well as it being more expansive, in comparison to SAM.

Conclusions
The current study explored the potential of the SAM classification for mapping land cover in the Azrou forest ecosystem.The data used in this work are free and the software chosen for processing the data is open source, which allows other potential users to apply the same data and the same classification techniques to other locations.The results obtained from the application of the proposed method showed that the study area is mainly covered by four classes (holm oak, cedar forest, bare soil, and others-unclassified); they represent, respectively, 27, 11, 24, and 38% in the area of study.The OA of classification was estimated to be 99.72%.The results of this study show that the combination of Sentinel-2A data and SAM classification techniques may provide an inexpensive, easily implemented alternative to expensive, sample-based National Forest Inventories.Moreover, the results of this research, combined with other source data (e.g., field data) could be used in order to investigate the forest cover lost in the Mediterranean area.
• 1 25.700" N/36 • 59 27.554" N, and it covers approximately 24,588 ha.The climate is typically Mediterranean with annual rainfall varying from 563 mm to 1122 mm, while the maximum temperature is around 30.3 • C at Ifrane Meteorological Station and 43 • C at Tagounit Meteorological Station.July and August are the driest months in the study area.

Figure 1 .
Figure 1.Location of the study area (Central Middle Atlas of Morocco).

Figure 1 .
Figure 1.Location of the study area (Central Middle Atlas of Morocco).

A
Sentinel-2A image acquired on 28 August 2015 was used in this study.It consists of 13 spectral bands from the visible to the SWIR; of which 4 bands have a resolution of 10 m, 6 bands have a resolution of 20 m and 3 bands have a resolution of 60 m.

Figure 2 .
Figure2.Spectral angle formed between the reference spectrum and the image spectrum (adapted from[30]).

Figure 3 .
Figure 3. Land cover classification maps of the study area in 2015.

Figure 4 .
Figure 4. Spectral signatures based on the selection of pixels selected in the present study.

Figure 3 .
Figure 3. Land cover classification maps of the study area in 2015.

Figure 3 .
Figure 3. Land cover classification maps of the study area in 2015.

Figure 4 .
Figure 4. Spectral signatures based on the selection of pixels selected in the present study.

Figure 4 .
Figure 4. Spectral signatures based on the selection of pixels selected in the present study.

Table 2 .
Error matrix for SAM classification result in the study area.