The alteration of LULC has a significant influence on ecological components (Bhungeni et al. 2024; Z. Zhao et al. 2024). This provides essential ecosystem services that are advantageous to both society and biodiversity. These services include the supply of drinkable water, support for farming activities, recreational opportunities, and the conservation of natural habitats. However, human-induced changes in LULC pose a substantial risk to the preservation of natural land cover and freshwater resources in various countries, particularly in watershed areas. Examining alterations in the Earth’s surface is crucial for understanding ecological and societal transformations (Peterson and Levinson 2020). Traditional approaches and aerial imagery for data collection are time-intensive and often lack sufficient accuracy. The latest analysis techniques for satellite images now allow for more efficient detection of changes in LULC (Li et al. 2019). In recent years, remote sensing imaging has been increasingly utilised to identify changes in LULC and vegetation. The abundance of historical data and remote sensing imagery has made it convenient to study the impact of human activities on LULC (Ali and Johnson 2022). Classifying LULC is essential for detecting changes, and the approaches and strategies used for classifying LULC and extracting precise data from remote sensing images are highly adaptable. Mapping and monitoring LULC changes are essential for gaining insight into the spatial distribution of human activities on land and their impact on the natural environment (Dong et al. 2019). This knowledge is crucial for water resource managers and environmental health practitioners to formulate efficient environmental strategies.

Management plans and regulations. Remote sensing techniques have largely replaced conventional field observation methods for LULC mapping. They offer a cost-efficient and effective way to gather spatiotemporal data across large areas. Remote sensing data can be systematically recorded, stored, and shared, making it a valuable resource for studying different landscapes. Many studies have shown that machine learning methods are very effective at classifying LULC using remote sensing data (Andrew et al. 2023).

The extensive body of research on LULC dynamics in Morocco is primarily limited by the lack of a uniform spatiotemporal scale and a generalised methodology. Most of the existing literature on the topic is either localised or methodologically constrained or temporally disjointed. This is often combined with the use of conventional classification methods such as maximum likelihood or minimum distance. These limitations make it difficult to ensure reliability due to the possibility of misclassification and low transferability across sensors. Most studies did not combine multi-temporal Landsat datasets from both Landsat 7 and Landsat 8 with uniform preprocessing, spectral harmonisation, and atmospheric correction, especially in the Berrechid–Settat region. This region is characterised by high rates of urban expansion and agricultural intensification but is exposed to increased pressure on water and land resources. Most of these studies did not include a comparative assessment of machine-learning methods using the same training and validation datasets. They have almost focused on assessing a single model or did not report significant statistical differences in accuracy metrics, which include Overall Accuracy, Producer’s Accuracy, User’s Accuracy, and Kappa. Finally, all regional studies demonstrated a lack of data due to the absence of efficient GEE-based computing for extensive, cloud-based analysis and standard temporal mosaicking. Consequently, this study tried to respond to this lack by (i) combining imageries from Landsat 7 ETM+ and Landsat 8 OLI within the same GEE workflow; (ii) rigorously comparing three robust classifiers, SVM, RF and DT, implemented under the same conditions; and (iii) investigating land use and land cover transitions for two decades, i.e., between 2010 and 2023, in a key agro-industrial corridor in Morocco. Thus, the approach presents not only the most up-to-date and consistently mapped landscape of Berrechid–Settat but also an adaptable methodology for studying other semi-arid regions grappling with similar environmental challenges and urban sprawl (Chomani and Pshdari 2024a).

The primary goal of our study is to categorise changes in LULC and assess the effectiveness of three machine learning algorithms. The classification process involved five LULC categories: vegetation, built-up areas, forests, barren land, and water bodies. The efficacy of these models was evaluated using accuracy metrics such as Overall Accuracy (OA), User’s Accuracy (UA), Producer’s Accuracy (PA), and Kappa Coefficient (KC). Our findings indicate that the DT model performed better in 2010, achieving an overall accuracy of 87.36% and a Kappa coefficient of 0.79, highlighting a strong agreement between the reference and classified data. By 2015, the RF model had emerged as the top performer with an overall accuracy of 88.74% and a Kappa coefficient of 0.81, while DT continued to perform strongly. In 2023, RF remained the best-performing model, with an OA of 91.84% and a KC of 0.86, underscoring its robustness in classifying LULC changes over time. The objective of this study is to enhance the existing knowledge on LULC classification and offer valuable insights for environmental management and policy development in the Berrechid-Settat area. Furthermore, it showcases the efficacy of GEE in efficiently obtaining and manipulating vast quantities of satellite imagery data, providing vital insights for land planning and management in the region, fostering sustainable development (Achahboun et al. 2023; Laalaoui et al. 2024; Sellami and Rhinane 2023).

