Monitoring of forest cover dynamics in eastern area of Béni-Mellal Province using ASTER and Sentinel-2A multispectral data

Abstract

This study aimed to monitor and analyze the spatial and temporal dynamic of forest cover in Eastern area of Beni-Mellal Province (Morocco), using multispectral ASTER and Sentinel-2A MIS images acquired in 2001 and 2015, respectively. The supervised classification algorithm and NDVI were combined within a GIS environment to quantify the extent and density change of forest cover stands, i.e., Holm oak, Aleppo pine, Thya, Zea oak, Crops & others and Bare ground. The classification overall accuracy was 97.76 and 95.80% in 2001 and 2015 images, respectively. The result revealed an overall forest cover change with an increase in forested area. All species stands showed expansion at the expense of the bare ground and crops & others classes. The density maps showed a net density change with an expansion of dense forest class. The observed forest cover expansion may be due to the favourable climate in the examined period, the protection, the reforestation programs and the regeneration through clandestine cutting. These results constituted the first attempt at mapping and monitoring of forest cover change in the study region that used a remote sensing-based product. They will help authorities and forest managers for the development of sustainable forest conservation and management decisions.

Keywords:

• Béni-Mellal Province
• forest stand changes
• ASTER and Sentinel-2A data
• supervised classification

1. Introduction

Forest is a treasure trove of biodiversity and plays a vital role in maintaining the ecological balance and the health of our planet. The forest cover in arid and semi-arid ecosystems is important in support of global biodiversity, carbon sequestration, and economic activities. It provides such a key ecosystem goods and services to over one billion people living in the arid and semi-arid lands (Safriel et al., 2005). It also conserves many animal and plant species (Food and Agriculture Organization of the United Nations [FAO], 2015), characterized by many related functions such as timber production, climate regulation and recreation (Gamfeldt et al., 2013). Several scientific researches enumerated ecological, economic and social profit of forest for humanity, and above all for local populations (i.e. FAO, 2010; Karjalainen, Sarjala, & Raitio, 2010; Rotenberg & Yakir, 2010; Salehi & Karltun, 2010; Yang et al., 2012). The forest cover in arid and semi-arid regions (fragile environment) has undergone in the last few decades negative changes (de Waroux & Lambin, 2012). Its continued decline is due to rapid climate change and human pressures such as overexploitation, over-pasturage, conversion to agriculture, and urbanization. Even though this decline occurs in several arid areas in the world, and sometimes at a faster rate than to that of tropical forests (Seabrook, McAlpine, & Fensham, 2006; Zak, Cabido, & Hodgson, 2004), the forest cover in arid and semi-arid regions has received less attention (Grainger, 1999). International agreements on forest monitoring have begun to include forest degradation and deforestation (Kissinger, Herold, & De Sy, 2012). According to statistics from the Food and Agriculture Organization of the United Nations, global forest area decreased by 3.1% (129 million hectares) between 1990 and 2015 (FAO, 2016). The excessive exploitation of forest resources has led to depletion. Given this alarming situation of over-exploitation and degradation of forest heritage, it needs good management and sustainable development.

