Monitoring of forty years of agricultural expansion in the Oum Er Rbia valley (Morocco). The use of Google Earth Engine com- pared to Sentinel Application Platform

Agriculture in Morocco has been extensive until the middle of the 20th century due to the distribution of rainfall and the availability of water. In the middle of the last century hydraulic works were built that allowed the transition to intensive agriculture by the increase of irrigated areas, allowing that in the territories where there is water for irrigation and the climate allows it, the crops adapt to the demands of the market. The objective of the study is to assess by satellite images the land cover between 1985 and 2020, analyzing the changes in cultivation areas, as well as the changes in desert, sub-desert and forest areas of the Oum Er Rbia hydrological basin in Morocco. Landsat satellite images have been used since 1984 by the US government (Aerospace and Geological Agencies). A series of vegetation indices (NDVI, RVI, TNDVI and EVI) have been used; among which TNDVI (Transformed Normalized Vegetation Index) stands out for its better accuracy, which has allowed us to distinguish vegetation in cultivated and forest areas, as well as arid zones. In addition, the study has compared the use of two methodologies to calculate changes in the coverage of the Earth’s surface, has used local image processing from the Sentinel Application Platform tool and has also used the Google Earth Engine tool. The latter being the most optimal, although at the moment it has great limitations. In both methodologies and in the different indices it has been possible to observe during these 35 years as the cultivated area has increased (related to the availability of water by the construction of reservoirs and canals), how plant cover has improved in forest areas, and a range of variations in arid areas.