The Berrechid Settat area, situated in the Casablanca-Settat region of Morocco (Figure 1), has been chosen as the study location due to its diverse landscape (El Assaoui 2023; Elgendy and Abdelatif 2021). Berrechid Settat lies between the Casablanca region to the north and the Marrakech-Safi region to the south (Bouzekraoui et al. 2020).

The area is situated between longitudes 7° 43’ 48’ W and 7° 10’ 48’ W, and the latitudes of 32° 57’ 36’ N and 33° 25’ 12’ N. The total land area measures 2062 square kilometres. It has varied topography, including agricultural plains, urbanised regions, and forests (Koukal et al. 2020).

The Berrechid Settat region is mainly known for its large agricultural areas, where a considerable amount of land is used for growing crops. It also includes forested areas with different amounts of tree cover, from dense forests to lighter woodlands. The area is experiencing significant urban expansion, with much development taking place in the areas on the outskirts of cities (El Assaoui Fouad Amraoui and El Mansouri 2015; Laalaoui et al. 2024).

In this work, we employed Landsat satellite images, specifically utilising Landsat 7 for the year 2010 and Landsat 8 OLI for the years 2015 and 2023.

The Landsat satellites provide multispectral imagery with various resolutions (Pickens et al. 2020). For our analysis, we used the blue, green, red, near-infrared (NIR), and short-wave infrared (SWIR) bands (D. Chen et al. 2018). Additionally, several spectral indices were computed to improve classification accuracy.

All Landsat imagery was preprocessed in Google Earth Engine to ensure radiometric and temporal uniformity for the years 2010, 2015, and 2023. Specifically, for Landsat 8 OLI (2015, 2023) and Landsat 7 ETM+ (2010), we used the Collection 2 Surface Reflectance (SR) tier 1 products. Atmospheric correction was applied using LaSRC for Landsat 8 and LEDAPS for Landsat 7. Furthermore, masks were utilised to remove clouds, cloud shadows, dilated clouds, cirrus, and snow, based on the QA_PIXEL bitmask, and saturated pixels according to QA_RADSAT, retaining only clear observations. Additionally, the optical reflectance bands were adjusted using the Collection-2 SR scale and offset. Subsequently, these corrected images were clipped to the relevant year (from 1 January to 31 December) and Area of Interest (AOI). Composites were built using the median to generate a cloud-free image for this period, reducing residual cloud cover and avoiding Landsat-7 SLC-off striping artefacts. This cloud-free imagery dataset is spectrally consistent between Landsat sensors and across dates, serving as the input for LULC. Table 1 summarises the data, including the imagery acquisition and bands used (Achahboun et al. 2023; Bhungeni et al. 2024; Wachowska et al. 2018).

For each year, a composite image was generated by processing the entire Landsat dataset. Using Google Earth Engine (GEE), the median value for each pixel over the selected time period was computed. The code defines the Area of Interest (AOI) and retrieves Landsat image collections for the specified years. Filters are applied to match the AOI and relevant dates, and the appropriate bands are selected for analysis. The composite images are created by calculating the median pixel value for each year and clipping the result to the AOI. These images are then displayed on the map using GEE’s visualisation tools (Sellami and Rhinane 2023; Z. Zhao et al. 2024).

