Land Use / Land Cover ( LULC ) Using Landsat Data Series ( MSS , TM , ETM + and OLI ) in Azrou Forest , in the Central Middle Atlas of Morocco

The study of land use/land cover (LULC) has become an increasingly important stage in the development of forest ecosystems strategies. Hence, the main goal of this study was to describe the vegetation change of Azrou Forest in the Middle Atlas, Morocco, between 1987 and 2017. To achieve this, a set of Landsat images, including one Multispectral Scanner (MSS) scene from 1987; one Enhanced Thematic Mapper Plus (ETM+) scene from 2000; two Thematic Mapper (TM) scenes from 1995 and 2011; and one Landsat 8 Operational Land Imager (OLI) scene from 2017; were acquired and processed. Ground-based survey data and the normalized difference vegetation index (NDVI) were used to identify and to improve the discrimination between LULC categories. Then, the maximum likelihood (ML) classification method was applied was applied, in order to produce land cover maps for each year. Three classes were considered by the classification of NDVI value: low-density vegetation; moderate-density vegetation, and high-density vegetation. Our study achieved classification accuracies of 66.8% (1987), 99.9% (1995), 99.8% (2000), 99.9% (2011), and 99.9% (2017). The results from the Landsat-based image analysis show that the area of low-density vegetation was decreased from 27.4% to 2.1% over the past 30 years. While, in 2017, the class of high-density vegetation was increased to 64.6% of the total area of study area. The results of this study show that the total forest cover remained stable. The present study highlights the importance of the image classification algorithms combined with NDVI index for better understanding the changes that have occurred in this forest. Therefore, the findings of this study could assist planners and decision-makers to guide, in a good manner, the sustainable land development of areas with similar backgrounds.


Introduction
Land use/land cover (LULC) changes are among the most important applications of Earth Observation (EO) satellite sensor data [1][2][3].It provides a comprehensive and good understanding of Environments 2018, 5, 131; doi:10.3390/environments5120131www.mdpi.com/journal/environmentsecosystem monitoring and functioning, and responses to environmental factors [4][5][6][7][8][9][10].Remote sensing techniques have been recognized as a powerful means to obtain information on Earth's surface features [11] at different spatial and temporal scales [12].Since the beginning of the 1970s, scientists have used various types of remote sensing data, acquired by the Landsat series of satellites, in the identification of LULC changes [13][14][15][16][17].More recently, several researchers have concentrated on evaluating the potential of satellites data on land use classification, monitoring the land use or quantifying and analyzing land use changes [10].
Yuan et al. [18] used multitemporal Landsat TM data to map and monitor land cover change in the seven-county Twin Cities Metropolitan Area of Minnesota for 1986Minnesota for , 1991Minnesota for , 1998Minnesota for , and 2002.The study showed an increase of urban land from 23.7% to 32.8% of the total area between 1986 and 2002, while rural cover types of agriculture, forest, and wetland decreased from 69.6% to 60.5%.Demissie et al. [19] used Landsat data from 1973 to 2015 to assess LULC changes and their causes in Libokemkem District of South Gonder, Ethiopia.Their study showed that about 60.1% of the land experienced changes in LULC in 42 years.
In Morocco, several works have used various types of remote sensing data on LULC analysis.Barakat et al. [20] have used ASTER (Advanced Spaceborne Thermal Emission and Reflection Radiometer) and Sentinel-2A MSI images acquired in 2001 and 2015, respectively, for quantifying the changes in the Eastern area of Béni-Mellal province, based on the supervised classification algorithm and the normalized difference vegetation index (NDVI).The study showed that changes in land use have resulted in the increase of forest area.De Waroux et al. [21] estimated forest changes in argan woodlands (Morocco) using aerial photographs and satellite images between 1970 and 2007.They found that forest density decreased by 44.5% during the period of study.Hammi et al. [22] studied land cover changes in Aït Bouguemez Valley located in the middle of the Central High Atlas Mountains (Morocco) in the Azilal province, using aerial photographs dated from 1964 and Spot 5 satellite image from 2002.They found that, in the past 38 years, forest ecosystems have been affected by a relative decrease of 20.7% of the total forest area, and 8.7% for the mean canopy cover.
Several methods have been used for analyzing LULC, such as conventional image differencing, multidate image classification, image differencing rationing, vegetation index differencing, principal components analysis, and change vector analysis [23].For more details, see review papers from Lu et al. [24] and Coppin et al. [25].
The forest ecosystem in the Middle Atlas of Morocco has attracted the attention of several researchers (e.g., [26][27][28][29]), because of the enormous value it provides for the local population in their everyday needs (wood, timber, food, etc.), as well as for local authorities and urban communities.Despite those values, this ecosystem experienced major degradation on its biodiversity (fauna and flora) caused by both natural processes as well as human activities.This requires a depth of knowledge of ecosystem services provided by this forest ecosystem, as well as knowing the interaction between human and natural phenomena, to better manage and use resources in a good way.This is the role of the National Forest Inventories.However, they are costly and time-consuming [30].
Among other forests in Morocco, the Azrou forest has been gaining popularity as a tourist attraction, due to its rich biodiversity and ecosystem values.It is dominated by the Atlas cedar (Cedrus atlantica) known as a noble tree species [31], representing an important resource that underpins the progress of socioeconomic development and benefits to the local population and the country, in general.However, the study area is subjected to several environmental disequilibria, including soil erosion, overgrazing, cutting of trees, etc.Along with these, the cedar forest is altered and affected by an advanced stage of degradation, according to [32].For these reasons, the Azrou forest was selected for this research to create up-to-date LULC maps, which could be used as scientific data for efficient decision-making in the forest ecosystem.Keeping all the above reasons in view, we have integrated remote sensing data, Geographic Information System (GIS) analysis, and field surveys, in order to evaluate the potential of Landsat imagery to monitor spatiotemporal changes in land cover in the Azrou forest.Thus, to obtain the different land use classes, the NDVI was used to identify and to improve discrimination between LULC categories throughout the study area.Then, the maximum likelihood classification (MLC) supervised classification method was applied in this study because of its availability, and it does not require an extended training process [33,34].
The advantage of MLC method as a parametric classifier is that it takes into account the variance-covariance within the class distributions and, for normally distributed data, the MLC performs better than the other known parametric classifies (e.g., decision trees) [35].The results were validated with control points (ground truth).
The novelty of our study relies on the methodology proposed for the update of LULC information.Mapping and monitoring LULC, and their changes, is of critical and valuable importance for sustainable development and for monitoring environment change.Overall, the findings of this study could assist planners and decision-makers to guide, in a good manner, the sustainable land development of areas with similar backgrounds.

