Predicting land use and land cover change dynamics in the eThekwini Municipality: a machine learning approach with Landsat imagery

ABSTRACT

Monitoring and providing accurate land use and land cover (LULC) change information is vital for sustainable environmental planning. This study used Landsat imagery from 2002 to 2022 to create updated LULC change maps for the eThekwini Municipality. Random Forest (RF), Support Vector Machine (SVM), and Extreme Gradient Boosting (XGBoost) were used to conduct these LULC classifications, with XGBoost achieving the highest accuracy (80.57%). The generated maps revealed a significant decrease in cropland and an increase in impervious surfaces. As such, this research established a framework for continuous LULC mapping and highlighted Landsat 9’s potential in LULC classifications.

KEYWORDS:

• Land cover
• land use
• remote sensing
• machine learning
• Landsat

1. Introduction

The complex environmental, social, economic, and political matrix is the primary driver of a territory’s spatial arrangement (Trujillo-Jiménez et al. 2022). This matrix has also accelerated the natural Earth’s surface changes through increased resource consumption to compensate for increased human populations and accelerated economic development (Kafy et al. 2021, Digra et al. 2022). Land use and land cover (LULC) maps serve as excellent depictions of this matrix (Trujillo-Jiménez et al. 2022), primarily because producing up-to-date and accurate LULC change information provides an understanding of human-environment relationships and the adverse impacts of global change (Chakraborty et al. 2016, Nasiri et al. 2022). Moreover, this updated information also provides crucial insights for implementing sustainable land management practices, conserving natural resources, and enhancing resilience to climate-related hazards (Effiong et al. 2024). Therefore, LULC change information, especially at the regional scale, aids decision-makers in mitigating the persisting challenges of global change that results from LULC change through appropriate environmental management policies and strategies (Ai et al. 2020, Tadese et al. 2020). Thus, LULC maps are the key to sustainable environmental planning and development (González-González et al. 2022). This primarily pertains to impoverished rural regions, some of which are key areas for ecological and biodiversity significance (Musetsho et al. 2021).

The eThekwini Municipality has attempted to protect and manage its biodiversity and ecological hotspots that supply ecosystem services crucial for human well-being by establishing the Durban Metropolitan Open Space System (DMOSS) (Boon et al. 2016). This was due to the increased degradation and loss of the vegetated landscapes within the municipality (Zungu et al. 2020). However, the eThekwini Municipality is still susceptible to adverse environmental consequences due to the change driven by the municipality’s attempts to move past the apartheid’s discriminative spatial planning (Sutherland et al. 2014). Additionally, there is a lack of comprehensive knowledge and understanding regarding the rate of LULC change within the municipality and its ramifications for ecosystems and ecosystem services. Therefore, by providing updated LULC change maps, the eThekwini Municipality can implement appropriate environmental or spatial management policies and strategies. These maps will also inform an update to the DMOSS layer, which was last updated in 2018.

Generating LULC change maps previously relied on in situ monitoring techniques (Buchhorn et al. 2020). However, the prominence of remote sensing over the past decades, and the improved accessibility and development of various methods for extracting remotely sensed data has seen more LULC mapping studies utilizing remote sensing (Amini et al. 2022). Remote sensing presents the ability to effectively capture LULC change information (Nkundabose 2021) by making it possible to integrate spectral-temporal metrics and a diverse set of features to capture changes in various LULC classes (Liang et al. 2022). Most studies have utilized medium-resolution sensors such as Landsat (30 meters – m) and Sentinel 2 (10 m) given that these sensors can detect most human-nature interactions (Chen et al. 2015). These sensors also provide a short revisit time and vast spectral configurations (Nasiri et al. 2022). Hence, this research utilized imagery from Landsat 7 (Enhanced Thematic Mapper – ETM+), Landsat 8 (Operational Land Imager 1 - OLI 1), and Landsat 9 (Operational Land Imager 2 - OLI 2) with a spatial resolution of 30 m to analyze the changes in LULC within the eThekwini Municipality. The selection of these three sensors is also based on the consistency of the Landsat data record across the entire Landsat sensor series (Roy et al. 2016).

The utilization of remote sensing in mapping LULC change has attracted the application of machine and deep learning approaches. These approaches are applied to reduce spectral and spatial limitations associated with medium and low-resolution LULC classifications, which should result in more precise LULC maps (Talukdar et al. 2020). Some of the most popular algorithms in LULC change classification include Artificial Neural Networks (ANN), random forests (RF), Cellular Automata (CA), support vector machines (SVM), Markov Chain, and extreme gradient boosting (XGBoost) (Abdi 2020, Kafy et al. 2021). However, RF and SVM tend to produce more accurate LULC maps compared to the algorithms aforementioned (Abdi 2020, Talukdar et al. 2020, Ali et al. 2022).