In the classification process, to promote class separability, multiple spectral indices were derived from the ATM-corrected Landsat datasets. Each index was selected to highlight a specific surface characteristic relevant to the land use and land cover (LULC) patterns of the Berrechid–Settat area. The Normalised Difference Vegetation Index (NDVI) was applied to highlight vegetated areas. Near-infrared (NIR) represents band 4 of LANDSAT 7 (L7) and band 5 of LANDSAT 8 (L8), while Red denotes band 3 from L7 and band 4 from L8. The Normalised Difference Built-up Index (NDBI) was used to identify urban and impervious patch areas by contrasting short-wave infrared and near-infrared bands B5 and B4 for L7, and B5 and B6 for L8. Water bodies were enhanced using the Normalised Difference Water Index (NDWI), which incorporates the green band B2 for L7 and B3 for L8, and the near-infrared band B4 for L7 and B5 for L8. The Bare Soil Index (BSI) (equation 4) was computed to determine bare and sparsely vegetated areas by integrating reflectance information from the blue, red, NIR, and SWIR spectral bands B1, B3, B4, B5 from L7, and B2, B4, B5, B6 for L8 (Harfouche et al. 2020). All indices were computed within the GEE environment and included as additional predictor layers to the reflectance composites, thus increasing the spectral feature space used by the machine-learning algorithms SVM, RF, and DT for land use and land cover classification (Aydin and Sefercik 2025; Dong et al. 2019; Harfouche et al. 2020).

The processed composite images were subsequently used for further analysis and classification with machine learning models (Du et al. 2018; Lee et al. 2017). This allowed for the detection of changes in land use and land cover (LULC) within the Berrechid-Settat area across the selected years.

After choosing the set of Landsat 7 images for the year 2010, as well as Landsat 8 OLI images for the years 2015 and 2023, the selected images underwent a cloud-masking process to guarantee the clarity and accuracy of the data. Subsequently, five datasets were created, comprising the blue (B2), green (B3), red (B4), near-infrared (B5), and short-wave infrared (B6) bands (J. Chen et al. 2017). In order to train and validate our data, sampling and validation were used. Independent reference data were collected using a stratified random design with a minimum of approximately 60 points per class (total n = 300) interpreted from high-resolution imagery and local knowledge. Validation points were withheld from training and enforced a ≥300 m buffer from training AOIs to limit spatial autocorrelation. These points were distributed proportionally to mapped class area (Belgiu and Drăguţ, 2016). Accuracy was quantified using confusion matrices to derive OA, PA, UA, and Kappa, with 95% confidence intervals computed via bootstrap under a stratified estimator. This sample size provides stable per-class estimates for the five LULC categories: Vegetation, Built-up, Forest, Bareland, and Water Body.

The models were then examined by integrating the outcomes with testing data via accuracy assessments (Amin et al. 2024). Ultimately, the LULC maps were generated for the years 2010, 2015, and 2023. Figure 2 illustrates the Google Earth Engine (GEE) platform (Harfouche et al. 2020).

To conduct a pixel-based supervised classification, a distinct collection of training samples was acquired for each year (2010, 2015, and 2023) (Conrad et al. 2020). The assignment of each training sample pixel to a LULC class was determined using additional data sources, such as high-resolution photography and pre-existing LULC maps. The categorisation was executed with Support Vector Machine (SVM), Random Forest (RF), and Decision Tree (DT) classifiers within the Google Earth Engine platform (El Assaoui et al. 2023; Kamusoko and Gamba 2016).

The three selected algorithms have been chosen for their reliable performance and complementary strengths in remote sensing. The SVM classifier is particularly recognised for its resilience to high-dimensional spectral data and few training samples, making it appropriate for a heterogeneous environment. The RF has ensured liability against overfitting and offers excellent performance on multi-class tasks through ensemble learning. The DT, while the most straightforward method, does not suffer from the black box problem and is computationally accelerated. Therefore, it provides a perfect error metric with more complicated ensemble and kernel methods, offering interpretability and computational efficiency. Their combination enables a comprehensive evaluation of classification performance under differing algorithmic assumptions (Chomani and Pshdari 2024b; Feng et al. 2015).

The three classifiers used in this study for LULC categorisation are detailed below:

A Support Vector Machine (SVM) is a type of supervised learning algorithm used for both regression and classification. During the training phase, SVM classifiers construct an optimal hyperplane that effectively divides classes by minimising the misclassification of pixels from input datasets. Crucial factors for Support Vector Machines (SVM) include kernel functions, cost parameters, and gamma (Audebert et al. 2018).