Materials and Methods
The study area is the Azrou forest; it covers 24,588 ha in the Middle Atlas of Morocco (Figure 1).The vegetation consists of evergreen cedar (Cedrus atlantica)-pure stands or mixed with evergreen stands of holm oak (Quercus ilex L., and/or with Quercus canariensis).The soils in the study area are formed by limestone, dolomite, and basalt [36].The climate of the region is Mediterranean, with a dry season from May to October, and a rainy season is from November to April.Ifrane Meteorological Station documents average annual rainfall from 563 mm to 1122 mm, and the maximum temperature of 30.3 • C. At Tagounit Meteorological Station, the maximum temperature reached 43 • C [36].The study area is characterized by an altitude changing from 856 to 2422 m, where the cedar forests are located in an altitude higher than 1600 m, and the holm oak in an altitude ranging between 1200 and 1600 m.Despite the diversity of services ecosystems of Azrou forest, the Atlas cedar decline symptoms go back to the 1940s, when scenarios of severe drought have been caused a mortality of this species [26,32,37].According to [26], an estimated 23% of the Atlas forest surface area (3400 of 15,000 hectares) showed some decline symptoms after the 1999-2000 drought events.
supervised classification method was applied in this study because of its availability, and it does not require an extended training process [33,34].
The advantage of MLC method as a parametric classifier is that it takes into account the variance-covariance within the class distributions and, for normally distributed data, the MLC performs better than the other known parametric classifies (e.g., decision trees) [35].The results were validated with control points (ground truth).
The novelty of our study relies on the methodology proposed for the update of LULC information.Mapping and monitoring LULC, and their changes, is of critical and valuable importance for sustainable development and for monitoring environment change.Overall, the findings of this study could assist planners and decision-makers to guide, in a good manner, the sustainable land development of areas with similar backgrounds.