Introduction
Rainfed and irrigated agriculture in North Africa in general has been characterized by human presence throughout history. The low population density has shaped an agrarian or agricultural-livestock system, especially in modern and contemporary times. Agriculture has been extensive until the 20th century, when hydraulic works were carried out that allowed the irrigation of many areas and the transition from extensive agriculture to intensive agriculture [1]. The tribal organization of the population in the mountainous and sub-desert territory also contributed to the existence of collectivities that carried out agricultural and livestock development in local communities [2].
Water is an essential resource and in Morocco the management of this resource is even more vital, especially due to its semi-arid climate and limited water availability, in addition to the effects of climate change that are already being felt, since in the last three decades the country has suffered five periods of severe drought, exacerbating the problems of the region [3]. It is a country of contrasts between mountain ranges that are difficult to access and vast plains open to the Atlantic; between a relatively humid ocean front and large desert expanses to the east and south. The presence of the Atlas Mountains conditions the climatic variability of Northern Morocco. The territories located to the NW of it, open to the Atlantic, are characterized by a more humid climate, while those located to the SE, facing the Sahara, account for 60% of the national territory, have an arid climate and rainfall is less than 400 mm per year. The distribution of rainfall and the availability of water for irrigation determine agricultural use.
Because of its climatic characterization, the state considers agriculture as territorialized, indicating that there are regions where this activity is not possible and others where it is a major economic asset [4]. The territory of Morocco is divided into 16 regions, which in turn are divided into provinces and communes. In each province there is a ministerial delegation that manages the administrative aspects of agricultural activity. The next level of organization downwards are the Agricultural Valorization Offices, which are located in each of the areas with large irrigated territories and their nearby rainfed areas. In semiarid regions, with rainfall of less than 400 mm per year, agricultural development is carried out at the highest altitudes, to take advantage of the mild climate, while in areas with rainfall of more than 400 mm, the agricultural strategy aims at intensification. In areas where irrigation water is available, crops are adapted to market demands. In mountain areas, they are dedicated to pasture and forest use, with priority given to environmental conservation. Finally, in desert and oasis areas, the uses are unique and very site-specific [4].
In 1966, the U.S. Geological Survey proposed the creation of Earth Resources Observation Satellites, but it was not until 1972 that the first satellite, called Earth Resources Technology Satellite (ERTS-1), was in orbit. Later, in 1975, it was renamed Landsat at the time of the launch of the Landsat-2 satellite. From that time to the present, the Landsat mission had continuity, with Landsat-9 currently in its manufacturing [5]. This has provided a database of images of the planet unimaginable in its beginnings, whose continuity over time allows the elaboration of time series studies such as the one presented here.
The information obtained by the satellite from the air is related to the reflectivity that the planet emits in the upper atmosphere, known as TOA. The information is collected in different bands whose characteristics initially covered the visible spectrum in the three primary colors of light (blue, green and red) known as the RGB triad and in the near infrared (NIR). With the improvement of sensors in both quality and definition, the measurement bands have increased, as has the size of the measurement pixel, which is getting smaller and smaller and with greater spatial definition.
The treatment of the different bands could be done according to the researcher's criteria, producing the initial works about vegetation [6][7][8] and a few years later about water bodies [9,10]. When a sufficient database existed, different treatments of the bands were created to be used for the simplest uses, such as studies of the presence of water bodies, vegetation studies and studies on the earth's crust and its different coverages (ice, dry land, vegetation, water, urban areas, temperature, etc.), such as those of Dowdeswell et al. [11] on the distribution of ice masses using Landsat and SPOT images; the work of Toll [12] on land use using Landsat or that of Otterman and Tucker [13] on the relationship between temperature, albedo and vegetation also using Landsat images.
In our case, we have explored the possibilities of the already established vegetation indices to describe the variations of agricultural areas, considering the time period, the territory under study and the computer tools to process the available information. The need to compare Landsat-5 and Landsat-8 images has considered using indexes based on the satellite's green, red and infrared bands. On the other hand, the information available in Google Earth Engine at the time of the study (summer 2021) has also been used to compare the results obtained with this new system with those obtained by local processing of satellite images downloaded from existing distribution servers.
The most used current vegetation indices are known by the abbreviation of their name indicating the procedure used to calculate them. The following list presents those that use Landsat visible and infrared bands for calculation: Normalized Difference Vegetation Index NDVI [14] Ratio  [19] Our hypothesis and objectives were that that the increase of the agricultural area in certain provinces of Morocco has been due to the availability of regulated water and the construction of hydraulic works for its conduction to areas to increase agricultural production. For this purpose, the Oum Er Rbia river valley has been considered as a case study, evaluating by means of several vegetation indices applied to satellite images, how the agricultural area has been changing from the 1980s to the present. The contribution of this study lies in comparing for the first time the use of the Google Earth Engine tool with local processing by SNAP and whether there are differences in terms of the results obtained, presenting the advantages and disadvantages of each system.