Random Forest (RF) is an ensemble classifier that builds multiple decision trees using random subsets of the training data and attributes. The key input parameters for RF are the size of the training dataset and the number of trees created. In this study, 300 trees were used to improve the accuracy and performance of the classification.

A Decision Tree (DT) is a classification system that uses specified thresholds to make decisions. It splits nodes until it reaches terminal nodes. This approach involves grouping input data into distinct sets and then generating trees for the purpose of classification (Akar and Gormus 2021).

The LULC maps were produced entirely using Python-based workflows. Raster and vector datasets were imported and processed using libraries such as rasterio (v1.3.8), geopandas (v0.14.3), and numpy (v1.26.4). These tools facilitated image stacking, band alignment, and clipping of imagery to the Berrechid–Settat study area, preparing the data for classification. Only essential preprocessing tasks, such as cloud masking to remove pixels affected by cloud cover, were performed in Google Earth Engine (GEE). We used temporal aggregation methods, such as mean and median calculations, to fill in gaps left by cloud interference (Carlson et al. 2018; Wu 2020).

The LULC maps were produced entirely using Python-based workflows (He et al. 2020). High-resolution orthophoto images and pre-existing LULC maps were used to collect training and validation samples. These samples were saved as shapefiles and then imported into the Python environment to train the classifiers. The LULC classification in this study included five primary categories: forest, built-up land, barren land, bodies of water, and vegetation (Sharnagat et al. 2025).

For randomness in stratified sampling without attrition, maintaining the proportional representation of all LULC classes in both subsets was crucial for class balance. To reduce spatial autocorrelation and avoid overfitting, training and validation points were dispersed geographically within the AOI. A suitable minimum distance was also implemented between samples of the same class (Riche et al. 2024). This methodology enhanced the independence of the test data, which facilitated the model’s generalisation. The samples were divided using the same ratios. Before training the model, min-max normalisation was applied to the spectral bands and resulting indices. This brought all features to a uniform scale of 0–1. This was achieved using the StandardScaler and MinMaxScaler functions from the scikit-learn.preprocessing package (v1.4.2) (Saqr et al. 2025). This standardisation helped to lessen the variations in range across the reflectance of different bands. It assisted the algorithms in maintaining numerical stability, particularly in distance-based models such as SVM and RF. This normalisation ensured that each feature contributed equally to classification, thereby improving the performance of all considered models through convergence.

After preparing the training and testing data, the models were run with default settings. First, supervised classification was implemented using the Scikit-learn library (v1.4.2). Parameters were specified for the SVM model, such as the kernel type (e.g., ‘rbf’). After initialising the model, the ‘fit()’ function was used to train it on the training data. Predictions were then evaluated using the test data (Phan et al. 2020). The accuracy score function from Scikit-learn assesses the model’s accuracy for the predictions. The datasets used in this study were generated through stratified random sampling, with 70% allocated for training and 30% for validation et al. 2020). The accuracy score function from Scikit-learn assesses the model’s accuracy for the predictions. The datasets used in this study were generated through stratified random sampling, with 70% allocated for training and 30% for validation (Chomani and Pshdari 2024b; Lu et al. 2021).

Evaluation of classification accuracy was conducted by measuring overall accuracy (OA) and the Kappa coefficient (K), producer’s accuracy (PA), and user’s accuracy (UA). Their formulas are detailed below in equations 1, 2, 3, and 4 (Kobayashi et al. 2022). These metrics assess the quality of classification by quantifying the accuracy of the test data and the agreement between classified and reference data. Their calculations were performed using scikit-learn.metrics (V1.4.2), with functions such as confusion_matrix, classification_report, and cohen_kappa_score. Results were displayed and interpreted using matplotlib (V3.8.4) combined with seaborn (V0.13.2) for graphical display, statistics, and confusion matrices (Figure 5).

(1)

(2)

(3)

(4)

Where:

Kappa values range between 0 and 1, where:

The final step involved using the trained models to generate predictions and exporting the results as GeoTIFF files (Dong et al. 2019). While various tools such as GDAL or GeoPandas/Geocube are available, they are not well-suited for large-scale geographic operations (Drusch et al. 2017). An efficient approach using the rasterio library is presented below (Gerber et al. 2018).