Materials and Methods
The study area is the Azrou forest; it covers 24,588 ha in the Middle Atlas of Morocco (Figure 1).The vegetation consists of evergreen cedar (Cedrus atlantica)-pure stands or mixed with evergreen stands of holm oak (Quercus ilex L., and/or with Quercus canariensis).The soils in the study area are formed by limestone, dolomite, and basalt [36].The climate of the region is Mediterranean, with a dry season from May to October, and a rainy season is from November to April.Ifrane Meteorological Station documents average annual rainfall from 563 mm to 1122 mm, and the maximum temperature of 30.3 °C.At Tagounit Meteorological Station, the maximum temperature reached 43 °C [36].The study area is characterized by an altitude changing from 856 to 2422 m, where the cedar forests are located in an altitude higher than 1600 m, and the holm oak in an altitude ranging between 1200 and 1600 m.Despite the diversity of services ecosystems of Azrou forest, the Atlas cedar decline symptoms go back to the 1940s, when scenarios of severe drought have been caused a mortality of this species [26,32,37].According to [26], an estimated 23% of the Atlas forest surface area (3400 of 15,000 hectares) showed some decline symptoms after the 1999-2000 drought events.

Data Source and Pre-Processing
The Landsat satellite images (path 201; raw 037) include Multispectral Scanner (MSS), Enhanced Thematic Mapper Plus (ETM+), Thematic Mapper (TM), and Operational Land Imager (OLI).All images were geometrically corrected and acquired in level 1T (L1T).They were acquired free of charge from the United States Geological Survey (USGS) website [38].Table 1 shows the Landsat scenes used in the study.To properly compare the NDVI differences, all the images were acquired between July and September, as this is the time of year where different kinds of vegetation are all at a stable stage.In addition, the time gap between all the satellite images was more than 16 days, because of cloudiness or technically damaged scenes.Two criteria were followed to choose the satellites images in this study, according to [39]: (1) the satellite images must have less than 10% cloud coverage (if possible, cloud free); (2) the images satellite series should be available for a long time series.The image processing chain was performed using the open source Quantum GIS (QGIS) software (Version QGIS 3.0.0)[40].Atmospheric correction was applied through the Semi-Automatic Classification Plugin [41] based on the Dark Object Subtraction (DOS) algorithm [42].Google Earth images and field data based on map forest stand [28,35] were used to validate the proposed methodology.Field work was carried out using ground control points (GPSs) for ground truth surveys.In addition, baseline data were obtained from a published map of forest stand of the study area, and a previous knowledge and a deep discussion with local expert and stakeholders to recognize the LCLU classes.A flowchart with the methodology followed in this study is presented in Figure 2.
For all images data, the calibration was achieved by converting raw digital numbers (DN) into at-sensor spectral radiance [43,44], through Equation (1): where L λ is the spectral radiance at the sensor's aperture in mW/(cm

NDVI Classification
Based on the values of NDVI and in the interpretation of Google Earth images to identify and to quantify the vegetation change during the analysis period (1987-2017), we applied an NDVI thresholding [45,46] to classify NDVI into three classes: (1) NDVI values below to 0.2 were considered as low-density vegetation (2) NDVI values between 0.2 and 0.5 were considered as moderate-density vegetation and (3) NDVI values higher than 0.5 were considered as high-density vegetation.This classification depends on the soil type and, in this case, is based on the knowledge of the study area characteristics.
NDVI, proposed by Rouse Jr. et al. [47], is a widely used vegetation index and provides information on the amount and vigor of vegetation considering the near infrared (NIR) and visible red bands of the electromagnetic spectrum [33,48], calculated through Equation ( 3): Where NIR represents the near-infrared reflectance, and RED represents the reflectance in the red band.NIR and RED are the reflectance radiated in the near-infrared band 7 (0.8-1.1 μm) and the visible red 5 (0.6-0.7 μm) wave band of the satellite MSS Landsat 5, respectively.For TM and ETM+ sensors, the NIR band is band 4, and the RED band is band 3. Finally, for the OLI sensor, the NIR and RED bands are band 5 (0.85-0.88 μm) and band 4 (0.64-0.67 μm), respectively.The NDVI value ranges from −1 to 1.The highest value represents healthy vegetation, while the lowest NDVI value indicates non-vegetative cover.