Therefore, this research investigated the effectiveness of three machine learning algorithms, specifically, SVM, RF, and XGBoost in mapping changes in LULC. However, it should be noted that remote sensing based LULC change classifications can be limited by the requirement to access and store large remote sensing data and high computational costs. Platforms, such as the Google Earth Engine (GEE), which is cloud-based, have overcome these limitations (Tassi and Vizzari 2020). GEE allows scientists and independent researchers to access large data banks for change detection freely and quantify resources on the Earth’s surface (Mutanga and Kumar 2019). Therefore, Landsat data were obtained and processed in GEE.

Overall, the objective of the study was to map LULC change from 2002 to 2022 within the eThekwini Municipality using the three most recent Landsat sensors and machine learning algorithms. This will aid in assessing the present rate of LULC change and identifying the factors responsible for these changes. The absence of comprehensive knowledge on the rate of LULC change and the lack of understanding of the consequences of LULC change for ecosystems and ecosystem services within the eThekwini Municipality highlight the critical need for updated LULC change maps to inform decision-making and policy development. This research addresses these gaps by providing enhanced accuracy and efficiency in LULC mapping by employing machine learning algorithms. Such advanced methodologies bring significant advancement to the fields of remote sensing and spatial science, enabling more precise monitoring and analysis of land use and cover dynamics. Furthermore, including Landsat 9 imagery expands the temporal scope of the study. It also enhances the accuracy of LULC change detection, emphasizing its importance in improving the understanding of landscape dynamics over time. Consequently, the generated maps will not only facilitate the identification of areas for vegetation restoration within the eThekwini Municipality but also contribute to developing appropriate environmental or spatial management policies and strategies. Thus, the study’s findings will provide policymakers with invaluable insights to make decisions that are well-aligned with practical realities and sustainability goals.

2. Materials and methods

2.1. Study area

The eThekwini Municipality, established in 2000, is located in the South African province of KwaZulu-Natal (Figure 1), along the Indian Ocean’s eastern coast (Musvoto et al. 2016, Shivambu et al. 2020). The municipality serves as a pivotal economic hub in South Africa, largely due to its renowned port which is recognized as the busiest port, on the African continent (Zungu et al. 2020). The municipality is formed by 103 wards with an approximate population of 3.5 million, which grows by one percent annually (Zungu et al. 2019). Moreover, its mainland predominantly comprises urbanized zones, constituting the city of Durban, which is surrounded by peri-urban residential areas and informal settlements, along with rural communities (Hellberg 2014). With the population and rural-urban migration increasing and more residential spaces being required, the municipality established the DMOSS, a network of areas reserved for conserving native fauna and flora (Khumalo and Sibanda 2019, Zungu et al. 2019). DMOSS also aimed to limit the conversion of native vegetation for anthropogenic purposes, which include agriculture, residential areas, and roads. As of 2019, 53% of native vegetation had already been converted for human benefit (Zungu et al. 2020).

Figure 1. The eThekwini Municipality’s location in KwaZulu-Natal, South Africa.

2.2. Data acquisition and pre-processing

2.2.1. Reference data

The 2017 Finer Resolution Observation and Monitoring-Global Land Cover-Segmentation at 10 m (FROM – GLCS10) was used as a reference. FROM – GLCS10 offers a fundamental land cover scheme that incorporates continuous field layers for all primary land cover classes. These classes provide proportionate estimates of vegetation/ground cover for the respective land cover types. Classes in land-cover products generally utilize a hierarchical structure and, thus, this is also the case in FROM – GLCS10 (Zhu et al. 2021). The classes in FROM – GLCS10 found within the eThekwini Municipality are presented in Table 1. The LULC class descriptions are derived from the descriptions provided by Wei et al. (2020) and Buthelezi et al. (2024).

Table 1. From – GLCS10 land cover classes (Wei et al. 2020, Buthelezi et al. 2024).

Name	Description	Code
Cropland	Land under cultivation showing discernible patterns on satellite imagery.	1
Forest	Areas designated for commercial plantations, orchards, and native tree species.	2
Grassland	Areas that comprise open landscapes dominated by various species of grasses.	3
Shrubland	Land characterised by low-lying vegetation dominated by shrubs.	4
Wetland	Areas characterised by saturated or seasonally saturated soils	5
Water	Areas characterised by the presence of rivers, lakes, ponds, and reservoirs.	6
Tundra	Treeless landscapes in cold regions with vegetation adapted to harsh conditions, including mosses and lichens.	7
Impervious surface	Encompasses developed areas comprising both urban and rural residential zones as well as industrial sectors.	8
Bare land	Land devoid of vegetation cover.	9
Ice/Snow	Areas covered by frozen water, such as glaciers, ice caps, and ice shelves	10