First, the original input image is retrieved, and metadata attributes such as height, width, and CRS are extracted. The input image is then preprocessed to align with the data used during model training. If additional features like NDVI or elevation were included during training, they must be combined with the input image before proceeding. This process ensures scalability for generating predictions over large geographic areas (Jiang et al. 2021). The results are saved in GeoTIFF format for further analysis.

Once the model training and validation were complete, each machine learning classifier was applied to the entire Landsat imagery to predict the LULC category for each year of the study: 2010, 2015, and 2023. The prediction process was then performed using the ‘predict()’ function from the scikit-learn library (v1.4.2), which resulted in pixel-wise classification outcomes. Finally, the resulting arrays were restructured and converted to raster formats using rasterio (v1.3.8) to ensure the geospatial conformity of the final outputs with the imagery. The classified rasters were then exported as GeoTIFF files for further processing and map production (Peterson and Levinson 2020). High-quality theme map production and post-processing were carried out using matplotlib (v3.8.4), geopandas (v0.14.3), and rasterio.plot modules. A unique colour scheme was assigned to each land category, including vegetation, urban, forest, barren land, and water body (Gerber et al. 2018; Wu 2020). Additional map elements, such as the title, legend, north arrow, and scale bar, were included to improve map interpretation. All the final maps were exported at a resolution of 600 dpi, which is suitable for publishing and scientific sharing. This process helped to create accurate geo-referenced outputs that showed the precise spatial distribution and temporal changes of LULC in the area of interest (Sharnagat et al. 2025; Q. Zhao et al. 2021).

Fig. 3 shows the LULC classification maps produced by three machine learning algorithms: SVM, RF, and DT, over the 13 years of the study. For each classifier, these colour-coded images display the spatial distribution pattern of five LULC classes: vegetation, built-up, forest, bareland, and water body. SVM’s classification exhibits much smoother and more homogeneous class boundaries, whereas RF and DT show local heterogeneity. In contrast, the green colour representing bareland varies more across the years, indicating the classifiers’ sensitivity to spectral variation. These multi-model maps offer a suitable platform for visual comparison of spatial patterns and temporal dynamics, which highlighted the observation of continuous urban expansion and a slight reduction in vegetation cover between 2010 and 2023.

The classified area distributions, as shown in Figure 4 and Table 3, reveal distinct temporal shifts and model-based fluctuations in the AOI over 13 years of the study. Vegetation, covering approximately 1200 to 1400 km² across all classifiers, continues to dominate this class. However, it decreased slightly by 2023, while built-up areas increased from about 200 km² in 2010 to over 300 km² in 2023. This data aligns well with the area’s ongoing urban expansion and agricultural land modification (Du et al. 2018). The forest class maintained intermediate values between 250–350 km² with little temporal change across models and years. Bareland, however, showed a gradual reduction from 150–200 km² to less than 100 km² in 2023. Water bodies represented the smallest fraction, less than 70 km², with a barely noticeable reduction over time (Carrara et al. 2024; Congalton and Green 2019). DT classifiers produced slightly higher vegetation estimates in 2010 compared to other classifiers. RF outputs were more balanced and consistent in the following years, particularly for built-up and bareland classes. The SVM model performed between these two methods, yielding slightly lower vegetation and higher built-up estimates due to its effectiveness in distinguishing mixed spectral signatures (Ali and Johnson 2022). Generally, the three algorithms reported similar temporal trends, which increases confidence in our classification methodology. RF is, however, the most suitable method for multi-temporal LULC mapping in this region because of its stable intra-class ratio and consistent spatial predictions (Meghraoui et al. 2024).

