Monitoring Land Use and Land Cover Changes Using Remote Sensing Techniques and the Precipitation-Vegetation Indexes in Morocco ECOLOGICAL ENGINEERING & ENVIRONMENTAL TECHNOLOGY

The study of land use and land cover change (LULC) is essential for the development of strategies, monitoring and control of the ecosystem. The present study aims to describe the dynamics of land cover and land use, and specially the impact of certain climatic parameters on the distribution of vegetation and land cover. For this study, multi-temporal remote sensing data are used to monitor land cover changes in Morocco, using a set of Landsat images, including Landsat 7 (ETM+), Landsat 5 (TM), and Landsat 8 (OLI), captured during the period 2000–2020, those changes were determined by adopting the maximum likelihood (ML) classification method. The classification results show good accuracy values in the range of 90% (2000), 80% (2007), 82% (2010), 93% (2020). The LU/LC change detection showed a decrease of agricultural and forest areas in the order of 5% between the year 2000 and 2020, and an increase of bare soil of 5% to 6%, and a notable change in urban area from 97.31 ha (0.03%) in 2000 to 2988.2637 ha (0.82%) in 2020. The overall results obtained from LULC show that the vegetation cover of the study area has undergone major changes during the study period. In order to monitor the vegetation status, an analysis of the precipitation-vegetation interaction is essential. The normalized difference vegetation index (NDVI) was determined from 2000 to 2020, to identify vegetation categories and quantify the vegetation density in the Lakhdar sub-basin. The obtained NDVI was analyzed using climatic index SPI (Normalized Precipitation Index) based on rainfall data from five stations. The correlation study between NDVI and SPI indices shows a strong linear relation between these two indicators especially while using an annual index SPI12 however, the use of NDVI index based on remote sensing provides a significant result while assessing vegetation. The results of our study can be used for vegetation monitoring and sustainable management of the area, since it is one of the largest basins in the country.


INTRODUCTION
The use of satellite images and data to assess and analyze changes in land use and land cover Congedo, 2020;Chavez, 1988], has become the most recognized and powerful technique to obtain more accurate information on land surface characteristics at different temporal and spatial scales. These images also allow the analysis of vegetation conditions and their response to climate change. In Morocco, several researches have been conducted to study land use using remote sensing data.  used multispectral AS-TER (Advanced Spaceborne Thermal Emission and Reflection Radiometer) and Sentinel-2A data acquired in 2001 and 2015, to monitor forest cover dynamics in the eastern area of Beni-Mellal province, the study was based on the supervised classification algorithm and NDVI combined in a GIS environment to quantify the extent of change in the density of forest stands, the result revealed an overall change in forest cover with an increase in the wooded area. [Gurgel et al., 2017] identified forest changes in argan forests (Morocco) using aerial photographs and satellite images between 1970 and 2007. Their study revealed a decrease in forest density of 44.5% during this period.
For our case study, Landsat5-TM, Landsat7-ETM+, Landsat8-OLI satellite images are used to Monitoring Land Use and Land Cover Changes Using Remote Sensing Techniques and the Precipitation-Vegetation Indexes in Morocco Fatiha Ait El Haj 1* , Latifa Ouadif 2 , Ahmed Akhssas 2 1 L3GIE Laboratory, Mohammadia Engineering School, Mohammed V University in Rabat, Morocco identify changes and map land cover, it consists of a set of images at different dates (from the year 2000 to 2020) taken at the same period of the year, to capture the conditions of the study area and improve the classification of changes over different periods.
In order to classify and analyze the changes LULC in the Morocco sub-basin, we performed the maximum likelihood method. This is the most efficient algorithm for supervised classification [Bonn & Rochon, 1993;Chatelain, 1996], which seeks to build a function called the likelihood function and to maximize its logarithm considering the unknown parameters. The coherence of the maps obtained was evaluated using the value of the Kappa coefficient. To measure and research how vegetation responds to climate change, meteorologists and climatologists from all over the world have developed a variety of drought indicators. Indicators that are simple to calculate and statistically significant. McKee, Doesken, and Kleist, three American scientists, developed the Standardized Precipitation Index (SPI) in 1993 [Mckee et al. 1995]. SPI is a powerful, flexible, and easy-tocalculate index, precipitation is actually the only required input parameter. Furthermore, it works just as well for evaluating dry periods and cycles as it does for assessing wet ones. Landsat satellites have particular characteristics in terms of spatial, spectral and temporal resolution. The choice of this type of satellite favors the study and analysis of surface reflectance. Among the monitoring indices of vegetation dynamics that is based on the monitoring of surface reflectance is the normalized difference vegetation index (NDVI) used to better characterize the spatial extent of drought events and to monitor the vegetation status [Congedo, 2020; Abebe et al. 2022;Jiang et al. 2022]. In summary, the objective of this study is to explore the changes of LU/LC of the sub-basin in Morocco and to monitor the vegetation status that is affected by precipitation inputs and drought status in the area. The results of this study could be useful for vegetation monitoring and sustainable management of the area as it is one of the largest basins in Morocco.