Study site
The Oum Er-Rbia river is located in Morocco (Figure 1), in the north of the African continent. It is one of the most important rivers in the country, its source is located in the Middle Atlas Mountains, 40 km from Khenifra (1800 m a.s.l.), in the province of Meknès Tafilalet, the middle section runs through the province of Tadla Azilal and the riverbed is the geographical boundary of the Grand Casablanca region to the north on the right bank and Doukkala Abda to the south on the left bank. After approximately 550 km, it flows into the Atlantic Ocean in the town of Azemmour, in the province of Doukkala Abda. The hydrographic network is formed by the Oum Er Rbia river and its main tributaries which are the Tassaout, El Abid and Lakhdar rivers on the left bank, which constitutes the foothills of the Atlas Mountains. The right bank, on the other hand, is part of the socalled "phosphate plateau" or Ouardigha Plateau, and is only drained by temporary rivers. The amount of annual rainfall, which varies from 300 mm to 1100 mm and the melting of snow from the Atlas Mountains, mean that the average supply at the mouth of the river is considerable and regular and is estimated at 3680 hm 3 /year, with a minimum of 1300 hm 3 and a maximum of 8300 hm 3 [20][21][22][23].
Fifteen reservoirs are currently built on the river, but the three main dams from headwaters to mouth are Bin El Ouidane, Ahmed El Hansali and Al Massira, which are intended to produce hydropower and supply water for urban use and supply to irrigated agriculture, which is the main consumer of water in the river basin [24].
The study area corresponds to the Oum Er-Rbia river basin (Figure 2), which covers an area of approximately 38,000 km 2 . With such a large surface area, there is great internal variability, as the basin ranges from the Atlantic coast to the High Atlas (4,000 m a.s.l.). The relief are intracontinental folded belts, which form the Atlas Mountains, shortened with areas of plateaus that are poorly deformed tabular domains, generated by the clear dominance of the Alpine orogeny [25]. The climate of the basin is Mediterranean with cold, wet winters and hot, dry summers. Temperature varies between 10 ˚C and 50 ˚C with mean minimum values of 3.5 °C in January and mean maximum values of 38 °C August. Evaporation can reach 1600 to 1800 mm/year. The annual cycle in the entire basin has a clear variation, it is divided into two periods, a wet one (from November to April) and a dry one (from June to September, with July being the driest month), with the months of May and October being transitory. There is a gradient in the basin ranging from arid to semi-arid, with decreasing precipitation from east to west, due to the orographic effect of the Atlas Mountains, with the average precipitation of the basin being 500 mm [20,21,26]. The hydrographic network is formed by the Oum Er Rbia river and its main tributaries which are the Tassaout, El Abid and Lakhdar on the left bank, which constitutes the foothills of the Atlas Mountains. The right bank, on the other hand, is part of the so-called "phosphate plateau" or Ouardigha Plateau, and is only drained by temporary rivers. The amount of annual rainfall, which varies from 300 mm to 1100 mm and the melting of snow from the Atlas Mountains, mean that the average supply at the mouth of the river is considerable and regular and is estimated at 3680 hm 3 /year, with a minimum of 1300 hm 3 and a maximum of 8300 hm 3 [20][21][22][23].
The middle zone of the basin are plains formed by sedimentary rocks under a Quaternary cover. It is an area between the Settat plateau to the south of marl-limestone formation and the Berrechid plain to the north of sandstone-limestone with silt and clay soils. The area is of high Plio-Quaternary tectonic activity due to the existence of major faults [28]. The climate of the region is represented by Beni Mellal, at 500 m a.s.l. at the foot of the mountain. The city has a mild climate with an average temperature of 17.3 ˚C, according to the Köppen classification it is a subtropical Mediterranean climate (Csa) (Figure 3). Annually it precipitates on average 558 millimeters, rainfall occurs throughout the year, except for the summer season in the northern hemisphere. The upper Oum Er-Rbia basin is a mountainous area of sandstone-karstic character, geologically in the domain of the Middle Atlas, with an altitude variability ranging between 1,000 and 4,000 m a.s.l., generated by the convergence of the Eurasian and African plates in the Cenozoic (Alpine Orogeny). It presents a great diversity of relief forms, with large closed depressions, ravines, alluvial terraces and forested mountains which also include the central Hercynian massif in sedimentary Khenifra with intense volcanic activity due to this orogeny that generated granitic intrusions on this, and a high limestone plateau, the Ajdir plateau [22,23,29,30]. The climate of the region can be symbolized by the climate of Azilal at 1,356 m a.s.l., according to Köpen the city has a humid subtropical climate (Cfa) ( Figure 3). The average temperature is 14 ˚C, with cool winters and hot humid summers. Precipitation is spread throughout the year and rains an average of 700 millimeters. In the upper river basin, about 1,100 millimeters are precipitated annually, where it snows on average about 20 days per year [24].
The construction of reservoirs in the river basin begins in the middle of the 20th century with an initial plan for irrigation of one million hectares, which had its origin in the proposals of the French Protectorate in 1936 [31]. The large reservoirs act as regulators and the small reservoirs have the mission of distributing the water through the canals that take it to the irrigable areas. Of the fifteen reservoirs, seven have a capacity greater than 10 hm 3 , whose location is shown in Figure 2 and their data are presented in Table 1. Table 1. Reservoirs with a capacity greater than 10 hm 3 in the Oum Er Rbia river basin, indicating the year of commissioning and the planned extension of their irrigable area. For each one, its abbreviation in three letters is indicated, year of construction and planned irrigable area.

