Modeling of Urban Growth Using Cellular Automata and GIS Case of Benslimane in Morocco

In this study, our goal was to research land-use change by combining spatio–temporal land use/land cover monitoring (LULC (1989–2019) and urban growth modeling (1999–2039) in Benslimane, Morocco, to determine the effect of urban growth on different groups based on cellular automata (CA) and geospatial methods. A further goal was to test the reliability of the AC algorithm for urban expansion modeling. To do this, four years of satellite data were used at the same time as population density, downtown distance, slope, and ground road distance. The LULC satellite reported a rise of 3.8 km2 (318% variation) during 1989–2019. Spatial transformation analysis reveals a good classification similarity ranging from 89% to 91% with the main component analysis (PCA) technique. The statistical accuracy between the satellite scale and the replicated built region of 2019 gave 97.23 %t of the confusion matrix overall accuracy, and the region under the receiver operational characteristics (ROC) curve to 0.94, suggesting the model's high accuracy. Although the constructed area remains low relative to the total area of the municipality's territory, the LULC project shows that the urban area will extend to 5,044 km2 in 2019, principally in the western and southwestern sections. In 2019–2039, urban development is expected to lead to a transformation of the other class (loss of 1,364 km2), followed by vegetation cover (loss of 0.345 km2). In spatial modeling and statistical calculations, the GDAL and NumPy Python 3.8 libraries were successful.


Introduction
Urbanization is a significant symbol of the growth of science, technology, and the greater capacity of humanity to reform the natural world. It is also one of the required measures to modernize a country [1]. In the local natural environment, the rapid development of urban areas has an impact on the complex socio-economic and natural systems, such as deforestation, air pollution, and agricultural land reduction [2,3]. To ensure sustainable urban growth [4] and to explicitly understand the spatial distribution of urban areas and spatio-temporal patterns, it is important to find an effective method to simulate and model urban expansion [5][6][7]. A variety of techniques and models have been developed over the last three decades to explain the dynamics of urban growth processes [8,9]. The introduction of the geographic information system (GIS) into the modeling of urban expansion was important at the end of the 20th century [10,11], to understand the predictor variables of spatio-temporal changes. Regression models, including [12][13][14][15][16], artificial neural networks (RNA) [17][18][19][20][21][22] , agent-based models (ABM) [23,24], are used for their efficient spatial computational capabilities. Cellular automata (CA) models are widely considered to be the most functional instruments [25][26][27][28][29][30],Tobler first suggested the urban AC model in 1970, to simulate Detroit's urban growth. In San Francisco, Clarke et al. [31] simulated shifts in land use using the SLUTH CA model. Batty and Xie [32], using the CA model, researched urban sprawl in Baffalo. To research the urban growth of Guangzhou, Wu [33] developed a logistic regression CA model. In northern China, He et al. [34] simulated land use shifts by integrating the model of system dynamics and the model of CA. The urban growth of Dongguan was simulated by Li and Gar-On Yeh [35], using data mining and AC. Li et al. [36] applied the GPU technique to the CA model to simulate improvements in urban land use in the province of Guangdong. Mustafa et al. [37] merged cellular automatons and agent-based approaches in Wallonia, Belgium, to model urban growth from 1990 to 2000 [38]. The urban planning of the cities of the Pearl River Delta was analyzed by Liu et al. [39] by combining the Landscape Expansion Index (LEI) and the CA model. In the last three decades, research using CA urban growth models have proliferated and provided in-depth insight into urban growth processes, which contributes to the development of strategic strategies for urban development to ensure a sustainable urban future [29]. The urban modeling and forecasting capabilities of CA models have been demonstrated and has been used successfully over the past fifteen years worldwide [40]. In this study, we aimed to examine the spatio-temporal dynamics of land use/land cover (LULC) over the past decades in order to model urban development for the municipality of Benslimane, Morocco using the CA algorithm. In this analysis, the LULC of four separate years: 1989, 1999, 2009, and 2019 was evaluated in order to achieve this objective. We also considered other criteria and growth parameters, such as topography, road proximity, proximity to the city center, demographic statistics, as well as areas with growth constraints. We estimated of the building scope for the future. The prediction of the future was useful in extrapolating the environmental effects of the urbanization of the town of Benslimane, which could help improve sustainability. In our analysis, we also report on the CA algorithm's accuracy and reliability in modeling the urban development of the city of Benslimane in Morocco. Benslimane is medium in size and in recent decades has undergone urban development and essential population influx. This development has contributed to a range of urban, environmental, and socio-economic issues, including the depletion of air quality, groundwater quality, temporary ponds, and green spaces. In recent decades, the urban area of Benslimane has expanded, so it seems appropriate to examine the urban growth trend and deduce its potential impact on other land use characteristics.

Study Area
The analysis was carried out on the territory of the town of Benslimane (Figure1), the capital of the agricultural province and hinterland, bordered to the west and north by the municipality of Ain Tizgha, to the south by the municipality of Ziaida, and to the east by the municipality of Oulad Yahya Louta. The high altitudes in the north-west range from 175 to around 390 m in the south-east. The area of the town of Benslimane is 72 km 2 , the polar coordinates are 33°36N latitude and 7°06W longitude, the average Lambert coordinates are x = 342,000, y = 336,000. The population of the town of Benslimane is approximately 58,000 Hah (RGPH 2014).

Climate
Benslimane has the following characteristics: a temperate and humid climate along the coastal regions. A dry climate with cold and heat alternating, with a continental climate when entering the mountainous regions and in the south of the province [41].

Temperature
Benslimane's annual average temperatures are 23.7 °C for the maximum and 10.3 °C for the minimum. As one travels away from the Atlantic coast, there are considerable variations. The average annual temperature of the coastal zone is 17.5 °C, with a mean temperature not exceeding 32 °C. The low interior plateaus show very high thermal amplitudes, but without sudden variations. The average annual temperature on these plateaus is 18.5 °C, with a high of 40 °C. [41].

Rainfall
The average annual rainfall recorded in the province is 350 mm. The average annual precipitation measured for the period 1977-2008 for the Feddan Taba and Malleh Dam pluviometry stations is approximately 310 mm. The average monthly distribution of rainfall suggests the presence of two distinct seasons: the wet season from October to May, during which nearly all rainy episodes occurred (86 to 92% of the annual rainfall); the dry season from May to September, with just 8 to 14% of the annual rainfall [41].

Methods and Data Used
This research involved numerous thematic layers prepared using satellite images, Benslimane census data, and topographic sheets. Land use/land cover maps (LULC) were prepared to consider the dynamics of urban development from 1989 to 2019 using satellite datasets. Different factors favoring urban growth with LULC were prepared as raster layers for modeling purposes, namely population density, slope, proximity to road networks, and proximity to the city center (Figure2). A model was carefully tuned to keep the simulated LULC both statistically and spatially similar to the real LULC to achieve the best possible threshold values. For a development forecast for the year 2029, the pattern of the threshold values obtained was evaluated. In ArcGIS, the LULC, roads; other thematic layers were processed in Python 3.8.

Thematic layers
LANDSAT multi-temporal satellite images for the years 1989, 1999, 2009, and 2019 were downloaded from the USGS website (https:/earthexplorer.usgs.gov) and were used to prepare LULC maps, using the classification to obtain various groups-constructed, vegetation, water area, and others. Tables 1 and 2 present a summary of the satellite data used and of the major LULC groups. The goal of this research was to design the built feature class, which includes constructed areas such as residential areas, companies, institutions, and informal settlements. Land was classified as vegetation primarily if it was under any form of green cover. Streams, etc., were included in the water zones, and the remainder was listed as other, mainly agricultural land. Satellite maps classified by the LULC were validated using field checks and Google Earth photos. For each of the observation years (1989,1999,2009, and 2019), the average accuracy was measured and ranged from 80 percent to 89 percent, while the kappa coefficient ranged from 0.73 to 0.85. To examine LULC transition and urban sprawl over 1989-2019, and then to calibrate the model, the finalized LUC maps were adopted. To calculate the population density for 2024 and 2030 using the unit process, the thematic layers of population density by region were prepared. Demographic data from the Moroccan census (1994,2004, and 2014) were taken into consideration for the simulations. During the model calibration method, population density maps for 1994, 2004, and 2014 were used and, for future forecast purposes, the population density map for 2024-2030 was used. To fit the geometry of the other participating layers, the density was determined in the vector layer attribute of the different areas of Benslimane, which was later converted to a raster layer of cell size 30 x 30 m. In all three groups (<3, 3-5, >5), the digital elevation model (DEM) with a spatial resolution of 30 m was used to prepare slope maps.  The main road network was obtained from the satellite picture showing proximity to the main roads, grading from 250 m to 3 km intervals. Similarly, in 8 classrooms, multiple buffer circles across the city center location were sampled at 1 km intervals. The vector layers were rasterized to a cell size of 30 m and cut to match the geometry of the old layers with a constant magnitude. A layer for restriction zones was considered to enhance the model's ability to replicate the current growth system. The model did not produce pixels accumulated in restricted areas representing the forest and habous lands by this binary layer.

Model Algorithm for CA
To take into account all the variables that lead to urban growth in Benslimane, the script for the CA model was formulated. A kernel of 3*3 size kept the test pixel in its middle, as shown in Equation The model depends primarily on the current state of the test pixel, the current state of the immediately adjacent pixels, and the set of rules for transformation [42]. Matching layer geometry is important to ensure that in all raster layers, every random pixel (ai, j) represents the same piece of ground. According to the following equation [42], the dependence of the future state (t + 1) of a pixel on the passage rules (Ø) and the normal state of a pixel was examined: The passing rules (φ) were a function of a set of threshold conditional statements represented by: Ø=f (T,B) (iii) Equation (iii) suggests that the transition rules are a function of T, a set of threshold values for all the affected parameters, and B, a set of kernel count values associated with each set of T.
where the threshold values for road proximity, city center distance, population density, and slope value are TR, TC, TP, and TS, respectively; BR, BC, BP, and BS are the corresponding number of pixels incorporated in the test kernel for each element belonging to T. The laws of the adopted CA model comply with the existing general conditions in the real world: (a) Built-up land and water area classes were exempt from removal; (b) vegetation land or other classes can, depending on the threshold values (T) and the number of neighboring built-up pixels (B), move to built-up land if the test pixel falls into the unregulated area layer category [43].

Calibrating the Model
Using the LULC satellite data for 1989 to simulate the year 1999, the model calibration was adopted, then the LULC satellite data for 1999 was used to simulate the year 2009, and finally, the LULC satellite data for 2009 was used to simulate the year 2019. By simulating the LULC of time t2 using the LULC of time t1 and the driving parameters, the model deduced the best-defined threshold values, which represent the closest result to the real world [42], both statistically and spatially. The threshold values of the four factors (TR, TC, TP, and TS) and their corresponding integrated pixel count values (BR, BC, BP, and BS) were obtained through trial and error. Moreover, the script also tracked the contribution of each factor influencing the new integrated pixel generation. Threshold values were plotted after calibration, and trend lines were used to project thresholds for the future stage to estimate the extent of future accumulation [43].
Principal component analysis (PCA) was used to determine accuracy at each point of the simulation [43], where the difference between two images was determined by spatially subtracting the cumulative time pixels t2 from the corresponding time pixels t1. Here, from satellite and simulated LULC images from the same time, dichotomous made up layers were extracted. In the event of a change, the built-up area common to the two images was eliminated by subtraction and generated non-zero values (1 or −1). For a precise percentage measure, the proportion of these non-zero values relative to the total number of pixels accumulated from satellite images was considered.

Results and Discussion
To deduce land use/land cover shift and urban sprawl patterns in Benslimane for 1989-2019, LULC satellite maps were prepared and analyzed. Later, to calibrate the model and predict the potential scale of the agglomeration, the key factors leading to urban growth were analyzed. To provide the best possible outcomes, care was taken to calibrate the model.

Mapping and Urban Development of Land Use/Land Cover
For 1989-2019, multi-temporal satellite data were used to chart land use/land cover at 10 year intervals. The real (satellite) LULC was named after this. The data was then used for 1999-2039 to simulate the LULC-the simulated LULC. During the model calibration period (for 1999-2019), these simulations were correlated with the satellite LULC.

Satellite-Based Land Use/Land Cover Mapping and Urban Growth
We found that the real built-up land (by satellite) increased from 1,206 to 5,044 km 2 in 1989-2019 with a 5,305% variation. Periodic observations show that the built-up region in 1989 was 1,206 km 2 (1,667% of the total area), which rose to 2,147 km 2 (2,969%) in 1999 with a rise of 1,301% (Table 3). It increased to 2,999 km 2 (4,146%) in 2009 with an 1,177% rise. In 2019, the built-up area increased to 5,044 km 2 (6,972%) with a 2,827% growth. The overall variation in the vegetation region was observed at 6,464% between 1989 and 2019. In 1989, the vegetation cover covered 21,392 km 2 (29.57%) with a variation of 9,601%, which increased to 28,337 km 2 (39.17%) (1999). It then increased in 2009 to 31,223 km 2 (43.16%) with a 3.989% variance. In 2019, with a variation of −7.127%, the vegetation cover decreased to 26,068 km 2 (36.04%). In the class of water characteristics from 1989 to 2019, the cumulative variation observed was 0.173%. In 1989, the area protected by the water zones was 0.077 km 2 (0.106%), which rose to 0.266 km 2 (0.367%) with a 0.261% variation in 1999. This area increased to 0.313 km 2 (0.433%) in 2009 with a 0.066% variation. In 2019, with a variation of −0.154%, the total area under the water zones rose to 0.202 km 2 (0.279%). In the other function grades, the average decrease was −11,942%. In 1989, all the remaining objects classified as others occupied a total area of 49.66 km 2 (68.65%), which dropped to 41.59 km 2 (57.49%) in 1999 with a −11.164% improvement. In 2009, with a difference of −5,232%, the area of the other class decreased to 37.80 km 2 (52.26%), then to 41.02 km 2 (56.71%) in 2019 with a difference of 4.454%.
Our research illustrates the impact of urban sprawl from 1989 to 2019 on other LULC groups. Spatio-temporal mapping of LULC demonstrates that built-up land was largely concentrated in the central Benslimane area. Later, because of the availability of suitable sites and the advancement of the new development plan, the building was expanded to the south and southwest, resulting in a population increase. Since built-up areas determine the basic human-environment dynamics in urban environments, unbuilt elements (i.e., vegetation, water areas) are intrinsic components of built-up areas.

Urban Expansion Modelling
Urban growth modeling was carried out using an earlier LULC map and various contributing factors, close to major roads, city center, population density, and slope data, to forecast the future urban growth model for the year 2039.

Contributing Factor in Urban Development
The study indicates that the densification of the agglomeration and the formation of subdivisions in remote areas resulted from an irregular mixture of the impact of the road, the town center, and geomorphic parameters. The following sections discuss the different contributing parameters Land Use/Land Cover The probability of transformation of the various groups of LULC into built-up land is different. To this end, the group of others was given a higher rank, followed by vegetation. In modeling, constructed-up land was considered to be a non-transformable class.

Proximity to Major Roads
The road map shows the density of roads in the middle and south of the city, and the density in the northwest and southwest is lower. The density of roads in most suburban areas is poor ( Figure  4a). Given the role of main roads in the urbanization process, the high accessibility of areas closer to roads has given growth priority and vice versa.
Proximity to the Town Center As a major multi-urban functional area, the city center functions. To track non-uniform development in all directions, proximity to roads and the center was taken into account. Distance from the center has primarily affected urban growth in the central region, and as distance rises, its impact has diminished.  1999 (a, d, g), 2009 (b, e, h), and 2019 (c, f, i).

Slope
The analysis indicates that the relief ranges between 177 and 339 m relative to the Benslimane mean sea level (NGM). The Benslimane slope map shows that most of the region has a slope of less than 5°, but small parts have slopes greater than 5° to the north, middle, southeast, and southwest (Figure 4 (c)). The slope of the land is a factor that has considerable influence over the growth of buildings. High landforms and steep slopes cause greater runoff and make the land less suitable for construction, whereas low landforms and gentle slopes provide buildings with a comparatively more stable foundation. This could explain why the growth of low slope areas was favored.

Validation of the LULC simulated
An accuracy assessment revealed a very high statistical and spatial similarity (> 89%) between the simulated LULC and the satellite LULC for the known years (1999-2019), as shown in Table 4. The percentage of spatial accuracy and statistic was calculated with respect to the spatial location and built extent of built-up land in satellite LULCs. The calculated confusion matrix has an overall accuracy of 97,23% (Table 6), and the area under the curve (AUC) of the receiver operating characteristics (ROC) graph ( Figure 8) has been shown to be as high 0,94. The high precision of the confusion matrix and the AUC indicates high precision of the modeling. This is probably due to the use of an average prediction period (10 years) at regular intervals during the calibration as well as the proportion of the frame which represents only 7% of the overall area. The transition probability matrix was calculated using temporal land use / land cover maps (Tables 7 and 8

Influence of Potential Urban Growth on Various LULC
A net rise of 1708 km 2 in urban areas and an urban growth rate of 33.88 percent between 2019 and 2039 was demonstrated by the built-up expansion shown by the CA-based urban growth model (Table 5 and Figure 9). This will affect the other categories of LULC (net loss: 3.32 km 2 ) and the coverage of the plant (loss: 1.3 km 2 ). The research shows that, during the period 2019-2039, the urban growth rate in Benslimane will not be high. Table 4. The statistical variability as per satellite-based and simulated photos of built-up land.

Conclusions
The research covers the spatio-temporal monitoring of Benslimane's LULC (1989-2019) and urban growth modeling (1999-2039) to glean the past and future urban growth paradigm and its effect on the land use/coverage variety using the CA model. During the period 1989-2019, the simulated and satellite-based LULC satellite experienced rapid urban growth (net increase of 3,838 km 2 and 4,327 km 2 ), resulting in a land-use transition in Benslimane. The other groups were modified by population development-17.4%, followed by vegetation cover (20.13%) and water level (1,62% variation). The simulated LULC based on AC shows that the urban area will expand to 5,044 and 6,752 km 2 in 2029 and 2039, respectively, mainly in the western and southwestern sections. In the period 2019-2039, urban development will replace and transform the "other" category of LULC (loss of 1.36 km 2 ) and the category of vegetation cover (loss of 0.34 km 2 ).
In the simulated and satellite LULC over the described years (1989-2019), the dynamics of transition were parallel. For the years 1999 and 2019, the study of the spatial variance performed using the PCA technique showed high construction precision. The overall accuracy of the confusion matrix (97.23%) and of the area under the receiver operational characteristics (ROC) curve (0.94), measured using the satellite-built patch and simulated for 2019, underlines the high modeling accuracy.   An attempt was made to demonstrate the relationship between variables related to the phase of spreading over time , which highlighted the strong impact over time of proximity to built-up areas and the city center during the initial observation period, although the population density contribution increased later. Table 6. Confusion matrix for actual and simulated 2019 built-up pixels. Table 7. Transition matrix (satellite-based 1989 to 2019). The calibrated CA model replicated the process of sprawl in the real world statistically, while the spatial analogy was mild. In the case of the town of Benslimane, where urban sprawl is often spontaneous and unforeseen, this high spatial accuracy shows the limits. To calibrate the model to forecast urban sprawl for years to come, this work centered on increasing the precision of the CA model and used data from four uniformly spaced expansion times.