NDVI Classification
Based on the values of NDVI and in the interpretation of Google Earth images to identify and to quantify the vegetation change during the analysis period (1987-2017), we applied an NDVI thresholding [45,46] to classify NDVI into three classes: (1) NDVI values below to 0.2 were considered as low-density vegetation (2) NDVI values between 0.2 and 0.5 were considered as moderate-density vegetation and (3) NDVI values higher than 0.5 were considered as high-density vegetation.This classification depends on the soil type and, in this case, is based on the knowledge of the study area characteristics.
NDVI, proposed by Rouse Jr. et al. [47], is a widely used vegetation index and provides information on the amount and vigor of vegetation considering the near infrared (NIR) and visible red bands of the electromagnetic spectrum [33,48], calculated through Equation ( 3): where NIR represents the near-infrared reflectance, and RED represents the reflectance in the red band.NIR and RED are the reflectance radiated in the near-infrared band 7 (0.8-1.1 µm) and the visible red 5 (0.6-0.7 µm) wave band of the satellite MSS Landsat 5, respectively.For TM and ETM+ sensors, the NIR band is band 4, and the RED band is band 3. Finally, for the OLI sensor, the NIR and RED bands are band 5 (0.85-0.88 µm) and band 4 (0.64-0.67 µm), respectively.The NDVI value ranges from −1 to 1.The highest value represents healthy vegetation, while the lowest NDVI value indicates non-vegetative cover.

Maximum Likelihood Classifier
The maximum likelihood (ML) classifier (supervised classification) has been widely used [33,34,49] because of its availability, and its does not require an extended training process [33,34].
The ML classifier is based on the probability that a pixel belongs to a particular class [35].Detailed explanations can be found in the literature [50][51][52][53][54] The algorithm for computing the weighted distance, or likelihood D of unknown measurement vector X, belongs to one of the known classes, Mc, that is based on the Bayesian Equation (4) [35]: where D is the weighted distance (likelihood); c is a particular class; X is the measurement vector of the candidate pixel; M c is the mean vector of the sample of class c; a c is the percent probability that any candidate pixel is a member of class c (defaults to 1.0, or is entered from a priori knowledge); Cov c is the covariance matrix of the pixels in the sample of class c; |Cov c | is the determinant of Cov c (matrix algebra); Cov c − 1 is the inverse of Cov c (matrix algebra) ln = natural logarithm function; and T is the transposition function (matrix algebra).In this study, ML classifier was used for the classification of the Landsat images.Based on field data and high-resolution imagery from Google Earth as a reference data, three LULC (Table 2) categories were identified: (1) cedar forest; (2) holm oak; and (3) bare soil.In the case of the MSS image, the training and validation data were generated from the visual interpretation of Google Earth images, from 1987, of the study area, since no field surveys were available for this period.The reference data samples were used to create ROI (region of interest) polygons in ENVI 5.0 software on the Landsat data, for the generation of the training areas.
Considering the training dataset, the number of pixels per each class cover selected in this study area is shown in Table 2. To evaluate the class separability of the training areas, the Jeffries-Matusita (JM) distance [51] was computed.The Jeffries-Matusita (JM) distance is the average distance between two class density functions [55][56][57], and it is calculated through Equations ( 5) and (6): where i and j are the two classes being compared, m i and m j is the mean vector of class i and j, respectively, and C i and C j is the covariance matrix of the feature response of class i and j, respectively.The value of 2 suggest complete separability.
Analyzing the pairs' separability between classes, it is possible to conclude that the ROIs selected in this study present good separability, as shown on Table 3 (most of the values are higher than 1.9).The only case that not fulfil the good separability requirements is the MSS image (1987) for the ROI pair cedar forest/holm oak.This fact could be justified by the low spatial resolution of the MSS image and lower classification accuracy.