Assessment forest areas and their dynamics at local to regional scales has therefore become an urgent necessity given their vulnerability to anthropogenic and climatic pressures, importance for adaptive management and implementation of national directives and international treaties (United Nations Convention to Combat Desertification of 1994 (UNCCD); Convention on the International Trade in Endangered Species (CITES), Intergovernmental Forum on Forests (IFF), Evaluating their changes and species diversity over time and space is important to understand the effect of anthropogenic activities and climate change on forest cover. Despite the importance of forest inventories in obtaining information on changes in forest ecosystems, they are costly and time-consuming (Mohajane, Essahlaoui, Oudija, El Hafyani, & Cláudia Teodoro, 2017). An attractive and less-expensive method to carry out such monitoring is the remote sensing that has demonstrated its effectiveness for rapid assessment of forest dynamics and degradation over several decades and on a large-scale (Eva et al., 2010; Hansen et al., 2009; Nagendra et al., 2013; Wani, Joshi, Singh, & Shafi, 2016). For mapping and estimating forest cover change, a number of remote sensing techniques were used, such as the multi-temporal vegetation indices (NDVI, SVI, EVI) (del Castillo, García-Martin, Aladrén, & de Luis, 2015; Czerwinski, King, & Mitchell, 2014; Höpfner & Scherer, 2011; Mancino, Nolè, Ripullone, & Ferrara, 2014), the spectral band analysis (Röder, Udelhoven, Hill, del Barrio, & Tsiourlis, 2008; Sexton, Urban, Donohue, & Song, 2013) and conventional classification methods (Harrison, Rivard, & Sánchez-Azofeifa, 2018; Melville, Lucieer, & Aryal, 2018; Murray et al., 2018). In this context, remote sensing data from several satellite sensors (i.e., Landsat, ASTER, Sentinel-2A, MODIS), which is now available free and provides temporal datasets for few decades, has become increasingly used to identify the forest cover types and their relative changes.

Morocco, a semi-arid/arid country, benefits from a considerable floristic diversity represented by forest ecosystems occupying 31,482 km² of the total area (Lahlaoi, Rhinane, Hilali, Lahssini, & Moukrim, 2017). Forest stands vary due to climatic, orographic, lithological and anthropogenic criteria, but they are often linked to landforms. The forest provides important services to local and regional society such as global carbon storage, soil protection, hydrological and biochemical cycle, biodiversity, economy and especially local well-being. In the last five decades, these fragile ecosystems have begun to suffer greatly from anthropogenic impacts and climate change that remain, nevertheless, under-treated. Deforestation due to the intense use of woody resources by local rural societies, over-pasturage, ageing of forests and hydric stress due to climate change, have an impact on the tree regeneration and the forest production (Bahri, Haboudane, Bannari, Bonn, & Chillasse, 2007; Haboudane & Bahri, 2007; Hammi et al., 2010; HCWFCD, 2017; Rejdali, 2004). Despite efforts made by the Moroccan High Commissariat for Water, Forestry and Combating Desertification (HCWFCD), particularly through national forest programmes (PNF) promoting forest conservation and development of the forest sector, the reforestation and natural regeneration of forest resources remain insufficient to reverse the observed degradation trend, due mainly due to anthropogenic and arid constraints.

To assess their status, more studies have been conducted, but they mainly concerned monitoring the threatened or vulnerable vegetal species, namely regression of Argan trees (Occidental High-Atlas and Souss plain) (de Waroux & Lambin, 2012; Linares, Taïqui, Sangüesa-Barreda, Seco, & Camarero, 2013), Atlas cedar in the Middle Atlas forests (Bahri et al., 2007; Haboudane & Bahri, 2007; Linares et al., 2013), cork oak communities at Mâamora (Laakili et al., 2016; Saber, Rhazi, Rhazi, & Ballais, 2008), and Juniperus thurifera in the High-Atlas of Marrakech (Montès, Bertaudière-Montes, Badri, Zaoui, & Gauquelin, 2002), Carob trees in Azillal Province (Hajib, Lahlaoi, Alaoui, & Khattabi, 2015) and green oak in High Atlas (Slimani, El Aboudi, Rahimi, & Khalil, 2017). Many of them used an approach based on field data and visual image interpretation from different sensor data to detect vulnerable tree species change.

The eastern part of the province of Béni-Mellal, a study area, which is part of the Atlas mountain of Béni-Mellal, has a significant forest cover whose preservation of biodiversity requires sustainable management. However, it remains highly vulnerable to various pressures and aggression caused by multiple and complex interactions between natural and anthropogenic factors, and its monitoring is, therefore, necessary to check and estimate the change in forest cover and to propose appropriate conservation plan. Thereby, monitoring of this forest cover is, consequently, essential given their vulnerability to anthropogenic pressures and those associated with climatic fluctuations. While monitoring has often revolved around field data, remote sensing may play a fundamental role in the diagnosis of sustaining species and habitat diversity and the quantification of degradation or regeneration in the study region. The aim of this work was to map, classify and evaluate the spatial and temporal distribution of tree species and densities in the Eastern area of Béni-Mellal province using high spatial resolution and free satellite images data. ASTER and Sentinel-2A MIS images acquired in 2001 and 2015, respectively, constituted a principal source of data based on supervised classification and NDVI calculation. The field data, images of Google Earth and a local forest map were used to validate the accuracy of the derived forest stand maps.

2. Materials and methods

2.1. Study area

The study area corresponds to the Eastern mountain area of Béni-Mellal Province, and lies between 32°28′-32°24′N and 5°54′-5°52′W with a total area of 91819 ha (Figure 1). It is part of the Atlas chain of Béni-Mellal, and covers six mountain communes (El Ksiba, Aghbala, Naour, Tzi Inisly, Tanougha and Boutferda) that are underdeveloped and faces many socio-economic challenges resulting from the high level of poverty and marginality, which manifests in a high dependency on forest resources leading to overexploitation. Agriculture and livestock are the fundamental component of the rural economy throughout the region. In most cases, these are a traditional farming and extensive livestock herding. The total human population in 2014 was about 86869.

Figure 1. Location map of the study area. Source: Author.

The study area experiences sub-humid climate that is governed by seasonal changes. The average annual precipitation is between 22.7 and 1003.3 mm, with a mean of 438.7 mm. The wet season lasts from mid-October to May with the wettest month usually being November, while the dry season takes up the rest of the year. The bedrock consists mostly of limestone and dolomites, covered with soils that are generally fersialitic and brown forest types. There are other less well-represented soils such as poorly developed soils occurred on steep slopes and brown calcareous soils developed on soft limestones or on marl-limestone.

The altitude ranges from 869 to 5542  m, with an average of 1778  m. The forest ecosystem has much changed due to the land use change (such as conversion of forests and grasslands into cropland) and climate change in the past decades. Deforestation is prominent, and the demand for forest land for agriculture and pasturage is so considerable. Therefore, a detailed spatial and temporal analysis of the forest resources is critical and urgent for better forest resource management in the eastern areas of Béni-Mellal province.

2.2. Materials

In this study, multispectral data from two different satellite sensors, ASTER (Reflective Radiometer and Space Advanced Thermal Emission) and Sentinel-2A MIS, were used to map species distribution and forest density in 2001 and 2015, respectively. All images are acquired during sunny periods, without cloud cover. ASTER images were obtained from the U.S. Geological Survey (USGS) (http://earthexplorer.usgs.gov//), and those of Sentinel-2A MIS were collected from the official Copernicus Sentinel-2A MIS Data Hub (ESA 2016) (https://scihub.copernicus.eu/s2/#/home). These images were downloaded with a preprocessing level (geometric correction in UTM WGS84 projection), at the same time of the year to avoid problems of seasonal variation. Near-infrared and visible bands have been used knowing that many studies carried out in the world show that spectral information from satellites can be used to identify vegetation cover (Guyot, Baret, & Major, 1988; Warner, Skowronski, & Gallagher, 2017). The characteristics of satellite images and bands used in this study are summarized in Table 1.

Table 1. Satellite images and bands used in the analysis.

Image	Date d’acquisition	Resolution	Bands (μm)
ASTER	06 December 2001	15 m	Band 1: Green (0.52–0.60)
	18 December 2001		Band 2: Red (0.63–0.69)
			Band 3: Near infrared (0.78–0.86)
Sentinel-2A	16 December 2015	10 m	Band 3: Green (0.5425–0.5775)
	26 December 2015		Band 4: Red (0.65–0.68)
			Band 8: Near infrared (0.7845–0.8995)

All image processing, classification and GIS analyses were performed using QGIS 2.8.6, ArcGIS 10.2.2 and MATLAB R2013a softwares. The ASTER digital elevation model (DEM) with a spatial resolution of 30 m was used to extract topographic features of the study area, such as elevation and aspect.

Dominant higher plant species are stratified into four stands, with respect to classification map available from HCWFCD. This 1/20000 scale forest map was produced by the National Geographic Institute (IGN) (1962). The four stands are Holm oak (Quercus rotundifolia), Aleppo pine (Pinus halepensis), Thya (Tetraclinis articulata), Zea oak (Quercus canariensis). Some information on the height and diameter of the species is determined from our field observations. Holm oak is the backdrop for all forest stands in the study area, occurring at an elevation between 100 and 1819 m. Amongst the stands of oak, there are coppices of various sizes, and stands of young to adults, with different levels of density. Aleppo pine occupies a small extent, and occurs in the form of adult stands, generally located on the south and southeast slopes. The height of the trees reaches 12 m, and their diameter oscillates between 15 and 35 cm. Thya formation constitutes the main resinous species of the studied forest. It is found at the mean elevations of the southern and southeast mountainous slopes. The height of the Thya varied between 1 (the youngest trees) and 5 m. Zea oak formation covers a very small area of the study area and occurs as a moderately dense adult forest structure. It occupies the north and northwest mountains slopes with elevations between 1550 and 1770 m. Some Zea oak trees reach a diameter greater than 60 cm and a height over than 13 m. The cover land was also comprised of bare ground and crops & others (lower plants and shrub). Crops & others correspond to areas cultivated with annual crops, vegetables, or fruit, lower wild plants and shrubs. Most of the cultivated area has been recently managed for modern fruit cultivation, like apple or cherry. Bare ground constitutes to the exposed soils, urban and rural areas and roads.

A field campaign was targeted to ascertain the accuracy of the classification by field observations, and the accuracy of the land cover map generated. It should be noted that two field missions were carried out, one of reconnaissance in April 2016 and the other of validation of the results in May 2016. Forty-three control points showing complicated vegetation spatial patterns or homogeneous plant stands were established using Global Positioning System (GPS) instrument (Figure 1). Validation and quantitative analysis of accuracy were performed using confusion matrices, global precision, and Kappa coefficient. The forest map available from HCWFCD served as a reference to complete the validation by superimposing it and comparing it with the satellite images.

2.3. Methods

The methodological framework adopted in this study for monitoring changes in canopy extent and density of forest stands is presented in Figure 2. The spatial images of ASTER (2001) and Sentinel-2A MIS (2015) were selected free archive availability. These images were processed using QGIS 2.8.6 to correct atmospheric errors and resample the resolution of the ASTER image from 15 to 10 m. Atmospheric correction is made using the Semi-Automatic Classification Plugin extension (Congedo, 2016) based on the Dark Object Subtraction (DOS) model. Then the images mosaic method was applied because of the study area that is located between two successive images for the two satellite products.

Figure 2. Methodological flowchart of the procedure used in the study.

2.3.1. Image processing

Before the application of the classification method on processed images, the coloured compositions of each product were tested to select the best combination of bands. The best band combination is the coloured compositions of 3-2-1 (RGB = NIR-Red-Green) for the image ASTER and Sentinel-2A MIS (Desclée, Bogaert, & Defourny, 2006). The choice of the bands in the infrared (band 3 for ASTER and band 8 for Sentinel-2A) and in the red (band 2 for ASTER and band 4 for Sentinel-2A) allowed to increase the contrast and to discriminate the vegetation formations between them and to differentiate them from other land-use classes. Then we applied a classification by the maximum likelihood algorithm (Chutia, Bhattacharyya, Sarma, & Raju, 2017; Kavzoglu, 2017) on the different images, which assumes a normal distribution of spectral pattern for each class (Xie, Sha, & Yu, 2008). This classification method requires the selection of forest formations in the study area. For this reason, we used a regional forest map was used identify the forest cover types, and the field data were used to support the interpretation process by assigning each pixel type to a specific land cover class. Six different land cover types, namely Holm oak, Aleppo pine, Thya, Zea oak, crops & others, and bare land, were used in the classification process. In this classification, the pixel categorization process is done by specifying, a numerical description of land cover classes previously defined. To check if the training samples were different radiometric characteristics, the signature separability was used since the spectral signatures of similar features have similar shapes.

Since the classified maps often contain some error, accuracy assessment is an important step in the classification process. Accordingly, an assessment to ascertain the accuracy of our classification portrayed forest cover was conducted based on field data, local forest map used in training step, and Google Earth image. The overall accuracy, kappa statistics and user’s and producer’s accuracies (Congalton & Green, 2008) were determinate to evaluate the accuracy of the supervised classification used in this study for two dates (2001 and 2015) images. The Kappa values lie between 0 and 1, with Kappa values >0.80 represent a strong agreement and good accuracy, between 0.40 and 0.80 indicate a middle accuracy, and <0.40 indicate a poor agreement between classification and observation.

2.3.2. NDVI differencing

The method of thresholding NDVI (Slimani et al., 2017) was adopted to highlight the density of each forest stand. The NDVI considered as one of the common vegetation indices because of its known theoretical and empirical relationships with vegetation abundance and greenness (Goward, Markham, Dye, Dulaney, & Yang, 1991; Pettorelli et al., 2005; Wang, Adiku, Tenhunen, & Granier, 2005), was used to identify and quantify the vegetation. It was calculated from the red and near-infrared (NIR) band values, and its values range from -1 (unvegetated areas) to +1 (healthy vegetated areas). The NDVI extracted from our study area has values that range from −0.14 to +0.74. The two NDVI images generated in 2001 and 2015 were processed for extracting three forest cover density classes using image thresholding methods. According to the NDVI value, the image was classified into three distinct thematic classes: NDVI <0.1: without vegetation, 0.2 < NDVI < 0.4: slightly dense forest, 0.4 < NDVI < 0.74: dense forest. The identification of the best-fitting threshold value was based on the comparison between the Google Earth images and the NDVI results.

3. Results

3.1. Analysis of forest cover changes

3.1.1. Accuracy assessment

The overall accuracy, kappa statistics and user’s and producer’s accuracies (Congalton & Green, 2008) have been used to assess classification accuracy in the present study. The confusion matrix for 2001 and 2015 image’s classification, including the results of accuracy and kappa coefficient assessment are given, respectively, in Tables S1 and Table S2 in Supplementary material.

The overall accuracy of the supervised classification method was about 97.76% in 2001 and 95.80% in 2015, according to obtained overall Kappa values of 0.971 and 0.946, respectively. Knowing that the critical accuracy value of 75% beyond which a classification is deemed acceptable (Girard & Girard, 2010), these results indicated the best classification accuracy. In both dates, the bare ground, Aleppo pine, and Holm oak classes were the most discriminated according to the scores on the diagonals of Table S1 and Table S2 (Supplementary material). However, the great confusion occurred between the Thya, Aleppo pine and crops & others stands. This confusion results in close spectral responses. Regarding producer’s accuracy, all classes were above 85% in both years. The highest values are for the Aleppo pine stand followed by the crops & others class, and then by Holm oak in both years. The user’s accuracy was above 90% except for the Aleppo pine stand showing user’s accuracy of 83.33% in 2001. The strongest scores are also for the classes of Holm oak, bare ground, and Thya.

A spectral signature analysis between the identified land cover classes performed during the classification process to discriminate these classes, showed that it is impossible to separate between all species, namely the forest trees rarely represented and the shrub. For this, these species, shrubs, and trees seldom represented, were merged with crops into a single class termed “crops & others.”

3.1.2. Spatial distribution of forest stands

The results of the digital classification deemed reliable were used to assess the vegetation stands being changed over this period from 2001 to 2015. At the class level, the area per class of land cover and its expansion or retraction quantified for each map are represented in Figures 3 and 4, Table 2.

Figure 3. Forest stand map for 2001. Source: Author.

Figure 4. Forest stand map for 2015. Source: Author.

Table 2. Total area per class and temporal variation.

Land cover classes	Area in 2001		Area in 2015		Variation
	ha	%	ha	%	ha	%
Holm oak	44119.11	48.04	48311	52.61	4191.89	+4.57
Aleppo pine	2537.31	2.76	2638.51	2.87	101.2	+0.11
Thya	4837.00	5.26	6302.16	6.80	1465.16	+1.54
Zea oak	67.73	0.07	46.73	0.05	−21	−0.02
Crops & others	8135.91	8.86	7666.65	8.35	−469.26	−0.51
Bare land	32122.33	34.98	26854.37	29.24	−5267.96	−5.74
Total	91819.42	100	91819.42	100

Forest was the dominant land cover in 2001 with 56.15% of the entire study region of the eastern area of Beni-Mellal province, followed by the bare ground class (34.98%) and the crop & others class (8.86%). The forest area increased by about 6.25% of the whole area in 2015 (Table 2) or by about 5728.5 ha.

Focusing on the class level (Table 2), the crops & others class also decreased by 0.51% (from 8135.91 ha in 2001 to 7666.65 ha in 2015), and the bare ground class retracted by 5.74% (from 32122.33 ha in 2001 to 26854.37 ha in 2015). Within the forest area, Holm oak is the most expanded species occupying the large extension of about 48.04% of the total area in 2001 and 52.61% in 2015. Aleppo pine has established a small area of about 2537.31 ha and 2638.51 ha in 2001 and 2015, correspondingly. Thuya stand accounted for about 5.26% (4837 ha) and 17.1% (6302.2 ha) of the total study region, respectively, in 2001 and 2015. Zea oak species has the lowest area (67.73 ha) compared to other classes, which is about 0.07% in 2001 and 0.05% in 2015 of the study region. The bare ground area reduced from about 35% of the study region in 2001 to 29.24% in 2015.

From this analysis, an overall land cover change between 2001 and 2015 was identified with an increase in a forested area. All species stands showed expansion at the expense of the bare ground and crops & others lands.

3.1.3. Density change of forest stands

The density assessment maps for 2001 and 2015 were reclassified to three density classes separated into low (NDVI <0.1), moderate (0.2 < NDVI < 0.4) and high density (0.4 < NDVI < 0.74) (Figures. 5 and 6), using ArcGIS 10.2.2 software. The results of the spatial and temporal evolution of the density of the canopy in the study area during the 2001 and 2015 periods are reported in Table 3(a). The spatial pattern of classified forest cover density zones indicates that the areas with high canopy density are located mainly in the west of the study region, while the areas with low canopy density are located in the northwest, central east and southeast of the study region (Figures 5 and 6). Overall, the forest experienced a somewhat slight density change in the period analyzed (Table 3(a)). The dense forest increased from 24814.50 to 27160.87 ha between 2001 and 2015, signifying a rate of 0.02%. This increase in high-density areas occurred at the expense of low-density areas, and could likely to be related to tree development in these areas especially the Holm oak and Aleppo pine stands in 2015.

Figure 5. Forest density map for 2001. Source: Author.

Figure 6. Forest density map for 2015. Source: Author.

Table 3. Temporal variation in forest cover density (a) and density variation per forest stand (b).

Density		Area in 2001		Area in 2015		Variation
		ha	%	ha	%	ha	%
(a) Forest cover density
	Low dense	29818.70	0.33	27002.25	0.30	−2816.45	−4.96
	Dense	24814.50	0.27	27160.87	0.29	+2346.37	+4.51
(b) Forest stand density
Holm oak	Low dense	22271.33	24.25	21413.20	23.32	–858.13	+1.66
	Dense	21847.78	23.79	26897.80	29.29	+5050.02	+9.79
Aleppo pine	Low dense	258.24	0.28	102.93	0.11	−155.31	−0.30
	Dense	2279.07	2.48	2535.58	2.76	+256.51	+0.50
Thya	Low dense	1086.90	1.18	161.75	0.17	−925.15	−1.79
	Dense	3750.10	4.08	6140.41	6.68	+2390.31	+4.64
Zea oak	Low dense	12.15	0.013	13.09	0.01	+0.94	+0.01
	Dense	55.18	0.06	33.64	0.03	−21.54	−0.04

The final map with the density of each forest stand was also derived using the geometric intersection of the land cover density map and the spatial distribution map of forest stands by applying intersect geoprocessing function in ArcGIS. The overlay of the stand type map with the density map concerned only the Holm oak, Thya, Aleppo pine and Zea oak stands (Figures 7 and 8, Table 3(b)). The results indicate that the Holm oak species, moderately dense in 2001, became very dense in 2015, especially in the northwest of the study area, the north of Koumch Village and south of Boutferda Village. The Aleppo pine density increased south of the study region (Tizi n’Ait Ouirra) and west of the rural commune of Tizi-Nisly. The Zea oak formation density remained stable and dense. The Thya formation observed throughout the forest area except for southern part has undergone a significant density evolution in the period 2001–2015.

Figure 7. Forest stand density map for 2001. Source: Author.

Figure 8. Forest stand density map for 2015. Source: Author.

4. Discussion

Analysis of the dynamics of change revealed a slight increase in forest stands between 2001 and 2015. This increase in surface could mainly be explained by the effect of the enhanced protection by forestry authorities, the reforestation programs initiated by the High Commission of Water and Forests for the reconstruction of certain vegetation facies, namely pine, the climatic conditions that were more favourable during the period 2001-2015, and also the geographical and topographical position characterizing by a microclimate attenuating the effects of drought in dry seasons on the vegetation. The transformation of forest areas may be partly associated with the impact of a tree cutting promoting the stand regeneration. Nevertheless, this increase in surface and density of the forest tree stands remains relative as the rates of change remain low to moderate, and as the resolution of the ASTER and Sentinel-2A MIS products used could influence the accuracy of the spatial analysis. Such weak changes would also be linked to the high anthropogenic pressure that impedes the regeneration and dynamics of these ecosystems. The high anthropogenic pressure could also explain such weak changes in the region, which is hampering the regeneration and dynamics of these ecosystems.

The results of this study will be certainly a reference to the state of the forests and the evolution of the forest formations in the study area, which remain very vulnerable to anthropogenic activities as well as climate change. The dominant anthropogenic drivers impacting the forest cover in the study region, lie in the clandestine deforestation, overgrazing, expansion of agricultural land and urbanization. The causes of deforestation are mainly due to the habits of local populations and their extreme poverty, which leaves them too dependent on forest resources, and despite attempts to implement a management transfer system, the illegal exploitation remained considerable. Furthermore, the local population more accustomed to traditional methods in agriculture encroached the forest areas by eliminating stubble and shrubs to conquer the new cultivable land. One underlying cause of forest cover change is related to the excessive expanding pasturage. Despite these multiple pressures and aggression caused by interactions between natural and anthropogenic factors, the results showed that the study forest is doing quite well. This is, unfortunately, the case for many Moroccan regions suffering from intensive overgrazing; criminal harvesting of wood and parasitic diseases that are all factors that, combined with climatic uncertainties, cause abnormal evolution and regeneration in some parts of these forests (i.e., Bahri et al., 2007; de Waroux & Lambin, 2012; Haboudane & Bahri, 2007; Hajib et al., 2015; Hammi et al., 2010; Laakili et al., 2016; Linares et al., 2013; Rejdali, 2004; Saber et al., 2008). This forest state of the study area that appears to be normal would be due to favourable climatic conditions during the period 2001–2015, above all, to the ambitious approach followed by High Commission for Water, Forests and the Fight against Desertification – HCWFFD, where, in addition to protection, reforestation is the central action (HCWFCD, 2017).

The results also show that using and combining spectral signatures and NDVI index lead to interesting results on mapping land cover stands and estimating their densities. The results obtained constitute the first remote sensing-based product for mapping forest cover in the region. The present study demonstrates how satellite-based detection of vegetation change can provide reliable results in the assessment of natural forest dynamics. The study confirms thus how remote sensing can contribute to improved management of forest reserves, and highlights initiatives and challenges in this domain. Particularly since access to remote sensing data with better quality becomes free such as the Sentinel 1A and 2A satellites that provide high-resolution radar and optical imagery on the Internet.

5. Conclusion

This study explored the potential of remote sensing for forest stands mapping and density in the eastern part of Beni-Mellal province, using ASTER and Sentinel-2A MIS images acquired in 2001 and 2015 respectively. The principle of the methodology was the use of supervised classification to map the different types of forest stands. Also, the NDVI vegetation index was calculated to estimate the density of the forest. The field data, Google Earth images, and a local forest map were used to validate the accuracy of the derived forest stand maps. This assessment constitutes the primary study of the forest in the region using remote sensing tools and provides the first forest maps offering more information about a spatial distribution of forest stands in the study area. The resulting maps showed expansion in extent, and density of forest stands on bare ground class resulting from natural processes and management actions adopted by HCWFCD. Forest area and density increased to about 6.25% and 4.51% in 2015, respectively. Holm oak stand dynamic appeared to be most important with an expansion of 4.57% and a density of 9.79%. Aleppo pine and Thya stands also expanded by about 0.11% and 1.54%, respectively, in the examined period, and their densities showed significant growth in 2015. Zea oak, Crops & others and Bare ground classes displayed a net retraction.

The results of this study showed the potential of the remote sensing products to assess the large forest area. Long-term monitoring of forest cover becomes possible to forest managers using free multi-date satellite images due to their availability and the relative simplicity of the methodology. Nevertheless, to discriminate vegetable species using multispectral images alone remains a challenge but is possible with additional data and expert information.

This study allowed highlighting more information about forest stands and their dynamics that can be considered a useful starting point to further spatial and temporal analyses of forest cover change in other regional and national areas.

Disclosure statement

No potential conflict of interest was reported by the authors.

Funding

This work was supported by Sultan My Slimane University.

Supplemental data

The supplemental data for this article is available online at https://doi.org/10.1080/24749508.2018.1452478.

Acknowledgements

This work was conducted in collaboration with the Provincial Direction of Water, Forests and Desertification Control, Béni-Mellal, Morocco. We are grateful for their support, especially Mr Abdelfettah Naouri.