Each model operated on a confusion matrix derived from 30% of the entire figure 5. A confusion matrix is a statistical construct that illustrates the correlation between the actual class (ground truth) and the predicted or generated map (Nicolau et al. 2023). The confusion matrices derive overall accuracy (OA), Kappa coefficient, user’s accuracy (UA), and producer’s accuracy (PA) by land use/land cover type (Belgiu and Csillik 2018). Therefore, the confusion matrices yield the following OA and Kappa coefficients for the classification trees. The SVM confusion matrix shows an OA of 91.7% and a Kappa coefficient of 0.810. The RF confusion matrix shows an OA of 91.84% and a Kappa of 0.86. The DT confusion matrix shows an OA of 91.16% and a Kappa of 0.68. A Kappa of 0.845 was recorded for 2023. Table 4 and Figure 5 show the overall model performance over the years and classifications. However, the RF model had the highest overall OA and Kappa from all datasets. The classes performed across the models. The built-up land use class was achieved in all models with UA and PA exceeding 97% for 2023 RF. SVM and DT also had UA and PA exceeding 89%. On the other hand, the forest class performed less well with SVM, showing a PA as low as 75.2% in 2023 (Wang et al. 2018). Similarly, the numerical classifications for bare land and waterbody. The only inconsistency is the 81% PA for forest classification in 2020, which is still a decent score but 6% lower than the DT for forest classification accuracy. Where the forest category was not accurately represented was in built-up areas. This is similar to how the forest was perfectly classified in 2020 across all three models. However, all three models struggled with this over the years, with SVM again showing the highest omission errors (El Assaoui et al. 2015; J. Jiang et al. 2025). This resembles a replica of a city covered up by what should have been classified as non-built-up pixel values. Where SVM struggled, the RF model excelled, achieving the best classification for built-up in 2023 at 88.4%. The lowered classification accuracy for some classes in 2023 means that for bareland and waterbody classes, DT’s PA values were lower than RF and SVM. However, it still achieved fairly good accuracy for built-up and vegetation classes (Congalton and Green 2019).

The study has illustrated the effectiveness of various machine learning algorithms, including SVM, RF, and DT, in categorising land use and land cover (LULC) changes in the Berrechid Settat region for the years 2010, 2015, and 2023 (El Assaoui et al. 2021; Kussul et al. 2017). A comprehensive assessment of accuracy measurements, such as Overall Accuracy (OA), Kappa Coefficient, User’s Accuracy (UA), and Producer’s Accuracy (PA), shows that each approach has distinct advantages and disadvantages when used for different LULC classes.

Random Forest (RF) consistently achieved the highest accuracy throughout all years, particularly excelling in the classification of built-up, water body, and barren lands categories (Fentaw and Abegaz 2024), attaining elevated user’s accuracy (UA) and producer’s accuracy (PA) values. This indicates that RF is an effective model for land use and land cover (LULC) mapping, especially in intricate systems characterised by significant spatial variability. Conversely, the Support Vector Machine (SVM) exhibited commendable performance. Nevertheless, it faced difficulties in distinguishing certain classes, such as forest and bare ground (Baatz et al. 2025), resulting in more omission errors for these categories. The Decision Tree (DT) displayed similar accuracy across many classes but showed a slight decrease in performance, especially for smaller or more heterogeneous categories such as water body and barren land.

A notable discovery is that all models successfully distinguished the built-up class, which is essential for tracking urban expansion and land development in swiftly urbanising areas such as Berrechid-Settat. The forest class presented the most significant obstacle among all models, highlighting the necessity for enhanced spectral feature selection or the integration of supplementary data to improve forest classification accuracy (Andrew et al. 2023; Drusch et al. 2017).

In summary, RF has demonstrated the highest reliability for LULC classification in this area, especially for classes with distinct spectral signatures, such as built-up areas and aquatic bodies (Dharumarajan and Hegde, 2022). SVM and DT offered significant insights, particularly regarding processing efficiency and the management of intricate decision boundaries. This study highlights the importance of model selection in remote sensing applications and provides a foundation for future research focused on improving LULC classification accuracy through more sophisticated hybrid models or by including other data sources such as topography and climate factors (Alshari and Gawali 2021; Amazirh et al. 2024; El Assaoui et al. 2023).

This analysis offers a thorough understanding of LULC dynamics in the region, providing significant insights for land management, urban planning, and environmental monitoring.

Our functional framework for land monitoring could be useful in regional urban planning and water management initiatives. However, the classification accuracy could be improved by using Sentinel-2 imagery or higher-resolution data, or by incorporating an object-based approach to the research. Furthermore, future research should focus on including socio-economic and hydro-climatic variables to better explain the drivers of LULC changes.