2.2.2. Landsat imagery

Landsat 7, Landsat 8, and Landsat 9 were used to map LULC change within the eThekwini Municipality. This enabled an extensive assessment of factors that contributed to spatial changes within the municipality over 20 years (2002 to 2022). Landsat 7 (ETM+) was sent into orbit on 15 April 1999 with a fixed ‘whisk-broom’ multispectral scanning radiometer with eight bands (Moradpour et al. 2020). The ETM+ sensor can observe the Earth’s surface over a span of 185 km in width and operates on a revisit cycle of 16 days (Cao et al. 2022). However, on 31 May 2003, the scan-line corrector (SLC) of the ETM+ sensor experienced a permanent malfunction, leading to a situation where ETM+ images only captured 78% of the intended scanned pixels (Wang et al. 2021). Unscanned pixels in ETM+ images appear as a linear gap, measuring approximately 12 pixels in width. Therefore, these line gaps must be filled before Landsat 7 images can be used in any application.

Landsat 7 was followed by Landsat 8, which is equipped with the OLI and the Thermal Infrared Sensor (TIRS), providing a spatial resolution of 30 meters and a 16-day revisit cycle (Gorelick et al. 2017). OLI is a push-broom imager with nine spectral bands, and TIRS is also a push-broom imager with two bands (Vanhellemont 2020). The OLI sensor exhibits enhanced calibration and signal-to-noise properties, narrower spectral bands, a superior 12-bit radiometric resolution, and more accurate geometry in comparison to the ETM+ sensor (Roy et al. 2016). Landsat 8 was followed by Landsat 9, launched on the 27^th of September 2021. Landsat 9 has two instruments on board, namely, the OLI-2 and the Thermal Infrared Sensor 2 (TIRS-2). These two instruments improve on similar instruments found onboard Landsat 8 (Showstack 2022). The key distinction between OLI and OLI-2 lies in the data download process: OLI-2 downloads all 14 bits for each image pixel, whereas OLI only downloads 12 bits (Gross et al. 2022).

Furthermore, Landsat 7, Landsat 8, and Landsat 9 images for the study years were acquired through the GEE platform, utilizing the Landsat 7 Collection 2 Tier 1 calibrated top-of-atmosphere (TOA) reflectance, Landsat 8 Collection 2 Tier 1 TOA reflectance, and Landsat 9 Collection 2 Tier 1 calibrated TOA reflectance, respectively, as specified. The calculation of TOA reflectance values for the datasets employed in this research is found on the methodology established by Chander et al. (2009). These collections were filtered by the eThekwini Municipality boundary (as a shapefile imported to GEE) and cloud cover (<1%). The eThekwini Municipality boundary covered three Landsat scenes; thus, the filtering process was made uniform for all scenes, which were further clipped using the study area boundary. The clipped images were then mosaiced.

The initial dataset for classifications was composed of bands 1 to 7 from the mosaiced Landsat 7 images and bands 2 to 8 from the mosaiced Landsat 8 and Landsat 9 images. The reason for excluding other bands is that these bands do not have information about the surface; instead, they best characterize the atmosphere (Li et al. 2020). Additionally, to improve the LULC classification, five spectral indices, namely, the Enhanced Vegetation Index (EVI), Normalised Difference Vegetation Index (NDVI), Normalised Difference Water Index (NDWI), Modified Difference Water Index (MNDWI) and the Normalised Difference Built-up Index (NDBI) were calculated in GEE for all Landsat images and added to the data for classification (Table 1).

However, for Landsat 7, before spectral indices were calculated in GEE for the 2007 and 2012 mosaiced images, they still needed to be processed further to fill the unscanned pixels. Hence, the Geo-statistical Neighbourhood Similar Pixel Interpolator (GNSPI) was employed to address the unscanned pixels, with a detailed explanation provided by Zhu et al. (2012). The GNSPI code, scripted in the Interactive Data Language, is accessible online. GNSPI utilizes empirical or physical models and ordinary Kriging to predict the data lost due to SLC issues (Yin et al. 2017). After the 2007 and 2012 images, unscanned pixels were filled, these images were ingested as assets back to GEE, and the five spectral indices used in this study were calculated (Table 2).

Table 2. The spectral indices that were calculated for inclusion in the training data.