Study area
The sub-basin of Lakhdar [AHT Group, 2016] belongs to the hydraulic system of wadi Oum Er-Rbia Which belongs to the mountain area of the province of Beni-Mellal and has an area of 3503 km 2 . The Lakhdar sub-basin is located at the eastern end of the Haouz Mejjate basin, bounded to the south and south-east by the High Atlas mountains and to the west by the Tassaout sub-basin ( Figure 1). The topography [AHT Group, 2016] of Lakhdar sub-basin is relatively rugged, and characterized by the variation in altitude between 488 m at the Tensift wadi, to 4017 m at the High Atlas.The Lakhdar sub-basin is drained [AHT Group, 2016] by the main watercourse (Wadi Lakhdar) which drains the western and southwestern parts of the basin, its main tributary is Wadi Bernat which drains the northern and northeastern parts. It includes the Hassan I dam (built in 1986) and the Sidi Driss dam (built in 1984), which ensure the irrigation of agricultural perimeters and the supply of drinking water to the population of the area. The geological formations of the region [AHT Group, 2016] are composed mainly of conglomerates in a band of 2-3 km high aligned at the foot of the Atlas that extends northward in the axis of the current course of the wadi Lakhdar, and alluvial formations reworked in the quaternary, consisting of pebbles, gravels and sands with high permeability corresponding to former wadi beds, and permeable formations pass by going north to silty formations, sometimes encrusted surface. The rest of the northern part of the sub-basin Lakhdar is constituted of an alternation of Jurassic (Lias, Jurassic and Jurassic red sandstone of the High Atlas). In the southern part, formations of the older secondary (Triassic) and primary outcrop, including the unsubdivided Devonian and the unsubdivided Ordovician.
The studied area is categorized by an arid climate [AHT Group, 2016] with temperate winter covering the entire plain area of the sub-basin (18% of the sub-basin in terms of area), a semiarid stage with cool winter occupying the piedmont area (12% of the sub-basin), humid with temperate winter covers only 3% of the area of the sub-basin, and humid with cool winter which covers 46% of the area of the sub-basin and most of the mountain area.