Accuracy Assessment
The accuracy assessment is an important step for evaluation of the result of the classification process [58,59], as the user of land-cover output needs to know how accurate the result is, in order to use the data correctly [60][61][62].In this work, we consider that the minimum level of interpretation accuracy in the identification of land use and LULC categories from remote sensing data should be at least 85%, as recommended by different authors [52,63,64].A confusion matrix (or error matrix) [65], is a square array of numbers set out in rows and columns, which expresses the number of sample units (i.e., pixels, clusters of pixels, or polygons) assigned to a particular category relative to the actual category, as verified on the ground.An error matrix is widely used for classification accuracy [58,60,64,65], to derive a series of descriptive and analytical statistics [60,64].Several studies have explained, in detail, the meanings and calculation methods for overall accuracy (OA), producer's accuracy (PA), user's accuracy (UA), and kappa coefficient [58,60,[63][64][65].

Land Uses/Land Cover Changes Monitoring
LULC changes were estimated from 1987 to 2017 (Table 4).Holm oak was the dominant land class in 2017, with 37.1% of the occupancy.Cedar forest class changed slightly from 1987 to 2017.Bare soil decreased from 45.5% to 32.4% (Table 4).Classified maps, showing the LULC categories in Azrou forest, are shown in Figure 3.

Accuracy Assessment
The classification accuracy was evaluated through the confusion matrix.The classified images showed an OA of 66.8%, 99.9%, 99.8%, 99.9%, and 99.9% in 1987, 1995, 2000, 2011, and 2017 images, respectively, with a kappa statistic of 0.413, 0.974, 0.995, 0.998, and 0.999, respectively (Table 5).User's and producer's accuracies of individual classes for each year of the land cover map are presented in
User's and producer's accuracies of individual classes for each year of the land cover map are presented in Table 5, and indicate that all classes, with the exception of holm oak, have user's and producer's accuracies higher than 80%.The results presented in Figure 4 and Table 6 reveal that the class of moderate-density vegetation was the dominant NDVI land cover class.It covers approximately 43.8% of the region in 1987.However, it declined to 33.2% in 2017.Low-density vegetation occupied the smallest portion of the region, and is was predominantly situated on the boundaries of the study area (Figure 4).The area of high-density vegetation class increased from 28.7% in 1987, to 64.6% in 2017.

Discussion
The classification method that was employed in this study was very effective in improving LULC classification accuracy.
As expected, and as had been discussed in several works [66][67][68], for MSS classification results, cedar forest and holm oak did not show a better separability among the LULC types, resulting in classification accuracy.
Similarly, Pang et al. [69], who studied deforestation and resulting land use change patterns by analysis of Landsat satellite images from 1979 (MSS), 1992 (TM), 2001 (ETM+), and 2006 (TM) in Suan County, Hwanghae Province, Democratic People's Republic of Korea (DPR Korea), found that the 1979 Landsat MSS image was more difficult to classify than the other Landsat images.They explain that this is mainly due to the comparative lower spectral resolution of the sensor.
Ozesmi and Bauer [70] have confirmed that the MLC leads to higher accuracy than the decision trees (DT) classifier.
Moreover, our findings can be explained by a reinforced observation by Lu et al. [71] that the traditional approach to classification (such as ML) only distinguishes clearly between forest and non-forest land use and land covers.
Previous studies have successfully been carried out using remote sensing geospatial tools in a related study at the same profile area.For example, Barakat et al. [20] have used ASTER (Advanced Spaceborne Thermal Emission and Reflection Radiometer) and Sentinel-2A MSI images, acquired in 2001 and 2015 respectively, for quantifying the changes in the Eastern area of Béni-Mellal province based on the supervised classification algorithm and the normalized difference vegetation index (NDVI).The study showed that changes in land use have resulted in the increase of forest area.Haboudane and Bahri [72] have used multidate satellite images (Thematic Mapper (TM), Enhanced Thematic Mapper Plus (ETM+), and Advanced Spaceborne Thermal Emission and Reflection Radiometer (ASTER)) using maximum likelihood classification (MLC and spectral mixture analysis (SMA) for mapping major forest species in the Moroccan Middle Atlas.The results showed a decrease in cedar area of 12%, and an increase of 8% in cedar forest; in oak forest, the increment (21%) had exceeded the deforestation effect (17%); while deciduous was either degraded (11%) or remained stable (21%).
In summary, from the analyses of LULC changes, it could be found that the total forest cover remained stable and, thus, NDVI values were probably related with good climatic conditions during the period of study, as well as a result of huge efforts of the Moroccan High Commission for Water, Forests, and Combating Desertification (HCWFCD).The authorities have implemented extensive programs and a large portfolio of projects in the study area (e.g., development and protection of forest areas in the province of Ifrane and Participatory Rural Development in the Central Middle Atlas (HCWFCD 2017 [73,74].Moreover, the Moroccan High Commissariat for Water, Forestry and Combating Desertification (HCWFCD) have undertaken huge efforts for studying and regeneration of the forest ecosystem, especially via reforestation [36,75].The increase in forest area is mainly through plantations undertaken by the HCWFCD, for example; since 2011, the rate of reforestation and silvopastoral improvement achieve 42,988 ha per plantation, including (i) 24,464 ha of reforestation until 24 February 2012; (ii) 12,804 ha of forest regeneration (cedar-4038 ha, cork oak-2107 ha, cedar-2403 ha, argan tree-1851 ha, Acacia raddiana-1205 ha, Aleppo pine-500 ha, holm oak-400 ha, (Aleppo Pine + Thuja)-100 ha and (Thuja + Holm oak)-200 ha); (iii) 5720 ha of silvopastoral improvement.Moreover, the management treatments, which were introduced during 2003-2007, involved 167.500 ha, or nearly 33.500 ha/an, specifically, the regeneration of natural species: cedar, atlas cypress, cork oak, carob, and argan (for more information, please see [75,76] Additionally, since 2016, the Atlas cedar forests in Morocco are included in the world network Biosphere Reserve established by the United Nations Educational, Scientific, and Cultural Organization UNESCO 'S Man and the Biosphere Program since 2001.This will be an effective step to regenerate these species, while maintaining essential ecosystem functions.Overall, the methodology proposed in this work could be useful for planning the use of forest resources and for developing strategies to maintain the vitality of the studied forest.

Conclusions
The study set out to determine LCLU statuses over the period 1987 to 2017 for the Azrou forest in the Middle Atlas of Morocco.Based on remote data, MLC, GIS geospatial tools, and including the use of NDVI and ancillary data (field data), produced up-to-date reference information in the Azrou forest.
The objective of this study was to provide a recent perspective for LULC changes that have taken place, in the last 30 years, in Azrou forest.The results of this work can be summarized as follows: (i) it was found that a considerable increase in vegetation cover has taken place in the study area; (ii) the area of bare soil has decreased considerably during the period of study; (iii) our study achieved high classification accuracies (OA) for all the images, with the exception of MSS image.
The presented methodology was very promising for being used in such areas, based on remote sensing data and geospatial tools.With this, we hope to provide new insights for guiding planners in decision-making and ecosystem conservation.
This study demonstrates that exploiting a free archive of Landsat data, and its processing through open source software, provides an accurate approach to mapping and analyzing changes in land cover over time, that can be used as a guide to land management and policy decisions.In future research, we plan to compare the results of using Landsat data with other types of remote sensing data (e.g., Sentinel 2) to monitor LULC in the study area.However, Sentinel 2 data are only available since 2015.Furthermore, the proposed methodology will be applied to forestlands in other regions to advance a more comprehensive monitoring and mapping of forest change across the Middle Atlas of Morocco.

Figure 1 .
Figure 1.Location of the study area (elevation data).

Figure 1 .
Figure 1.Location of the study area (elevation data).

Figure 2 .
Figure 2. Flowchart of the methodology adopted.

Figure 2 .
Figure 2. Flowchart of the methodology adopted.

Table 1 .
Landsat scenes used in the study.
2•sr•µm); DN is the digital number of the quantized calibrated pixel value; L λmax is the maximum spectral radiance that is scaled to QCAL λmax ; L λmin is the minimum spectral radiance that is scaled to QCAL λmin ; QCAL λmax is the maximum quantized calibrated pixel value in DN (corresponding to L λmax ); and QCAL λmin is the minimum quantized calibrated pixel value in DN (corresponding to L λmin ).After this conversion, the radiance is converted into top-of-atmosphere (TOA) reflectance using Equation (2): Environments 2018, 5, x FOR PEER REVIEW 5 of 16

Table 2 .
The number of pixels used for classifications.

Table 4 .
Land-cover (in hectare and percentage) from 1987 to 2017.

Table 5 .
Accuracy assessment of the land cover maps generated.

Table 6 .
Normalized difference vegetation index (NDVI) land cover change (in km 2 and percentage) from 1987 to 2017.

Table 6 .
Normalized difference vegetation index (NDVI) land cover change (in km 2 and percentage) from 1987 to 2017.