Index	Equation	Label	Reference
NDVI	$\frac{\left(\mathit{NIR}-\mathit{RED}\right)}{\left(\mathit{NIR}+\mathit{RED}\right)}$	1	Tucker and Sellers (1986)
EVI	$\frac{[\left(2.5\mathrm{X}\left(\mathit{NIR}-\mathit{RED}\right)\right]}{[\left(\mathit{NIR}+\left(6\mathrm{X}RED\right)-\left(7.5\mathit{XBLUE}\right)+1\right]}$	2	Huete et al. (2002)
NDWI	$\frac{\left(\mathit{NIR}-\mathit{SWIR}\right)}{\left(\mathit{NIR}+\mathit{SWIR}\right)}$	3	Gao (1996)
MNDWI	$\frac{\left(\mathit{GREEN}-\mathit{SWIR}\right)}{\left(\mathit{GREEN}+\mathit{SWIR}\right)}$	4	Xu (2006)
NDBI	$\frac{\left(\mathit{SWIR}-\mathit{NIR}\right)}{\left(\mathit{SWIR}+\mathit{NIR}\right)}$	5	Zha et al. (2003)

EVI and NDVI are useful indices for classifying vegetation; however, it should be noted that NDVI saturates when classifying densely vegetated areas (Mutanga et al. 2012, Buthelezi et al. 2020). The shortcomings of NDVI are compensated by EVI, which performs better when classifying areas with high biomass and reduced atmospheric influence (Qiu et al. 2018). NDWI and MNDWI are best suited for identifying water bodies (Raut et al. 2020, Shahfahad et al. 2022). NDBI has been revised multiple times and is critical for improving the built-up area classification (Santra et al. 2021). It achieves this by separating built-up areas from barren land (He et al. 2010).

Qu et al. (2021) proposed using different auxiliary features to improve the accuracy of LULC classification that utilizes remote sensing images with medium resolution. Therefore, this study incorporated SRTM Digital Elevation Data at 30 m resolution.

2.3. Sample data composition

To compose the training dataset, the seven bands from the 2017 Landsat 8 image, computed spectral indices, elevation data, and the 2017 FROM – GLCS10 image were stacked in Google GEE, resulting in a single multiband image. Stratified sampling, where a maximum of 10,000 points was set for each 2017 FROM – GLCS10 class, was then carried out on the multiband image containing the information for the classification. The sample data were exported to the R computing platform (version 4.3.2) through the RStudio integrated development environment (IDE) (version 2022.07.2 + 576) in a comma-separated values (CSV) file. It should be noted that all R computations in this study were conducted in RStudio. However, before the models were trained in RStudio, the data were divided into 70% training and 30% test data. Both training and test data were then resampled using the K-fold cross-validation technique, which is part of the caret package in RStudio. Given the large size of the dataset, K = 5 was used and repeated ten times.

2.4. Machine learning algorithms

2.4.1. Random forests

RF was adopted to conduct the LULC classification using Landsat 7, 8, and 9 images. RF is a group of regression trees that grow from the random selection of training data samples (Ali et al. 2012), as outlined in Figure 2. Ma et al. (2017) and Masiliūnas et al. (2021) posited that RF, which improves from the CART algorithm, performs well when classifying remotely sensed imagery. RF can also solve classification and regression problems (Nghia et al. 2021). Thus, it has become the most utilized machine learning algorithm in classifying remote sensing images (Luo et al. 2021).

Figure 2. The formation of RF, as outlined in Chen et al. (2021).

2.4.2. Support vector machines

SVM is one of the most robust pattern recognition methods, which operates by determining an optimal separation line (hyperplane) that accurately separates two or more classes (Bennett-Lenane et al. 2022). Determining the optimal hyperplane maximizes the generalization ability of the SVM model (Cervantes et al. 2020). However, if the data is not linearly separable as in this study, the generalization of the model is hindered. Thus, a kernel function is employed (Saygılı 2022). This study selected the radial basis function (RBF) to perform the LULC classification because it is effective in handling non-linearly separable data as was the case in this study (Anyanwu et al. 2023).

2.4.3. Extreme gradient boosting

XGBoost is an extension of gradient boosting machines (GBMs), a widely employed technique in various remote sensing applications (Buthelezi et al. 2020). Its application has been limited to structural engineering, hence, its limited application in remote sensing (Nguyen et al. 2021). XGBoost model employs additive learning to develop a ‘strong’ learner through the integration of several ‘weak’ learners (Fan et al. 2021). The XGBoost algorithm was developed to a achieve low risk of overfitting and high accuracy (Nguyen et al. 2021).

2.5. Parameter optimisation

RF in RStudio requires parameter optimization, which significantly improves the algorithm’s accuracy (Buthelezi et al. 2020). The performance of RF is primarily influenced by two main parameters: the number of features randomly selected for constructing each tree (Mtry) and the total number of trees to be generated (Ntree). However, it should be noted that only Mtry can be adjusted manually. Therefore, to obtain the optimal Mtry, the expand.grid function in RStudio was used. In the case of SVM, the optimized parameters were the cost (C) and gamma (γ), both of which are independent of each other. A lower value of C places emphasis on maximizing the margin, whereas a higher value indicates a preference for minimizing errors in the SVM algorithm. The γ relates to how much curvature is required in a decision boundary. When the γ parameter is low, it means the single training sample is far-reaching, whereas when it is high, the opposite is true. Optimized parameters for RF and SVM are presented in Table 3.

Table 3. Optimized parameters for both RF and SVM.

RF Parameter	Value	SVM Parameter	Value
Mtry	1	Cost (C)	1
Ntree	1	Gamma (γ)	0.083

XGBoost requires multiple parameters to be optimized, among them; the number of trees, learning rate, and regularization parameters (Nguyen et al. 2021). The Bayesian optimization method, outlined in Bergstra et al. (2011), was utilized to obtain optimal parameters for XGBoost, and the resulting parameter values are presented in Table 4.

Table 4. XGBoost optimal parameters that were obtained using the Bayesian optimization method.

XGBoost Parameter	Value
Subsample Ratio (rs)	1
Learning Rate (τ)	0.3
Maximum Tree Depth (Dmax)	6
Minimum Child Weight (ωmc)	1
Column subsample ratio (rc)	0.9471

2.6. Variable importance

Variable importance for RF and XGBoost LULC classifications was determined using the importance and varImp functions in RStudio, respectively. RF’s variable importance can be presented using the permutation-based metric and a node impurity metric. The node impurity metric (mean decrease in Gini coefficient) was used for this study. It is defined as the average of a variable’s total decrease in node impurities after the splitting of the variable, averaged over all trees (Luo et al. 2021). For XGBoost, the relative importance (RI) score, which is also based on the Gini coefficient importance, was used (Kardani et al. 2022). RI for this study was scaled from 0 to 100. Given that SVM used RBF to perform the LULC classification, it was impossible to determine variable importance for the algorithm. That is because, in kernels other than the linear kernel, the separating plane exists in another space because of the kernel transformation of the original space.

2.7. Accuracy assessment

The primary evaluation of LULC classifications involved assessing overall accuracy, user’s accuracy, and producer’s accuracy derived from the confusion matrix, along with calculating of the area under the curve (AUC). The accuracy measures can be calculated in RStudio. The overall accuracy assesses the number of pixels correctly classified within the image and was employed as the primary metric in this study. The producer’s accuracy measures the precision with which reference pixels are classified, while the user’s accuracy assesses how accurately the classified map reflects the actual conditions on the ground (Patel and Kaushal 2010). The AUC is a measure of sensitivity and specificity over all possible threshold values (Halimu et al. 2019).

Furthermore, quantity and allocation disagreements were computed in RStudio and used to investigate the accuracy of the classifications further. Quantity disagreement measures the disparity in the quantity of a specific class of objects between the classified and reference maps, while allocation disagreement compares the spatial arrangement of pixels between the classified and reference maps (Pontius and Millones 2011, Fassnacht et al. 2014, Lottering et al. 2020, Verma et al. 2020). The allocation disagreement in this study uses the exchange and shift measures. Exchange occurs when a pair of pixels is classified as category A in the reference map and as category B in the classified map, and vice versa (Pontius and Santacruz 2014). Whereas shift is a measure of allocation disagreement used when comparing maps with more than two categories, which identifies situations where the spatial distribution of classes differs between the maps. Still, these differences are not due to a simple exchange of categories (e.g. forest classified as water in the reference map and water classified as forest in the classified map) (Bontempo et al. 2020).

2.8. Percentage area for each class

The LULC classification outputs were investigated in QGIS Desktop 3.26.1 for the quantity of the pixels belonging to each class using the raster layer unique values report tool. The values from the report were used to compute the percentages for the area occupied by each LULC class within the study area. The overall methodology is outlined in Figure 3.

Figure 3. The overall framework used to perform RF, SVM, and XGBoost LULC classifications.

3. Results

3.1. Variable importance

Figure 4 shows the variable importance for RF LULC classifications based on the 2017 Landsat 8 image classification. The figure indicated that the red and blue bands were the most important variables, whereas slope had the least importance. When comparing the importance of spectral indices used in RF classifications, NDVI was the most important index.

Figure 4. Variable importance based on the mean decrease in the Gini coefficient for RF LULC classification on the 2017 Landsat 8 image.

When XGBoost was used to perform LULC classification on the 2017 Landsat image, the red band and MNDWI were the most prioritized variables (Figure 5). Like in RF, the slope was the least prioritized variable.

Figure 5. Variable importance based on the RI for XGBoost LULC classification on the 2017 Landsat 8 image.

3.2. Accuracy assessment

3.2.1. Confusion matrix

All three algorithms successfully classified the 2017 Landsat 8 image in RStudio. However, these classifications were executed at varying accuracies (Figure 6). Based on the overall accuracy for RF, SVM and XGBoost LULC classifications performed on Landsat images, it was observed that XGBoost classifications (80.57%) performed slightly better than RF (80.42%) and SVM (79.59%). Based on the user’s accuracy RF and XGBoost (97.16%) outperformed SVM (97.05%). However, RF (92.27%) did perform fractionally better than XGBoost (92.26%) and SVM (91.75%) based on the AUC measure. RF (79.77%) also performed better than XGBoost (79.74%) and SVM (78.49%) based on the producer’s accuracy.

Figure 6. RF, SVM and XGBoost classification accuracies based on the overall accuracy, user’s accuracy, producer’s accuracy, and AUC measures from the 2017 Landsat 8 image.

3.2.2. Quantity and allocation disagreements

Figure 7 shows the results for disagreement metrics which represent the percentage of disagreement between LULC maps classified based on the 2017 Landsat 8 images and the reference map. The exchange and shift represent the percentage of allocation disagreement. RF LULC classification had a quantity disagreement of 15.72%, an exchange of 22.79% and a shift of 17.28% (Figure 7(a)). The quantity disagreement observed in the RF classification was significantly less than that observed in the SVM classification (Figure 7(b)), which was 70.66%. However, SVM performed better than RF when looking at the exchange and shift measures, which were 5.50% and 6.22%, respectively. Overall, XGBoost performed better than RF and SVM with a quantity disagreement of 11.92%, an exchange of 28.73% and a shift of 10.29% (Figure 7(c))

Figure 7. Disagreement metrics for the (a) RF, (b) SVM and (c) XGBoost LULC classification maps based on the 2017 Landsat 8 image.

Overall, XGBoost produced slightly more accurate classifications than RF and SVM. Hence, it was selected to map LULC change within the eThekwini Municipality.

3.3. Class area distribution

Figure 8 illustrates the changing landscape dynamics within the study area from 2002 to 2022. Over this period, notable trends include a consistent decline in cropland area from 31.5% in 2002 to 7.5% in 2017, followed by a slight increase to 11.8% in 2022, indicating potential shifts in agricultural practices or land use conversions. Forest cover remained relatively stable until 2017, experiencing a significant increase to 20.5% in 2022, possibly reflecting afforestation efforts or natural regeneration processes. Grassland areas fluctuated, peaking at 38.5% in 2012 before declining to 26.4% in 2022, suggesting changes in land management or ecological processes. Shrubland habitats remained consistent until 2017, with a slight expansion to 24.7% in 2022. Wetland and water areas showed minimal changes, indicating stable hydrological conditions. Impervious surface areas peaked in 2012 at 22.4% before declining to 13.8% in 2022, possibly reflecting urbanization patterns. Bare land areas remained relatively stable, suggesting minimal changes in this class.

Figure 8. Comparison of the area covered by each class after classifying all the Landsat images using XGBoost.

However, it was observed that when XGBoost classified Landsat 7 images, forest vegetation was underestimated compared to its classification of Landsat 8 and 9 images. Conversely, Landsat 8 and 9 classifications exhibited an underestimation of impervious surfaces compared to Landsat 7. Given the improved radiometric resolution in Landsat 9, the final 2022 classification was accepted and presented in Figure 9. Impervious surfaces are found in the eastern and northern parts of the eThekwini Municipality. Vegetated landscapes are found mainly in the western part of the municipality. It was also observed that forested areas are mostly fragmented. Fragmentation of forested areas was observed through the discontinuous distribution of forest patches, with numerous small and isolated patches scattered across the landscape.

Figure 9. The spatial arrangement of LULC classes within the eThekwini Municipality from 2002 to 2022.

4. Discussion

Improved access to medium and high-resolution remote sensing data coupled with free cloud-based geospatial analysis platforms such as GEE has enabled large-scale LULC classifications to be conducted more efficiently (Tassi and Vizzari 2020). Machine learning algorithms have also gained more interest in remote sensing studies because they have outperformed most conventional parametric techniques (Ngolo and Watanabe 2022). Therefore, this study utilized Landsat data from GEE to map LULC change using RF, SVM and XGBoost from 2002 to 2022 for the eThekwini Municipality. The performance of these algorithms was compared using mainly the overall accuracy extracted from the confusion matrix and quantity and allocation disagreements where XGBoost produced the highest overall accuracy.

This is consistent with conclusions in recent research findings such as Pan (2018), who compared the performance of XGBoost in predicting the concentration of 2.5 µm aerosol particles (PM2.5) in the air with artificial neural networks (ANNs), RF, multiple linear regression (MLR) and SVM. XGBoost was found to have outperformed all the other algorithms based on the root-mean-square error (RMSE) (17.298) and R² (0.9520). Yu et al. (2021) investigated the performance of XGBoost, RF and SVM in mapping the vertical forest structure, where it was observed that XGBoost (0.92) performed better than RF (0.90) and SVM (0.90) based on the F1 Score. In their study, Geng et al. (2021) applied RF, SVM, ANN, and XGBoost to predict seasonal maize biomass using field observation data and imagery from the Moderate-Resolution Imaging Spectroradiometer (MODIS). They found that XGBoost (0.77) and RF (0.77) performed better than ANN (0.72) and SVM (0.72) based on the R².

The advantage of XGBoost over other algorithms is based on its ability to reduce inaccuracies resulting from over- and underestimation and corrects the residual error, which results from the generation of a new tree based on the previous tree (Li et al. 2019, Geng et al. 2021). XGBoost can also accurately classify noisy data due to its loss function flexibility (Izadi et al. 2021). However, it should be noted that the performances of both RF and SVM were not significantly lower than that of XGBoost. Furthermore, when the performance of the three algorithms was investigated using the AUC, RF performed better than both XGBoost and SVM. This meant all algorithms used in this study could perform LULC classifications on a regional scale.

This study further evaluated the importance of variables used to classify LULC change. However, it should be noted that variable importance could not be extracted in RStudio for SVM, due to applying the RBF kernel. That is because, in a non-linear kernelised SVM, the separating plane exists in another space because of the kernel transformation of the original space. In RF classifications, the top four variables in terms of importance were the red band, blue band, green band and NDVI. The red and green bands provide critical information for detecting vegetation, given their sensitivity to a wider range of chlorophyll concentrations (Odebiri et al. 2022). Landsat’s green band is positioned above the green edge, a range of wavelengths where vegetation’s reflectance behavior changes, similar to the red band and red edge, both of which shift toward longer wavelengths in areas with high chlorophyll concentration due to chlorophyll’s absorption properties, serving as a crucial indicator of healthy and dense vegetation in remote sensing (Pastor-Guzman et al. 2015). Pastor-Guzman et al. (2015) added that as chlorophyll concentration increases, the red spectral band would reach minimum reflectance, whereas the green band continues to be sensitive. The blue band is useful in differentiating deciduous from coniferous vegetation and overall vegetation from bare soil (Li et al. 2019). NDVI is calculated using the red band and NIR and is a measure of vegetation greenness and health of vegetation in each pixel in a remotely sensed image (Kumar et al. 2022). This meant RF was more efficient at discriminating vegetation from other classes when performing LULC classifications.

The red band was also the most important variable in the XGBoost classification, followed by MNDWI, NDVI and the blue band. The MNDWI improves on the NDWI, which helps detect hydrological changes (Teng et al. 2021). It was successfully employed by Pal and Pani (2016) to detect changes in the Ganga River water surface coverage in India. MNDWI uses green and SWIR bands to detect open water features and to discriminate them from impervious surface features as these are often correlated in other indices (Zhang and Liu 2022). NDVI and the blue are critical for vegetation detection. This meant XGBoost was more efficient at discriminating vegetation and water bodies from other classes.

Given its better performance, the images classified using XGBoost were investigated for the area covered by each class within the study area boundary. The results showed that cropland, grassland, shrubland and impervious surfaces dominated the eThekwini Municipality landscape. Classifications also revealed that cropland, shrubland, and water exhibited a decreasing trend from 2002 to 2017, while grassland and impervious surface classes experienced an increase from 2002 to 2012. In 2022, cropland accounted for 11.8%, forest for 20.5%, grassland for 26.4%, shrubland for 24.7%, wetland for 0.7%, water for 1.3%, and impervious surface for 13.8% of the total land cover within the eThekwini Municipality. However, it should be noted that when XGBoost classified Landsat 7 images, it was observed that it underestimated forest vegetation compared to when it classified Landsat 8 and Landsat 9 images, which also underestimated the impervious surface compared to Landsat 7. This meant it was impossible to determine the extent of the increase in impervious surfaces after 2012. However, the improved radiometric resolution of Landsat 8 and Landsat 9 over Landsat 7 provided enhanced spatial and radiometric details (Mushore et al. 2022). Therefore, it was concluded that outputs from classified Landsat 8 and Landsat 9 images were more representative of the spatial arrangement of the eThekwini Municipality.

The 2022 LULC map output showed that impervious surfaces were most prevalent along the coast, which was along the eastern side of the eThekwini Municipality. Cropland was prevalent in the northern part of the municipality, where Tongaat Hulett, a sugar and maize company is, a major landowner (Moffett and Freund 2004). However, during the 1990s, Tongaat Hulett diversified its operations to include land management and property development, resulting in notable projects like gated townhouses, holiday complexes, shopping centers, and the establishment of the King Shaka International Airport in the region (Todes 2014). Hence, impervious areas have been expanding toward the north.

It was also mentioned in this study that the eThekwini Municipality is susceptible to adverse environmental consequences due to accelerated spatial change resulting from attempts to rectify inequalities and socioeconomic distortions created by apartheid through infrastructure projects (Sutherland et al. 2014, Musvoto et al. 2016). These government infrastructure projects further contributed to the increase in impervious surfaces and the decrease in vegetated landscapes. These areas are surrounded mainly by croplands, forests, and shrubland. Furthermore, with the projected increase in population numbers within the eThekwini Municipality and infrastructural expansion, it is probable that impervious areas will expand further and, therefore, will reduce the area covered by vegetation (Otunga et al. 2014).

The loss of vegetated landscapes means that natural ecosystem services that are vital for climate change mitigation, such as carbon storage, flood attenuation, and food production from these landscapes, are also lost (Tew 2019). The 2022 LULC map revealed that forested areas exhibit a high degree of fragmentation. Fragmentation in forested areas was observed by identifying discontinuities or breaks in the forest cover across the landscape. This was evident through the presence of isolated patches of forest separated by non-forested areas. Forest fragmentation is often associated with species loss, bird nest parasitism, and increased edge effects, which also result in increased predator density (Slattery and Fenner 2021). Increased predator density significantly impacts biodiversity through increased direct killing of prey and fear of predators, which changes some prey species’ behavioral patterns, physiology and reproduction rate (Yousef et al. 2021).

This study also underscores the significance of Landsat 9 in LULC classification, highlighting its utility and contribution to enhancing the accuracy of classification outcomes (You et al. 2022, Jombo and Adelabu 2023). By incorporating Landsat 9 imagery alongside Landsat 7 and Landsat 8, the study leveraged the improved radiometric resolution of Landsat 9 to refine the classification process, thereby providing more precise and detailed insights into LULC dynamics within the eThekwini Municipality. This aligns closely with findings by Shahfahad et al. (2023) and Ghasempour et al. (2023), who similarly compared the effectiveness of Landsat 8 and Landsat 9 for LULC mapping across a heterogeneous landscape. The successful integration of Landsat 9 underscores its potential as an asset in remote sensing applications, offering enhanced capabilities for monitoring and assessing environmental changes over time.

Therefore, by updating LULC maps for the eThekwini Municipality, this study aids decision-makers in mitigating the persisting challenges of climate change through appropriate environmental management policies and strategies. By applying applicable environmental policies, there will be fewer vegetated landscapes lost. This study also highlighted the applicability of remote sensing, machine learning and GEE data in LULC classifications. The utilization of GEE, in this study, is very significant given the finding by Chaves et al. (2020) that big data processing is still challenging for analysts. GEE eliminates this challenge by allowing analysts to access, analyze and visualize geospatial big data in powerful ways without needing for specialized coding expertise or supercomputers (Tamiminia et al. 2020). This study also employed Landsat 9, which was launched toward the end of 2021. Therefore, this study becomes one of the few studies that have used it for regional LULC classifications.

5. Conclusion

This study successfully classified LULC change for the eThekwini Municipality using machine learning algorithms and Landsat data. The findings illustrated that XGBoost classifications using Landsat imagery could perform better classification than RF and SVM. Therefore, the eThekwini Municipality’s LULC change was mapped using XGBoost and Landsat imagery. The produced maps will aid decision-makers in making socioeconomically sound plans that do not infringe on the environment. It is also recommended that the eThekwini Municipality update its LULC maps annually. Future studies can adopt the methodology used in this study to perform accurate LULC classifications. Also, studies that are not financially constrained can utilize high-resolution imagery to perform LULC classifications to improve accuracy.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Data availability statement

The data supporting this study’s findings are available on request from the corresponding author.

Additional information

Funding

This work was supported by the [National Research Foundation of South Africa (NRF) Research Chair in Land Use Planning and Management] under Grant [84157, 127354 and 114898].