Materials
The majority of the data used for this study was gathered from two sources: satellites sources, and rainfall data collected from five ground stations. Satellite data -multispectral images from Landsat satellite sensor, and depending on the date search and image quality used to map land use from the year 2000 to 2020, several types were used: Landsat 7 -Enhanced Thematic Mapper Plus (ETM+), Landsat 5 -Thematic Mapper (TM), Landsat 8 -Operational Land Imager (OLI) sensor images, which were obtained from the U.S. Geological Survey (USGS) (http://earthexplorer. usgs.gov//). All images were acquired during sunny periods, without cloud cover, and were pre-processed (geometric correction with a UTM WGS84 projection and atmospheric correction) in order to compare changes in land cover and land use, all images were acquired at the same period between June and September, which is the period with the stable vegetation stage.
The characteristics of the satellite images and resolution used in this study are summarized in Table 1. Image processing, classification and analysis are performed in ArcGIS 10.3 and ENVI 5.3 software. The processed images have a resolution of about 30 m with multi-spectral bands which facilitates the extraction and characterization of both vegetation cover and bare ground in the study area. The vegetation cover extracted according to these images and according to the field survey, is divided into two main types, forests and cultivated land (arboriculture of different sizes, cereals, olive trees...), and for bare soil it is mainly rocky outcrops bare spaces and / or eroded, and exposed soil (urban areas, rural roads, tracks...). In order to verify and analyze the accuracy of the established maps, precision studies were conducted by determining the confusion matrices, the overly precision and the

Meteorological data
The annual and monthly rainfall data made by the hydraulic agency of Oum er-rbia ABHO for the period of 2000 to 2020 were used. The unit of rainfall is the millimeter (mm). The data are organized in a grid. Figure 2 shows the location of the weather stations (5 stations) in the Lakhdar sub-basin district.

Methods
The methodology adopted to study land cover at the different dates and to identify vegetation changes is summarized in Figure 3. The Landsat images used are mainly chosen according to their availability for the study area [Faruque et al. 2022;]. The downloaded images followed radiometric corrections (Radiometric calibration extension) and atmospheric (semiautomatic extension by the Dark object subtraction model DOS) on ENVI 5.3 software. for the extraction of the study area according to the basin boundaries we used the image mosaic method on Arcgis in order to identify the limits of the area of study.

Image pre-processing
The processing of the images was conducted under ENVI 5.3 software. The images had initialy undergone all geometric and radiometric where: L λ -spectral radiance at the sensor aperture for a single band mW/(m 2 ·sr·µm), L λmax and L λmin -scaled spectral radiance (provided in a header file of image in- The images are also influenced by the effects of the atmosphere, hence the need for atmospheric corrections to facilitate the extraction of information from the signal independently of the atmospheric effects that are variable in time and space. The atmospheric correction has been applied by the semi-automatic classification plugin [Chowdhury et al. 2020] based on the Dark Object Subtraction (DOS) algorithm [Congedo, 2020]. 2022] is to delimit at the level of the Raster image a set of polygons that will be assigned to a given class of LULC. The method is based on the classification of pixels surrounded by these polygons to create homogeneous spectral signatures, which will present the types of land cover recorded at the image, each class is differentiated from the other according to their characteristics and their spectral response. The training sites obtained will be validated from the visual interpretation of Google Earth images of the study area with the available field surveys.
Based on the field data and using Google Earth high resolution images as reference data, five categories of LULC were identified ( Table  2): (1) water body, (2) agriculture, (3) forest, (4) bareland, (5) urban area. These categories were identified and validated using GPS control points and reference data samples to generate training areas of each category.

Accuracy assessment
Accuracy assessment is a very important step in checking the accuracy of the classification re- The accuracy study is based on reference data that present the reality of the land cover, and the error evaluation is represented by a confusion or error matrix generated by ENVI 5.3 software. This matrix is in the form of a square table that gathers the pixels of the image where the row represents the classified category [Anand 2012;Jaouda et al. 2018;Purwanto et al. 2016] and the column represents the reality of the field, and the non-diagonal row indicates the values that are misclassified or that are classified in another category. The matrix is also used to calculate the global accuracy,

NDVI index
Based on Landsat images downloaded in Google Earth and according to the characteristics of the vegetation in the study area and soil type, the NDVI index was mapped to quantify the changes in vegetation over the period (2000-2020).
The NDVI is determined from information on the quantity and density of vegetation by con- where: NIR -represents the reflectance in the near infrared, RED -represents the reflectance in the red band.
NIR and RED represent the radiated reflectance in the near-infrared band 7 (0.8-1.1 μm) and visible red band 5 (0.6-0.7 μm) of the MSS Landsat 5 satellite, respectively. For the 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-vegetation cover.
Reference data was obtained on the area from field surveys and maps published by the hydraulic agency of the study area, and according to extensive discussions with stakeholders, three density classes were determined, namely: The NDVI class with values below 0, presents areas without vegetation cover and according to the reference data they were considered as water bodies and bare soil; the positive values between 0.15 and 0.3 were considered as medium density vegetation (cultivated area); and the positive values above 0.3 were considered as high density vegetation (forest).  [Bordi et al. 2011] the SPI method is mainly used because it can give a reliable comparison and is relatively easy to use in various climatic conditions and locations. The rainfall data was provided by ABHO which covers a period of 20 years (2000 to 2020), this data is used to calculate the SPI in R studio software. The SPI is used to measure drought for time scales ranging from 3 months to 28 months. The authors of this index define a drought event when the SPI becomes less than or equal to -1.0 ( Table 3 where: P -total rainfall of a period (mm), Pm -historical average rainfall of the period (mm), σp -historical standard deviation of rainfall for the period (mm).

SPI is an index developed in 1993 by Mc
This index defines the severity of drought in different classes (Table 3) (Table 4) [OMM, 2012]., these indices offer temporal flexibility in assessing rainfall conditions in relation to water supply.

Land uses/land cover changes monitoring LULC
The LULC classification map were conducted using Landsat multi-date images from the years If SPI 9 < -1.5 then this is a good indication that substantial impacts may occur in agriculture (as well as the possibility of other sectors).

SPI12 months Forms of long-term precipitation
Possible flows related to reservoir levels, and also groundwater levels.  2000 to 2020 by the ML maximum likelihood method, as shown in Figure 4. The results of this classification (Table 5 and Figure 4)

Accuracy assessment
The accuracy assessment using the confusion matrix in the post classification phase. The results of this classification give values of Overly Accuracy of the order: 90%, 79%,82%,93% for the years 2000, 2007, 2010, 2020 respectively and a kappa coefficient of 0.8473, 0.7091, 0.7492, 0.8814 respectively ( Table 6). The Kappa coefficient [Anand, 2012;Baral et al. 2011;Xie & Ren, 2011] is a statistical measurement technique of the gap between the reference value (real image) and the identified classification value (classified image), the acceptable values which present a good correspondence between the classified pixels and the reference pixels set between 0.61 and 1 [Anand, 2012;Baral et al., 2011], the blue cells of the matrix (diagonal) present the acceptable classified pixels between the reference image and the classified image. The user and producer accuracies of the classes for the period under study are presented in Table 6, and the results obtained indicate that almost all the classes present user and producer accuracy values higher than 80% proving a good accuracy of the classification method.

Spatial distribution of NDVI over 2000 to 2020
The NDVI measures the balance between the energy received and emitted by ground objects and thus characterizes the vegetation mass present in a given environment [Rouse et al., 1973]. In principle, positive values correspond to denser vegetation and negative values to bare ground, clouds, lakes and rivers. In this study the NDVI is determined based on satellite images taken in the same period of the year (June & July). The values obtained vary from -0.55 to 0.75 over a period of 20 years ( Figure 5). According to the results obtained ( Figure 5), each map shows variations in NDVI values for the study area. NDVI varies from a positive value of 0 (very low to no photosynthetic activity) to 0.77 (very heavy and dense vegetation). The values between 0.15 and -0.28 are often assimilated to bare soil, while the very negative values below -0.28 present other forms of land use, they are attributed in this case and according to field surveys to water bodies or snow. In terms of variation over time, the state of vegetation has undergone very significant changes from one year to another, where we can see from Figure 5

SPI Index
Being standardized, the SPI has the advantage of comparing drought conditions over various time periods. The values of SPI are calculated using R studio software. Comparison of the 3-months, 6-months and 9-months seasonal SPI indices highlights the seasonal characteristics of rainfall and drought conditions.

DISCUSSION LULC changes from 2000 to 2020
The analysis of land change (Table 5 and Figure 5) shows an increase in the urban area from 0.03% to 0.83% between 2000 and 2020 caused by the population density increase and a remarkable decrease in vegetation, mainly explained by the climate change impact causing frequent dry seasons. Among the strategies adopted for this area by the Hydraulic Agency of Oum-Erbia, building several dams to collect surface water aiming to improve the sustainability of these resources. That improvement is clearly noticeable in our mapping according to Table 3 and Figure  5, a significant increase in water resources from 0.63% in 2000 to 1.03% in 2010, which generate an increase in the vegetation area that can be seen in Table 3 (jumping from 10.25% in 2007 to 11.14% in 2010). Unfortunately, from 2010 other constraints cause the considerable decrease in these resources such as the accentuated impact of dry seasons and the important decrease of precipitation which result into a decrease of water cover from 1.03% in 2010 to 0.74% in 2020.The overall accuracy and Kappa coefficients for the LULC maps from 2000 to 2020 shows a good accuracy, which proves that the Maximum Likelihood Supervised Classification method used in this study was very effective in improving the land use classification. The conditional Kappa statistics of each LULC class were all between 0.70 to 0.80, except for barren land. Barren land can be confused with urban areas where the accuracy of the producer and user (Table 4) does not exceed 60%, this limitation of the obtained accuracy values is mainly due to the rather similar reflectance effect between the two categories.
Similarly Chooi et al. [2010], studied the urbanization and the resulting land use change by analyzing Landsat satellite images during the period from 1999 to 2007 in Penang Island, Malaysia, they chosed the maximum likelihood classification method, their study reveal that the urban area increased dramatically, and the grassland area increased moderately. Conversely, barren land decreased obviously, and forest area decreased moderately. They found according to their study that the maximum likelihood classification produced superior results and achieved a high degree of accuracy and that The remote sensing technique used in their study proved to be effective; it reduced the analysis time of urban expansion and proved to be a useful tool for assessing the impact of urbanization on LST. Ozesmi & Bauer, [2002] confirmed that the MLC leads to higher accuracy than the decision tree (DT) classification. In a related investigation at the same profile region, remote sensing geospatial technologies have been used with great effectiveness. For instance, based on the supervised classification algorithm and the normalized difference vegetation index NDVI, Barakat et al. [2018] used Sentinel-2A MSI images and ASTER (Advanced

Vegetation monitoring NDVI and SPI indexes
The revelation of a relationship between NDVI and SPI was tested by calculating correlation coefficients "r" according to a multiple regression model and the coefficient of determination "R²" between the different variables. The scattering of points on the SPI3/NDVI and SPI6/NDVI scatterplot illustrates a weak correlation between these two indices (Figure 7). The coefficients of determination R² between NDVI and SPI vary according to the time scale considered, which means that NDVI is not significantly correlated with SPI3 and SPI6 which show seasonal variations. On the other hand, we can notice a slightly significant correlation with SPI9, which has an R 2 of about 0.489, and it is more correlated with SPI12 where we have an R 2 of 0.61 (Figure 7).
Similar studies have been conducted, by Lei & Peters [2004] who confirms that the response of vegetation to precipitation is with a time lag, and the impact of water deficits on vegetation is cumulative. Rokhmatullah et al. [2019] studied the problem of dryness that occurs in several regions in Indonesia (Bekasi), where they monitored this phenomenon by using SPI and NDVI index to find the relationship between these two, the result of this research shows that there is a close relationship between a region with high NDVI value and SPI value. But in a region with a low NDVI value, several regions have high SPI values.

CONCLUSIONS
The objective of our study was to determine the LULC status over the entire period from 2000 to 2020, to characterize the vegetation cover that will influence the land erosion at Lakhdar sub-basin in Morocco. Thanks to processing and analyzing remote sensing and satellite images under computer tools (ENVI software and GIS), and field surveys a detailed study was conducted on the occupation and use of land in the area. The use of multi-date images and different satellites and sensors allowed to update the LULC over a period of 20 years.
During the 20-year period, several changes are observed in the Lakhdar sub-basin. The results of the supervised classification by the method of maximum likelihood, can be summarized as follows: a significant decrease in the vegetation cover of (5%), and a considerable increase in the bare land and urban environment during this period. This classification accuracy is checked using a confusion matrix to determine the margin of error between the classified images and the reference images with control points based on the field survey.
The accuracy values obtained show high values for all images where the value of OA shows a good classification and presentation of reality. The obtained results of accuracy using the value of OA and Kappa coefficient, prove that the supervised classification by the method of maximum likelihood is an effective method to identify the changes of LULC and it is also a strong tool in decision making to have a better management and preservation of this basin. The results of correlation between the point index based on climatic data (SPI) and the spatial index based on satellite images (NDVI), gave limited correlations with seasonal variations (SPI3 and SPI6), and a more significant correlation with the annual and average variation (SPI9 and SPI12). To characterize and monitor the state of vegetation, the SPI remains very limited in these applications, as it is based on defined locations data from weather stations and therefore does not present the spatial impact on the entire area, and it is an index that is calculated based on only rainfall data, and it neglects the influence of other phenomenon on the water status of the area, such as evaporation and temperature. On the other hand, the NDVI gives a better modeling of vegetation status and with more accuracy on the areas affected by drought. It is an index based on the use of satellite images, which allows a better mapping of daily luminance levels and it also clarified the inter and intra-annual variations of the greens which can be influenced by climatic conditions. The advantage of this approach is providing an impact analysis of the different spatial and climatic components on the state of vegetation, moving from the location scale SPI to the regional scale NDVI while keeping the same scientific rigor.