Reservoir
Abrev.  Tassaout  1969  175  Hassan I  HAS  Lahkdar  1986  263 40,000 The largest reservoir is Al Massira and downstream is Imfout which is the closest to the mouth ( Figure 2) and distributes through the canals of the lower coastal area. Near the headwaters is El Hansali and then Ait Messaoud which is the distributor in the middle section in the irrigable area of Beni Amir on the right bank. On the El Abid tributary (on the left) is Bin El Ouidan, the first large reservoir to be commissioned in Morocco, which sends the water through a canal to the irrigable area of Beni Mellal. Of smaller size are the Hassan I reservoir on the Lahkdar river (tributary on the left) and the Moulay Youssef reservoir on the Tassaout river (tributary on the left of the Lahkdar).

Methodology
To make the location map of the region, free access layers of the continents, African countries, cities and provinces of Morocco were searched, the layers were downloaded from Natural Earth (https://www.naturalearthdata.com/). To build the location map (Figure 1), the ArcGIS tool was used. The base layer of the map is the countries of the world, from which the country of Morocco was extracted to highlight the area. This process was carried out from a selection by attributes of the country and exported as another layer, in order to give it a different color from the rest of the countries. After this, the same was done with the layer of provinces of Morocco, the profile of the basic Oum Er Rbia river was placed, and the provinces through which the river passes were selected, after this the provinces were exported and a layer was generated. To mark the most important cities within the provinces, a clip was made, the layer of points of the most important cities of Morocco was crossed with the layer of the provinces that was created previously. From this process, a new layer of dots was generated, which constitutes the most important cities of the provinces through which the Oum Er Rbia river flows. In this map, the river profile corresponds to a discharge of a vector of lines representing the most important rivers in the world. From this layer an extraction was performed to obtain in one layer only the Oum Er Rbia river.
In addition to these, a Digital Elevation Model (DEM) was obtained to make the profile of the basin and the Oum Er Rbia riverbed, with a resolution of 250 meters. To obtain the basin, the first thing that was done was to perform the process called fill, which fills the basin, i.e., the pixels of the DEM are not left without values. After this, a flow direction is made, which estimates how a drop of water would move, in each pixel a direction is drawn following the criterion of greater slope. After this, a point was created at the mouth of the river to obtain the drainage basin. With the point and the flow direction layer, the watershed profile of the Oum Er Rbia river drainage basin was created with the watershed tool. With this process a map of the basin profile with the altitude of the basin was obtained. To create the line representing the river on this map, a photointerpretation was made with the ArcGIS base map.
To conduct the study in the region, once known exact location and the surface of the basin, the Landsat images corresponding to each year were downloaded from the USGS server. Five Landsat scenes are needed to complete the Oum Er Rbia basin area. The study starts in the year 1985 with Landsat 4 images and ends in the year 2020 with Landsat 8 images. It was decided to be done in five-year intervals, so a total of 40 images from Landsat 4, Landsat 5 and Landsat 8 satellites were needed. The downloaded images with dates can be seen in Table 2. All images were taken in the same season, between April and September.
When all the images were downloaded, they were processed with the SNAP tool (Brockmann Consult Gmbh, Germany), with the mosaicking tool the five images that make up the study area for each year were joined. Thus, there is one product for each year, i.e. eight products, made up of five satellite scenes for each product.
After having all the products for each year, we proceeded to calculate the vegetation indices. In this case, the three most common indices, known by their abbreviations NDVI, TNDVI and RVI, were tested. NDVI (Normalized Difference Vegetation Index) is a classic indicator of the greenness of biomes, and is useful for understanding vegetation density and assessing changes in plant health. It is the most widely used index of vegetation globally, as it compensates for changes in illumination, land surface slope among other things. NDVI is calculated as a ratio between red (R) and near infrared (NIR) values [14]. Its formula is as follows: This determines values between -1 and 1, with values closer to 1 being those with a high density of vegetation and those closer to -1 being water, snow, and around 0 being bare soil or rocks.
TNDVI (Transformed Normalized Difference Vegetation Index) is the transformation of NDVI so as not to work with negative values. TNDVI indicates a slightly better correlation than the original one between the amount of green biomass found in a considered pixel [32]. Its formula is: Landsat 8: As a result, the lowest values are found in built-up areas and water bodies [32]. The RVI (Ratio Vegetation Index) is an index that highlights the contrast between vegetation and soil, because it is sensitive to the optical properties of the soil, but not to lighting conditions. It is mostly used to estimate leaf area biomass [33]. Its formula is as follows: When the values are high they reflect the presence of vegetation, and when they are low they represent soil, water, etc. [33].
With respect to the Google Earth Engine (GEE) tool, it is accessed through a web service (https://explorer.earthengine.google.com/#workspace) in which after registering one can access various types of products already calculated from Landsat 5 images in the period 1985-2010 and Landsat 8 in the period 2013-2020. One of the products is the vegetation index EVI, calculated according to the following formulas [19]: Landsat 8: After performing the calculation of each index for each annual product, 32 products remain in total, i.e., four products from each year. Then the same intervals were adjusted for each index, i.e. the products of each year of the same index have the same intervals to be able to compare with each other.
To calculate how much vegetation is in the watershed we need to insert the vector layer that was created that delimits the watershed. Once inserted in each product, we proceeded to create the histogram of the image using the vector layer, we calculate the number of pixels of each type that are in the watershed area. When the histogram was calculated, the data were extracted from which a table was created showing the number of pixels that made up each interval. Thus it was possible to calculate the area in square kilometers of each class, since the resolution of the image is known, which is 30 meters, so each pixel corresponds to 900 square meters.
After having the results, the next step was to present the figures created from the indexes. This presentation was made from ArcGIS, the SNAP figures were exported to GeoTIFF format, so that they were georeferenced. They were opened in the ArcGIS tool and the extract by mask tool was used, which extracts the data from the raster image that we had exported through a vector layer, in this case the watershed profile.
For GEE products, the area of interest is exported in GeoTIFF format directly in the web application interface, selecting the region of interest (ROI) in order to reduce the file size and obtain the definition as a 30 x 30 m pixel Landsat. The file is then treated as those obtained with SNAP.

Results
From the location of the study site, the Landsat satellite images needed to construct the watershed were obtained. In total there are five for each year, corresponding to the 201-37, 201-38, 202-37, 202-38 and 203-37 strikes and passes. Table 2 shows the dates on which the images were taken for each of these five in each year. The images were downloaded between 1985 and 2010 from the Landsat 4-5 TM satellite and those from 2015 and 2020 from Landsat 8 OLI.  Figure 4 shows the result of joining the downloaded images after the mosaicking process. The result shows the union and overlapping in some areas, as well as the irregularities due to the composition of several images from different days. Figure 5 is the cutout of Figure 3 with the basin profile created in false true color, i.e. the Landsat 8 red, green and blue bands corresponding to band 4, band 3 and band 2, respectively, have been applied. The colors of vegetation, water bodies, cities and the phosphate mining area can be seen. Figure 6 is the same as Figure 5 in which the bands have been combined to create a false color. In this case the near infrared (NIR), red and green bands have been combined, which in the case of the Landsat 8 satellite corresponds to bands 5, 4 and 3 respectively. In this combination, the bands respond sensitively to green vegetation due to the infrared reflectivity, which is very high for these bands, as opposed to the visible reflectivity, which is low. The red colors in Figure 6 represent healthy vegetation, in which the vegetation of the mountain area stands out in a darker red color, while the lighter red color highlights the crop areas that can be observed in the areas near Beni Mellal and Beni Amir. The pinker areas correspond to areas of less developed vegetation and the brown or golden colors are areas of scrub or grassland. A white color can also be seen, indicating areas with little or no vegetation, and, finally, a very dark blue or black color indicating the presence of lakes, as can be observed in reservoirs [34].   7 also corresponds to a false color, in this case it is a combination of the shortwave infrared, near infrared and blue bands, which in the case of Landsat 8 are bands 6, 5 and 2 respectively. The purpose of this figure is to highlight even more the cultivated areas, which as they are the same as those discussed above, and which in this one stand out with a bright green. While the mountain vegetation is dark green, and again stand out in golden, bluish and pink colors the areas with little or no vegetation [34]. The agricultural area of Beni Mellal (left bank of the Oum Er Rbia) and Beni Amir accounts for about 3,000 km 2 approximately, estimated visually according to the Google Earth measuring tool.

Results of the vegetation indices
Once all the images of each year joined together, we proceeded to calculate the four indices for each year. From the total of the eight dates calculated, the first and last dates (1985 and 2020) and one in between (2000) were selected. Therefore, the images of each index are presented in order to appreciate the changes, all on the same scale of representation. Figure 8 shows the results obtained. The year 2000 is the year with the smallest vegetated area, while 1985 is in an intermediate coverage and 2020 presents the maximum coverage of the dates studied. The visual coincidence of the results of the four indexes is highlighted. However, the TNDVI and EVI figures appear as those that best represent the situation.
The indices clearly show the three vegetated zones. In the center of the image, the irrigable zone of Beni Mellal, and to the south of it the forest zone of the Atlas Mountains, whose summit areas appear without vegetation above 2,500 m. To the west is the basin near El Jadida, close to the Mediterranean Sea.
In 1985, areas of sparse or almost no vegetation dominate, but on the contrary, there are a large number of green regions representing vegetation. In the TNDVI index, the color of the reservoir is not so clear and is blurred with the color representing no vegetation. But the cultivated areas stand out very well, because it represents them with intensity. The cultivated area near the mouth at El Jadida is well marked. It should be noted that in the year 2000, there is an increase in the areas of no vegetation. There is a clear decrease in the region representing vegetation in general, but especially the cultivated area disappears from the area south and east of El Jadida. There remain areas of vegetation in the province of Kenifra and Beni Mellal, among which are included the cultivation area and the Atlas Mountains forest area.
Finally, in 2020, the area of no vegetation continues to increase, especially in the western part of the basin. But there is also a decrease in the intermediate zone and a clear increase in the green surface and therefore in the vegetated area. The irrigable zone of the middle section of the river increases in vegetation intensity and also in the forest zone.
In order to know the exact area of vegetation in each year, a histogram has been calculated for each index of each year, using the histogram analysis function of SNAP. Thus, the histograms of each year have been calculated to obtain the numerical information and the pixels have been grouped in class. Figure 8 shows the histograms for the period 1985 -2020 of the EVI index. Figure 9. Histograms of the EVI index on the dates studied.
From the data of the histograms of the four indices, it has been possible to calculate the areas of each class in each index, and therefore to know the vegetated area in each year (tables 3 to 6).    Observing the four indices, despite their obvious differences, they are correlated with each other and the representation in classes is similar, as in the case of 2020 for example (Figure 10a). The most correlated indices are TNDVI and EVI (Figure 10b) whose relationship between them presents a coefficient of determination of 0.9706.
Regarding the periodic series between 1985 and 2020 in the river basin, in the intermediate period the vegetation area decreased and the null vegetation zone began to increase. At present, this null zone continues to increase, but so does the vegetation area. Figure 9 shows the distribution of EVI values, where values lower than 0.07 indicate arid zones, while values higher than 0.14 indicate vegetated zones. In drought years, distributions are observed more pointed in the low values, while in favorable years the distribution loses pointing and widens the tails, indicating that although there are areas of aridity, there are also areas with abundant vegetation.

Expansion of agriculture
The expansion of agriculture from satellite images in the middle basin of the river, corresponding to the province of Tadla Azilal, is clustered around the towns of Beni Mellal on the left bank and Beni Amir on the right bank, currently comprising an agricultural area of about 4000 km 2 . The increase from 1985 to the present has not been what was expected in the initial irrigation plans (about 10,000 km 2 ). In Figure 11a and 11b it can be visually observed how the irrigated area has increased due to the increase in vegetated area with high EVI values. The scrub and Mediterranean forest area appear in the mountainous part of the Atlas and has also increased its vegetation density.  While the aridity zone has been maintained in these years, the first vegetation class has decreased dramatically from 1985 to 2020, and this area has contributed to the increase of the other classes of the index already related to the presence of vegetation.
In the Beni Amir environment, landscape change has been significant in this study period. The ranges of values of the EVI index make it possible to distinguish the cultivated areas according to the existing vegetation. Figure 11a shows a regular distribution of fields around the town of Beni-Amir, where regular mosaic farming was carried out in 1948, interspersing plots of arable crops with woodland (mainly citrus) and fallow land [2]. At that time, irrigation was about 134 km 2 and was limited to water pumped from the river in the absence of the regulation of the El Hansali reservoir. When it was transformed into irrigated land with full availability of water with the construction of the reservoir and the extension of the irrigation canals, the regular mosaic was lost as shown in Figure 11b; the fallow land disappeared and became an area of varied herbaceous crops. In addition, new cultivated areas of herbaceous crops are created with pivot irrigation to the west, destined to forage plants for local use, whose circles are perfectly visible, some of which are larger than a kilometer in diameter; in the northern zone other large extensions of crops have also been created.

Discussion
Satellite images provide a temporal view of the planet's surface that is changeable over the years, and remain as a record of a past time, allowing to elaborate studies in long time series of changes in uses or natural changes [35]. The existence of the Landsat series is an unquestionable utility in current studies, as well as the possibilities of processing in large territories, being used by numerous published studies, such as Zadbagher et al. [36] who studied the change of land use and land cover by Landsat images in Turkey; the studies of forest changes by Banksota et al. [37]; the studies of temporal variations in wetlands by Kayastha et al. [38]. Specifically on agrarian issues, work has been developed all over the planet, such as studies of abandonment of cultivated areas in the Caucasus [39], studies on crop cycles in the Amazon in Brazil [40], or the spatio-temporal reconstruction of the development of irrigation in Uzbekistan [41].
Vegetation indices have been used in the study of agricultural areas since they first appeared, due to their economic importance. The first works with the first satellites that came into service, such as studies of wheat growth using the LAI leaf area index [42], water stress in wheat cultivation [43] and the effect of pest attacks on forests [44]. The TNDVI index has shown its usefulness in studies in agricultural ecosystems and among recent works, it has been used to study the sustainability of agriculture in Iraq using remote sensing and geomatics techniques [45], monitoring carbon sequestration in orchards [46], or spatiotemporal studies on aridity in the Punjab of Pakistan [47]. The use of EVI has limitations due to the influence of relief and topography [48], where the presence of shaded areas may give lower values of the index than would be appropriate. In our case, since the agricultural area is in a plain, this problem does not exist. Many studies with the EVI index have been performed with MODIS images [49] where index values were found for an area of Matto Grosso in Brazil varying in the range of 0.20 to 0.65 for agricultural areas, 0.25 to 0.50 for grasslands and scrublands and 0. 50 to 0.60 for forests; these values are coincident with those obtained in our study, although in the agricultural zone of Tadla Azilal the agricultural zones start at values of 0.15, at times when certain crops are in their initial stages. These lower values are in agreement with those obtained by Chen et al. [50] in cornfields in Mexico, where the minimum values were as low as 0.14 and the maximum values exceeded 0.60. Actually, the minimum values are debatable without field work in which vegetation measurements are made and the EVI index is estimated to observe the values in the absence of vegetation; in that sense, the work of Tran et al. [51] points out the value of 0.20 as the lower limit of vegetation when using Landsat-8 images. But the study by Gosh et al. [52](2021) who made comparisons of EVI index between Sentinel-2 and Landsat-8 imagery in agricultural areas of India, indicates as low values as 0.10 for fallow land, with normal values between 0.15 and 0.25, while values in cultivated land are between 0.15 and 0.35 mostly. Therefore, we consider our class separation values to be adequate for our study.
The changes that occurred in the agricultural extension in Tadla Azilal during the 20th century changed a territory dedicated to pastoralism and subsistence agriculture into agricultural areas dedicated to extensive cultivation first from 1940 [2] as a result of the French colonization plans. With independence of the country in 1956, the Moroccan new government continued the policy of agricultural expansion in view of the agreements with the European Union, continuing with the construction of reservoirs and new irrigable areas. At present, of the slightly more than 10,000 km 2 planned for irrigable area in the middle section of the Oum Er Rbia valley, about 4,000 km 2 have been deployed, and in the total basin about 25,000 km 2 are vegetated, either by crops or areas of scrubland and forests in the Atlas area and the cultivated areas near the coast in El Jadida and Azemmour.
The changes in land use in Morocco are comparable to those in other Mediterranean areas such as the Campo de Cartagena (Murcia, Spain), where a process of transformation to irrigation has taken place with surface water inputs (although of smaller dimensions) affecting about 200 km 2 . However, other agricultural territories suffered this transformation more than 1000 years ago in the Iberian Peninsula. The Huerta de Valencia was transformed into irrigated land in the 10th century by the channeling of surface water from the River Turia, as were the irrigated lands of the Segura in Murcia and Alicante (Spain). The banks of the Júcar, however, are later and began in the 13th century and ended in the 18th century as consolidated irrigation [53]. However, there is an important difference between these areas, namely the great extent of agricultural transformations in Morocco. While all the combined irrigated areas mentioned in Valencia and Murcia account for about 1,300 km 2 , in Tadla Azilal county there are already about 4,000 km 2 . In addition, the extent of irrigation in Morocco, considering the other river basins, could consolidate up to 100,000 km 2 if it reaches the maximum expected from the hydraulic constructions carried out in Morocco since the mid-twentieth century [4].
Regarding the usefulness of the vegetation indices in this case study, we consider the TNDVI as the best index to evaluate the vegetation, mainly because it allows you not to work with negative numbers, and because it highlights the vegetation better as seen in the figures of the indices, since it is the only one in which the cultivated area near El Jadida does not disappear. However, when the processing of the images with SNAP is long and costly and requires sufficiently powerful computer equipment, its use would not be advisable. In that sense, the use of the Google Earth Engine platform has the advantage of having the data already calculated without having to locally process the information, but it has the disadvantage that it only offers the EVI index for the moment (summer 2021) and besides you cannot choose the images specifically. This index is also suitable especially in flat terrain, where there is no influence of topography. Contrast with field data is necessary in future work to verify the minimum values of vegetated terrain.

Conclusions
The increase of the agricultural area in Morocco and specially in the Oum Er Rbia river valley has been due to the availability of regulated water and the construction of hydraulic works for its conduction to areas to expand the irrigation fields and increase agricultural production. We evaluated by means of several vegetation indices applied to satellite images, how the agricultural area has been changing from the 1980s to the present, comparing for the first time the use of the Google Earth Engine with local processing by SNAP software. We found that there are not differences in terms of the results obtained, but the local downloading and processing of Landsat images is time-consuming and the advantages of the cloud download and processing by Google is better if the need is covered. The increase of agricultural area is important is the study site and has changed from a production related to national demand to an intensive production for export to European Union.

Conflicts of Interest:
The authors declare no conflict of interest.

Appendix A
The appendix is an optional section that can contain details and data supplemental to the main text-for example, explanations of experimental details that would disrupt the flow of the main text but nonetheless remain crucial to understanding and reproducing the research shown; figures of replicates for experiments of which representative data is shown in the main text can be added here if brief.