Predictive Modelling of Land Cover Changes in the Greater Amanzule Peatlands Using Multi-Source Remote Sensing and Machine Learning Techniques

1. Introduction

In an era characterised by the escalating ramifications of global climate change, including severe weather events and species extinction, the management of peatland emerges as a pivotal concern [1]. Peatlands are unique ecosystems formed in waterlogged environments where slow decomposition of plant material leads to the accumulation of peat, a dense organic matter [2]. Such accumulations, often extending over millennia, establish them as significant carbon sinks, playing an essential role in climate change mitigation efforts. Globally, they are estimated to hold 650 billion tonnes of carbon on 3% of the Earth’s land surface, the equivalent to more than half of the carbon in the atmosphere or the carbon stored by Earth’s vegetation [1]. They also provide a wide range of ecosystem services including improvement of water quality and runoff regulation, and the provision of a landscape with cultural, recreational and livelihood values [3,4,5].

Despite their importance, peatlands face significant threats from anthropogenic land-use changes, including logging, plantation development, urban and industrial expansion, and accidental burning [6,7,8]. These contribute to the destruction of approximately 500,000 hectares of peatlands annually [6]. The United Nations Environment Programme [6] reports that the current rate of loss and degradation of healthy peatland is ten times faster than the rate of expansion over the last 10,000 years. It is noted that around 12% of the estimated 500,000,000 ha of global peatlands are severely degraded to the extent that peat formation has ceased, leading to carbon stock losses [6]. Such land use modifications decrease biodiversity and transform peatlands into net sources of greenhouse gases (GHGs), accelerating carbon dioxide emissions from drained peat soils [5]. It also triggers a loss of hydrological integrity and peat subsidence, amplifying the risk of flooding and altering nutrient storage and cycling dynamics [9].

In Ghana, the Greater Amanzule peatlands (GAP) reflect these global trends, experiencing ecological pressures due to agricultural plantation development, urbanisation and the burgeoning oil and gas industries [10,11,12,13]. Since the onset of oil and gas activities in 2011, these pressures have intensified, threatening the peatland’s ecological health and resilience [11]. This underscores the urgent need for accurate land cover information and predictive modelling to support effective management and conservation strategies for tropical peatlands like the GAP.

Previous research has extensively covered land cover mapping and change detection in tropical peatlands. Various remote sensing techniques have been employed to monitor and analyse land cover changes over time, contributing to a better understanding of these ecosystems’ dynamics [14,15,16,17]. Remote sensing data has been utilised to classify land cover types, detect changes, and assess the impacts of human activities on peatland environments [18,19,20,21,22]. Despite these advancements, there remains a notable gap in the predictive modelling of future land cover scenarios, which is crucial for proactive environmental management. This study therefore employs an integrated geospatial approach that combines multi-source remote sensing data with advanced machine learning techniques to enhance the accuracy of land cover classification and provide robust predictions of future land use scenarios in the GAP. Specifically, this study (1) utilises combined optical, Synthetic Aperture Radar (SAR), and digital elevation data for tropical peatland mapping and change detection in the GAP, and (2) predicts future land cover scenarios up to 2040.

2. Methodology

2.1. Study Area

GAP is in the Western Region (WR) of Ghana (Figure 1a). WR is bordered to the west by Cote D’Ivoire, to the east by the Central Region, to the north by Ashanti and Western North regions, and to the south by the Gulf of Guinea (Figure 1a). The peatland lies within the Coastal savanna and Rain Forest zones of Ghana (Figure 1a) and traverses the four coastal districts of Jomoro, Ellembele, Ahanta West and Nzema East (Figure 1b). The terrain predominantly consists of swampy lowlands and slightly elevated areas, with heights ranging from 0 to 196 metres above sea level (Figure 2).

The peatland’s hydrology is characterised by a network of small, interconnected watercourses and tributaries that flow into the three major river bodies in the area — Amanzule River, River Ankobra and the Tano River (Figure 2) — ultimately discharging into the Gulf of Guinea [23]. The hydrological processes are influenced by tidal fluctuations, rainfall, and groundwater dynamics, which collectively contribute to the formation and maintenance of the peat deposits [24].

The climate of the GAP is tropical wet, with an average annual precipitation of approximately 1500 mm and relative humidity of 87.5% [24]. The area has two rainy seasons: a major rainy season from April to July, and a minor rainy season from September to November. The dry season occurs from December to March, with the Harmattan winds bringing cooler, drier air from the northeast. Mean annual temperatures range from 24°C to 28°C, with slight variations due to the influences of the Gulf of Guinea and the various water bodies. These climatic variations may have implications for the hydrological processes and ecosystem dynamics, affecting the peatland’s structure and function [23].

Biodiversity assessment within the area has identified 400 species, of which 41 are presently recognised as threatened according to the International Union for Conservation of Nature (IUCN) Red List criteria [25,26]. The IUCN Red List evaluates species based on their risk of extinction, categorising them within a framework that ranges from ’Least Concern’, ’Vulnerable’, ’Endangered’, and ’Critically Endangered’, with the last category indicating the highest level of extinction risk for species that are still found in the wild [27]. In the GAP, four species, including the Slender-snouted Crocodile, White-thighed Colobus, Roloway Monkey, and the Hooded Vulture are identified as ’Critically Endangered’ [26]. Eight species are also classified as ’Endangered’, including the Chimpanzee, White-napped Mangabey, Diana Monkey, Western Red Colobus, White-bellied Pangolin, Krefft’s River Frog, Grey Parrot, and Cassava Croaker [25,26].

Records of biological groups such as butterflies, ants, spiders, fungi, reptiles and damselflies are currently non-existent for the GAP, suggesting that the published list of species is likely underestimated, and that the actual biodiversity in the area may be even richer than currently documented. The species counts mentioned here should therefore be regarded as absolute minimum estimates, emphasising the need for further research and surveys to assess the area’s biodiversity more accurately.

The socio-economic dynamics of the GAP is intrinsically linked to its ecological resources, with local livelihoods predominantly based on subsistence farming, fishing, and livestock farming. The region’s fertile soils and favourable climatic conditions support a variety of crops and cash crops, integral to both the local and national economy.

2.2. Land Cover and Classification Scheme

Classifying tropical landscapes from remote sensing data is complex due to structural complexity, high heterogeneity and the absence of a universal classification scheme for peatlands [28,29]. Defining appropriate thematic classes based on the characteristics of the study area and technical specifications of the imagery is therefore important. A classification system was adopted based on extensive field observations of the study area and a careful study of relevant literature [6,29,30,31]. Peatland classes were adapted from Lawson et al. [29], due to its detailed and comprehensive classification of peatland ecosystems, which mirrors the diversity found within the study area. The adapted categories comprised mangrove swamp, mixed swamp, palm swamp and bog plain; plantation comprised coconut, rubber and oil palm; artificial and bare classes included built-up land and bare surfaces, respectively. Other classes included sparse vegetation, natural forest, and water. In total, twelve classes were identified (Table 1).

2.3. Data

2.3.1. Satellite Imagery

A combination of optical, radar, and SRTM satellite data was used for the landcover classification (Table 2). This multi-source approach was adopted to enhance the accuracy of the land cover change analysis by leveraging the strengths of each data type [31]. All data were downloaded and processed using Google Earth Engine (GEE), a freely available machine learning and cloud computing platform capable of handling large datasets and having a comprehensive collection of preprocessed satellite imagery [32].

Optical data were obtained using Landsat-7 and Landsat-8 Collection-1 Tier-1 calibrated top-of-atmosphere (TOA) reflectance products [32]. TOA products are preferred due to their consistent calibration and pre-processing, which ensures uniformity across different scenes and sensors, and their simplicity for direct use without the need for additional atmospheric correction. Landsat data offers high temporal resolution and consistent coverage, enabling the detection of land cover changes over time [33,34].

ALOS/PALSAR and ALOS-2/PALSAR-2 satellites provided L-band Synthetic Aperture Radar (SAR) mosaic datasets [35]. SAR is important for its all-weather observation capability, especially relevant in the tropics where cloud cover can significantly hinder optical data acquisition. SAR’s ability to penetrate cloud cover ensures the collection of consistent and reliable data across all seasons [35,36]. The SAR data had a spatial resolution of 25 m and included both horizontal transmit-horizontal receive (HH) and horizontal transmit-vertical receive (HV) polarizations [35]. To ensure the accuracy of the backscattering coefficient layers, JAXA has pre-processed them using their Sigma-SAR processor, which involved radiometric and geometric calibration, ortho-rectification, slope correction, co-registration, and intensity tuning of neighbouring strip data before ingesting into GEE [37]. This pre-processing step is also important for minimising errors and ensuring data consistency [35].

Additionally, terrain features from SRTM at 30 m resolution were incorporated as ancillary data to account for the influence of topography on land cover distribution. SRTM data adds valuable contextual information to the analysis, enabling a more comprehensive understanding of the factors driving land cover changes [38].

Table 2. Satellite remote sensing data used for the landcover change analysis of the GAP

Data Type	Satellite /Source	Description	Collection Periods	Advantages	Source
Multispectral Optical	Landsat-7 and -8	TOA reflectance products	2010, 2015, 2020	Consistent coverage, enables detection of land cover changes over time [33,34]	U.S. Geological Survey
Synthetic Apperture Radar (SAR)	ALOS/ PALSAR and ALOS-2/PALSAR-2	L-band SAR mosaic datasets	2010, 2015, 2020	Penetrates clouds, all-weather observations [35,36]	Japan Aerospace Exploration Agency (JAXA)
Digital Elevation Model (DEM)	SRTM	Elevation data	2000 (used as ancillary data)	Provides contextual information on topography [38]	U.S. Geological Survey

Table 2. Satellite remote sensing data used for the landcover change analysis of the GAP

Data Type	Satellite /Source	Description	Collection Periods	Advantages	Source
Multispectral Optical	Landsat-7 and -8	TOA reflectance products	2010, 2015, 2020	Consistent coverage, enables detection of land cover changes over time [33,34]	U.S. Geological Survey
Synthetic Apperture Radar (SAR)	ALOS/ PALSAR and ALOS-2/PALSAR-2	L-band SAR mosaic datasets	2010, 2015, 2020	Penetrates clouds, all-weather observations [35,36]	Japan Aerospace Exploration Agency (JAXA)
Digital Elevation Model (DEM)	SRTM	Elevation data	2000 (used as ancillary data)	Provides contextual information on topography [38]	U.S. Geological Survey

The Landsat and L-Band SAR images were acquired for the years 2010, 2015, and 2020, corresponding to a decade of land cover observation. This timeframe provides adequate temporal range for observing changes in land cover patterns. The year 2010 also serves as a baseline before the onset of oil and gas activities in GAP, which began in 2011 [11]. The selection of 2015 allows for the observation of land cover changes 5 years into these activities. The inclusion of 2020 provides insight into land cover alterations a decade later, and importantly, after the emergence of small-scale gold mining activities, which became notable in 2016. This ensures a comprehensive analysis of the land cover changes in response to these significant anthropogenic interventions.

2.3.2. Reference Data

A combination of visual interpretation of very high-resolution Google Earth Pro imagery and field survey was employed to obtain the reference data. Prior information, including peatland inventory maps and pre-existing land cover classification [12,13,26,28] , guided the establishment of strata [28], ensuring representative sampling across the varied landscape. To ensure accuracy, stakeholders who have extensive experience and familiarity with the landscape were involved in the visual interpretation process as their expertise contributed to a better understanding of the ground conditions.This approach enhanced the representativeness of the reference data, accurately reflecting the different land cover types in the study area.

Reference data were represented as polygon features, which were then used to extract image pixels for training the classifier. This approach has been found to produce better classification outcomes than other methods, such as points (single pixels), point buffers (average pixel values), and image objects (area statistics) [39,40]. For the 2010, 2015, and 2020 image analysis, a total of 26,000, 25,700, and 30,000 image pixels were extracted as reference data, respectively. 30% of these reference data were randomly selected per class and set aside for testing. Image statistics were obtained based on the Landsat, L-Band SAR, and SRTM image stack using the corresponding region of interest. Reflectance plots were constructed to visualize the distributions of reflectance per land cover type against each optical band (See Figure A1 in Appendix). This approach ensured that the land cover classification scheme was grounded in both visual interpretation and quantitative analysis, leading to a more accurate and robust representation of the study area’s land cover characteristics [41,42].

2.4. Preparation of Image Features

Annual composite images for the years 2010, 2015, and 2020 were generated from the Landsat series using pixel-based compositing in Google Earth Engine (GEE). This method employed median statistics for composite generation, a robust technique for maintaining the representativeness of surface conditions by mitigating anomalies due to cloudiness and shadows [33,43,44]. Pixel-based compositing method acts as a countermeasure against challenges such as insufficient data availability, pervasive cloud cover, disruptions in image archives, atmospheric disturbances, and variations in radiometric readings induced by seasonal transitions or alterations in solar angles [45,46,47]. For the radar data, the digital numbers (DN) were first converted to power-scaled intensity values. Subsequently, a Refined Lee filter was applied to each channel (HH, HV) to reduce speckle noise. These values were then converted to sigma-naught (σ⁰) values, ensuring consistent and comparable backscatter values [35,48].

After these preprocessing steps, we calculated a series of indices highlighting greenness, moisture, soil properties, and structural variations of the landscape to increase the utility of the spectral and radar information contained in the original image bands (Table 3). These indices are derived from arithmetic combination of spectral reflectance measurements from Landsat imagery and backscatter coefficients from the SAR channels. Indices exhibit high correlation with vegetation characteristics such as phenology, biomass, moisture content, and leaf area index (LAI) [42,49], providing a comprehensive assessment of the ecological state of the landscape. To augment the analysis, texture features were extracted using the grey-level co-occurrence matrix (GLCM) from the SAR data. This texture analysis focused on the radar backscatter to discern the structural variations within the landscape, offering insights into the physical arrangement and condition of the vegetation and other surface features. Texture metrics included contrast (CON), dissimilarity (DIS), inverse difference moment (IDM), angular second moment (ASM), mean (SAVG), correlation (COR), entropy (ENT), and variance (VAR). These metrics capitalise on SAR’s unique ability to reflect the geometrical and structural properties of surface elements, enriching the land cover analysis with a level of detail unattainable through spectral data alone. The size of the neighbourhood included in each GLCM was set to 4, and the kernel size was $3\times 3$ square using the ‘glcmTexture’ function in GEE. The choice of kernel size followed literature in similar contexts [31,33,43].

Terrain features were generated from the SRTM data using the Terrain Analysis in Google Earth Engine (TAGEE) package [50]. The TAGEE package consists of a multi-band image containing the same data properties as the digital elevation model (resolution, data type, and coordinate reference system), with 13 bands (Table 3). Utilising all the terrain attributes of the package ensures a comprehensive analysis of the study area and provides the opportunity to test the relative importance of the features for tropical peatland (TP) delineation.

In total, 21 optical features, 28 radar features, and 13 terrain features served as input data in the machine learning model and analysis workflow (Table 3;Figure 3).

2.5. Feature Selection and Classification

The image features (Table 3) were standardized by subtracting the mean and scaling to unit variance prior to classification. To deal with the challenge of overfitting, recursive feature elimination (RFE) algorithm was used to reduce the number of features. RFE is a well-established iterative process that ranks and selects image features by eliminating the least important features based on stratified 2-fold cross-validation scores [51]. The algorithm creates a model with all the features, computes their importance scores, removes the least important features, and repeats the process until a desired subset of features is achieved. This process was employed via the Scikit-learn Python library [52], a widely used and robust library for machine learning applications. This was executed using Google Colab as the computational platform.

Identified important features were used to retrain a random forest (RF) algorithm for the classification. RF is an ensemble classifier that combines decision trees, bootstrap aggregation (bagging) and random subspace methods for classification and regression [53]. The combination of many weak learners in an ensemble contributes to RF achieving higher accuracy compared to machine learning algorithms based on a single classifier [53,54]. RF has become increasingly important in land cover classification in recent times because of its nonparametric nature, ability to limit overfitting and its flexibility [19,41]. It is also known for its high performance and efficiency in dealing with large input datasets with many different features [53]. Assessment of classification capabilities of 179 classifiers revealed RF as the best classifier among others that included support vector machine (SVM), decision trees, and neutral networks [55]. In their analysis of machine learning classifiers, Kaszta at al., [56] and Awuah et al., [39] also identified RF and SVM as the best performers among others that included k-nearest neighbours (kNN) and classification and regression trees (CART).

In the application of the RF classifier, two important parameters require careful selection to optimise the classifier’s accuracy and computational efficiency: the number of trees (ntree) and the number of features considered for splitting at each node (mtry). The ntree and mtry were set at 100 and 2, respectively. This configuration is grounded in findings by Nomura et al., [43], who demonstrated that classification accuracy does not significantly increase beyond 100 trees in a similar study, noting that a higher number of trees contributes diminishing returns to model performance. Also, the mtry value of 2 enhances the diversity of the model by limiting the number of features considered for each split, thereby reducing the model’s susceptibility to overfitting, and improving its robustness.

Feature importance scores were also computed to determine the relative contribution of each image feature in the classification of the different land cover types. Some of the most frequent approximations for feature selection include Gini Index [57], gain-ratio [58], and Chi-square test [59]. Feature importance scores for the overall classification were estimated using the RF-based Gini criterion.

RF uses Gini Index for feature selection at each node [31,39]. When assigning an input pixel to a class (Ci), for a given training set (T), the Gini Index measures feature impurity with respect to the different classes and is expressed as:

$$\begin{array}{c}\hfill \begin{array}{c}\sum \sum _{i\ne j}\left({\displaystyle \frac{f\left(C_i,T\right)}{\left(T\right|}}\right)\left({\displaystyle \frac{f\left(C_j,T\right)}{\left(T\right|}}\right),\hfill \end{array}\end{array}$$

(1)

where f (Ci, T) denotes the frequency of class Ciin the training set T, and ∣T∣ is the total number of elements in T. This formula iterates over all pairs of classes i and j, calculating the product of their relative frequencies, and then summing these products to get the overall Gini Index. A lower Gini Index value indicates a higher purity of the node, contributing to more accurate classification outcomes in the RF model [39,60]. Feature importance scores for each land cover class were subsequently estimated from the product of the overall feature importance estimates and the standardized mean value of each feature split for the given class.

A post-classification majority filter was applied to improve class homogeneity using a 3 x 3 pixel moving window [31,41,43]. This approach examines and potentially modifies the classification of each pixel based on the prevailing class among its immediate eight neighbours. By aligning the class of a central pixel with the majority class found within its surrounding context, this filtering technique effectively reduces the presence of isolated or anomalously classified pixels.

Accuracy of the classified images were evaluated using a confusion matrix [61], from which Precision, Recall, F-score and Overall Accuracy (OA) metrics were calculated usingEquations (2), (3), (4) and (5). Accuracy-adjusted estimates of tropical peatland area coverage were obtained following Olofsson et al., [62].

$$\begin{array}{c}\hfill \begin{array}{c}{\displaystyle Precision=\frac{tp}{tp+fp},}\hfill \end{array}\end{array}$$

(2)

$$\begin{array}{c}\hfill \begin{array}{c}{\displaystyle Recall=\frac{tp}{tp+fn},}\hfill \end{array}\end{array}$$

(3)

$$\begin{array}{c}\hfill \begin{array}{c}{\displaystyle F-score=\frac{Precision\times Recall}{Precision+Recall},}\hfill \end{array}\end{array}$$

(4)

$$\begin{array}{c}\hfill \begin{array}{c}{\displaystyle OA=\frac{tp+tn}{tp+fp+tn+fn},}\hfill \end{array}\end{array}$$

(5)

where tp, fp, tn and fn represent the number of true positive, false positive, true negative and false negative cases, respectively.

2.6. Predicting Future Land Use and Land Cover Changes with Cellular Automata and Artificial Neural Network (CA-ANN)

One aspect of this study was to predict land cover scenarios for the short term (2030) and long term (2040). This can be achieved through the application of LULC transition models, which help predict future land use scenarios by analysing historical data to understand the patterns and frequency of land cover changes [63,64,65]. Several models have been developed to forecast LULC and change detections, such as Conversion of Land Use and its Effects (CLUE), Markov chain, Clue-S model, Agent-based model (ABM), and grid-based LULC (GeoMoD) [66,67]. However, many of these models are limited by their theoretical underpinnings and methodological approaches, affecting their efficiency, sensitivity to policy changes, and ability to accurately simulate complex LULC transitions [68].

Cellular Automata (CA) on the other hand has been found to generate better spatial distribution of LULC categories than methods such as Clue-S, GeoMoD, and Markov chain [69], making it well-suited for short to long-term land cover simulations under various scenarios, such as socioeconomic and physical factors. CA offers dynamic simulation capabilities, high efficiency, simple calibration, and the ability to simulate multiple land covers and complex patterns spatially and temporally [65,70].

Despite its strengths, the CA model has some limitations, such as the inability to include driving forces in the simulation process. This limitation can however be minimised by integrating CA with other quantitative models, such as Artificial Neural Networks (ANN). The combined CA-ANN model is a robust approach for forecasting potential future LULC changes and has successfully been applied in subtropical and tropical areas to predict future LULC distributions [65,70,71]. The model covers both static and dynamic aspects of land cover changes and ensures higher prediction accuracy [72].

Given this robustness, the CA-ANN model was chosen for simulating future LULC changes using MOLUSCE (Modules of Land Use Change Evaluation), a plugin in the open-source QGIS Desktop. The application of this plugin involves several steps, including input data, area change analysis, modeling method, simulation, and validation.

2.6.1. Model Input Data and Selection of Explanatory Variables

The first step in applying the CA-ANN model is to input historical LULC maps and associated explanatory maps, from which patterns will be studied to make future predictions. In this study, the historical maps used were the classified maps for 2010 and 2015. Explanatory maps are spatial variables hypothesized to influence the changes in the territory. It was hypothesized that terrain features (i.e., DEM, slope, and aspect), population density, proximity of resources to roads, and proximity of resources to water influence land cover change. The choice of explanatory parameters was based on a review of similar studies. For instance, Kamaraj and Rangarajan [73] reported that proximity to roads and water influence significant variations in land use and land cover changes, often guiding the patterns of human settlement, agriculture, and resource exploitation. These variables serve as vital indicators that can delineate zones of varying vulnerability and propensity for change, thereby aiding in more accurate and predictive modelling of future land use trajectories. Kamaraj and Rangarajan [73] and Kafy et al. [74], also noted that terrain features such as elevation, slope, and aspect influence land use changes. The role of population density in shaping land use dynamics is also crucial. Within this context, distances related to roads and waterways were computed using vector representations of these features with the Euclidean distance function in ArcGIS software. The Euclidean distance measures the shortest path between two points in a plane or three-dimensional space using a straight line [75]. In two dimensions, the Euclidean distance between two points (x₁, y₁) and (x₂, y₂) is calculated as:

$$\begin{array}{c}\hfill \begin{array}{c}\sqrt{{(x_2-x_1)}^2+{(y_2-y_1)}^2},\hfill \end{array}\end{array}$$

(6)

2.6.2. Quantifying Magnitude of LULC Changes and Transition Potential (Modelling)

The model was run on QGIS to analyse the input data, identifying locations where changes have occurred and quantifying the magnitude of these changes. This step is crucial for understanding the spatial and temporal dynamics of LULC changes within the study area. Following this, the CA-ANN model was employed for LULC change transition potential modelling, which aimed to determine all possible changes that occurred on the landscape. CA spatial filter was then applied to the transition potential maps to simulate future LULC scenarios. The parameter requirements for the model, such as the maximum iteration (500) and neighbourhood pixel size (3 × 3), were set based on studies conducted in similar contexts [73,74]. This choice of parameters ensures a reliable and accurate representation of LULC dynamics, considering the spatial complexity of the study area. Using the LULC data from 2010 and 2015, explanatory variables, and the transition matrices, LULC projections were generated for 2020, 2030, and 2040.

2.6.3. Model Validation

The reliability, functionality, and acceptance of the model were evaluated using the kappa coefficient, inherent to the MOLUSCE plugin in QGIS. The kappa coefficient is an established statistical parameter for quantifying the level of agreement between predicted and observed LULC maps, beyond what would be expected by chance [64,76,77]. The coefficient is calculated with the following formula:

$$\begin{array}{c}\hfill \begin{array}{c}{\displaystyle Kappa\left(k\right)=\frac{P_o-P_e}{1-P_e},}\hfill \end{array}\end{array}$$

(7)

where $P_o$ represent the observed agreements, and $P_e$ is the expected agreement by chance. The observed agreement ($P_o$) is the sum of the probabilities of each category being correctly predicted:

$$\begin{array}{c}\hfill \begin{array}{c}{P}_o={\sum }_{i=1}^c{P}_{ii},\hfill \end{array}\end{array}$$

(8)

$$\begin{array}{c}\hfill \begin{array}{c}{P}_e={\sum }_{i=1}^c(P_{i+}\times P_{+i}),\hfill \end{array}\end{array}$$

(9)

P_iidenotes the probability of both raters agreeing on category i, P_i+is the sum of the probabilities across row i (the total predicted probability for category i), and P_+iis the sum of the probabilities down column i (the actual probability for category i). The subscript c denotes the number of categories.

The kappa value (k) ranges from -1 to 1, interpreting the agreement level: 0 indicates that the level of agreement is what would be expected purely by chance. A k value of 1 means there is complete agreement between the raters. Conversely, a negative k value suggests that the agreement is worse than random chance, indicating a systematic disagreement between the raters.

3. Results

3.1. Feature Selection for Enhanced Land Cover Classification

After executing the RFE, the optimal number of features was obtained, as shown in Figure 4. This represents subsets of feature variables that combine to produce the best classification accuracy. The results show that out of the 62 features evaluated for the landcover classification, only 27 were deemed important for the classification process. Consequently, the highest accuracy for the dataset was achieved when these 27 features were employed. Beyond the optimal point (27, 0.9852), the accuracy score fluctuated as the number of features increased, indicating the presence of irrelevant or redundant features that did not enhance classification accuracy and necessitated removal.

The features retained for the classification included six original image bands (radar and optical), twelve spectral indices, four texture features and five terrain features.Figure 5 shows the relative importance scores of the selected final input features in differentiating land cover categories. The top 10 most important features comprised five spectral indices, two original spectral bands, two terrain features, and one texture feature. Elevation emerged as the most importance feature, followed by the texture feature HH_savg (Figure 5). Among the original spectral bands, SWIR1 displayed the highest importance, followed by SWIR2 and NIR. The most relevant spectral indices were SAVI, GNDVI and MSAVI2, while the most influential terrain variables were Elevation, Slope and Hill shade. Texture features with the highest importance were HH_savg, HV_savg and NLI, respectively.

3.2. Classification Accuracy

The overall classification accuracy was high, achieving 93%, 94% and 94% for the 2010, 2015 and 2020 classifications, respectively (Table 4). High accuracy was consistently observed in the Natural Forest category, with F-score approaching 1, across all years, indicating robust classification performance. Mixed swamp and mangrove classes also showed improvement, particularly in F-score metrics over time. In contrast, categories such as sparse vegetation and coconut exhibited lower F-scores, highlighting areas where classification performance was less accurate and could benefit from further refinement. Consistent patterns of misclassifications were observed (See Tables A1.1-A1.3 in Appendix). Notably, sparse vegetation was often misclassified with natural forests and rubber plantations across all maps. Bare surface and built-up areas were frequently confused in the 2015 and 2020 maps, while mangroves were occasionally misclassified as mixed swamp. Error matrices for each classification are presented in Tables A1.1 - A1.3 in the Appendix.

3.3. Land Cover Changes

The study observed extensive land cover changes between 2010 and 2020 (Figure 6). A total of 114,028 ha, accounting for 21% of the entire GAP, experienced change from 2010–2015. A similar extent of change was observed from 2015-2020, affecting another 113,014 ha (21% of the area). When areas classified as water bodies are excluded from this analysis, the proportion of GAP that underwent land cover transitions—either through reduction or expansion of particular cover types — increases to approximately 44% for 2010-2015 and 43% for 2015-2020.

The spatial distribution of land cover and patterns of change across 2010, 2015, and 2020 are shown in Figure 7, while Table 5 presents land cover class areas and their respective changes throughout the study period. Although palm swamp forest experienced a substantial decline from 9647.91 ± 35.07 ha in 2010 to 6,325 ± 24.68 ha in 2020, peatlands in general increased from 54,502 ± 198.10 ha in 2010 to 61,072 ± 238.29 ha in 2020 (Table 5). Mangrove initially contracted by 297.05 ± 7.91 ha between 2010 and 2015, before expanding by 209.21 ± 7.93 ha between 2015 and 2020. Bog plain also expanded moderately, from 1,881.89 ± 6.84 ha in 2010 to 2,410.98 ± 9.41 ha in 2020.

Sparse vegetation experienced a significant reduction from 45,064 ± 163.79 ha in 2010 to 29,424 ± 114.81 ha in 2020. Rubber plantations on the other hand experienced considerable growth (85%), increasing from 30,530 ± 110.96 ha in 2010 to 56,617 ± 220.90 ha in 2020. Coconut and oil palm plantations diminished in extent across the three time periods (Table 5).

Artificial surface exhibited steady growth from 4,584 ± 16.66 ha 2010 to 4,946 ± 17.98 ha 2015 to 5252 ± 20.49 ha 2020. These findings underscore the dynamic nature of land cover changes in the GAP and highlight the importance of monitoring and understanding the driving forces behind these changes for better land use planning and sustainable development.

3.4. Relative Landcover Transitions

Examination of the landcover transitions show a landscape characterised by dynamic shifts (Figure 8). Establishment of rubber plantations are primarily occurring at the expense of forested areas, especially natural forest. Although there is evidence of peat swamp forest being converted to rubber (i.e., 4751.8 ha overall between 2010 and 2020), a more prominent trend of natural forest to rubber conversion is observed (i.e., 20,786.23 ha in the 2010-2015 period alone) (Tables A2.1 and A2.2 in Appendix).

The results also suggest that palm swamps act as pioneer habitats in the ecological succession towards mixed swamps by colonising the landscape before other tree species become established to form a mixed-species swamp (mixed-swamp). This succession is reflected in the observed decrease of palm swamp areas (3,149 ha from 2010-2015 and 2,920 ha from 2015-2020) and a corresponding increase in mixed swamp areas (Table A2.1 and A2.2 in Appendix).

An observed agricultural shift highlights a transition from areas characterised by sparse vegetation (including smallholder food crop farms) and coconut plantations, towards rubber cultivation (Figure 8). About 27,091 ha of sparse vegetation cover transitioned to rubber from 2010 to 2020, with an additional 1,482 ha moving from coconut to rubber within the same period. This shift likely reflects farmers’ and landowners’ responses to economic incentives, adjusting their crop choices to maximise financial returns.

Urban expansion was gradual, with built-up areas increasing by 360 ha (7.89%) from 4,561 ha in 2010 to 4,921 ha in 2015, followed by a further increase of 286.28 ha (5.82%) between 2015 and 2020, culminating in a total of 5,207.28 ha. A considerable amount of this urban development occurred on lands previously classified as sparse vegetation (Tables A2.1 and A2.2 in Appendix).

3.5. Projected Land Cover Changes

GAP is projected to remain largely stable, with minor modifications expected for 2030 and 2040 (Table 6). Modest expansion is forecasted for mangrove and mixed swamp areas, with each expected to increase by 0.73 ha by 2030. This positive trend may be attributed to conservation efforts by non-governmental organisations and their donor partners (Hen Mpoano, 2016a, 2017). Palm swamp and natural forest areas are anticipated to exhibit minor declines of 0.09 ha and 0.18 ha, respectively (Table 6), indicative of ongoing habitat loss and fragmentation due to agricultural expansion, logging, and infrastructure development.

Rubber plantations and built-up areas are projected to expand modestly, with a projected increase of 1.57 ha and 15.25 ha, respectively (Table 6). This expansion reflects the growing demand for agricultural commodities and the need for urban development to support increasing human populations [78,79]. Conversely, coconut plantations, bog plains, and sparse vegetation areas are projected to contract marginally, with a decrease of 4.28 ha, 4.45 ha, and 8.94 ha, respectively (Table 6). Projecting further to 2040, mangrove, mixed swamp, palm swamp, rubber, oil palm, and built-up areas are expected to display net increases of 4.90 ha, 34.42 ha, 6.60 ha, 2.79 ha, 13.55 ha, and 23.53 ha, respectively (Table 6). These trends suggest that pressures driving land cover changes, such as agricultural expansion and urbanisation, are expected to persist. On the other hand, natural forest, sparse vegetation, coconut, and bog plains are expected to manifest net decreases of 11.59 ha, 10.42 ha, 4.66 ha, and 13.47 ha, respectively (Table 6). This indicates that the loss and degradation of these ecosystems may continue, with potential adverse consequences for biodiversity, ecosystem services, and local livelihoods.

3.6. Evaluation of the CA-ANN Model for Land Cover Simulation

The overall Kappa value obtained for the CA-ANN model validation was 0.70, indicating a substantial agreement between the observed and simulated 2020 land cover maps. The accuracy rate was 80%, reflecting the percentage of land cover pixels that the model correctly classified when compared to the actual observed data during the validation process. The Kappa histogram (Category-Specific Kappa) value was 0.96, demonstrating an exceptionally high degree of consistency across different land cover categories. Kappa location value was 0.72, indicating a good level of spatial accuracy in matching the locations of observed and simulated land cover classes.

These results demonstrate the effectiveness of the CA-ANN model in accurately simulating land cover changes in the study area. It is however important to note that the accuracy of the CA-ANN model may be influenced by various factors, including the quality of the input data, the choice of model parameters, and the complexity of the underlying land cover dynamics [69]. Future research should therefore continue to refine and improve the model, as well as explore the use of alternative modelling approaches and data sources, to enhance the accuracy and robustness of land cover simulations [80].

4. Discussion of Results

4.1. Optimal Features for the Enhanced Land Cover Classification

The results showed that out of the 62 image features initially evaluated for the assessment, only 27 were considered essential for the classification process (Figure 4). These 27 features included six original image bands (radar and optical), twelve spectral indices, four texture features, and five terrain features, representing a diverse set of information that captures the unique spectral, spatial, and terrain characteristics of the land cover types in the study area.

Feature importance assessment identified elevation as the most important element for the land cover classification (Figure 5). This is consistent with our earlier findings [31], reinforcing the important role of terrain data in tropical peatland classification. Terrain characteristics offer insights into the physical and environmental factors affecting land cover distribution and are key in identifying the suitability of specific area’s activities such as agriculture, infrastructure development and conservation initiatives [50,79,81,82,83]. Radar-derived texture feature HH_savg was the second most important image feature in the classification, confirming the advantage of multi-source remote sensing data. This approach can be particularly effective for capturing the complex spatial heterogeneity of land cover types and offer a valuable solution in regions prone to cloud cover or where high-quality optical imagery is scarce [84,85].

Among the original Sentinel-2 spectral bands, SWIR1 band exhibited the highest importance, followed by SWIR2 and NIR bands (Figure 5). These findings align with research that have highlighted the enhanced sensitivity of SWIR and NIR bands in discriminating vegetation structures, moisture content, and land cover transitions [42,49,86,87]. Their wavelengths are adept at capturing distinct spectral signatures associated with different biophysical properties and land surface conditions [19].

The most relevant spectral indices included the SAVI, GNDVI and MSAVI2, which are known to provide robust measures of vegetation condition and stress, as well as minimise the influence of soil background and atmospheric effects [39,81,88]. The high importance scores of these indices confirms their effectiveness in differentiating between land cover types characterised by varying levels of vegetation cover and density, such as forests, plantations, and bogs [89,90,91].

Texture features with the highest importance were HH_savg, HV_savg, and NLI, respectively. According to Haralick et al. [92], texture features derived from the spatial arrangement and variability of pixel values within a given window, can effectively capture the structural and contextual information of land cover types and improve the classification accuracy, especially for classes with similar spectral characteristics but distinct spatial patterns. The inclusion of these texture features in the optimal subset of image variables demonstrates their utility in enhancing the discrimination power of the classifier and the robustness of the land cover classification process [93,94,95].

4.2. Classification Accuracy

The results revealed high overall classification accuracies of all land cover maps (Table 4). This can be attributed to the optimal selection of image features, the robustness of the random forest classifier, and the meticulous digitisation process that ensured the accuracy of reference data used for validation [40,61,62,96]. The satisfactory classification performance of individual land cover classes further confirms the effectiveness of the proposed methodology in capturing the spatial and spectral variability of different land cover types in the GAP. In particular, the peatland classes, including mangroves, mixed swamps, palm swamps, and bog plains, achieved F-scores exceeding 0.7 in most cases, except for mangroves, which displayed lower accuracies due to their misclassification with mixed swamps (Table 4). This may result from the spectral similarity between these two classes (See Figure A1 in Appendix), as well as their spatial proximity and the presence of mixed patches of mangroves and swamps in transitional zones [97,98]. However, the overall classification performance of peatland classes can still be considered satisfactory, considering the complexity and heterogeneity of these ecosystems and their susceptibility to various anthropogenic and natural disturbances that may alter their spectral characteristics [99,100].

Consistent patterns of misclassifications were identified among the other vegetation types (Tables A1.1 - A1.3 in Appendix). For example, the term ’sparse vegetation’ encompassed a range of vegetation within the study area, such as young plantation trees, rainfed croplands, and small-scale agriculture (Table 1).Consequently, misclassifications of sparse vegetation and other landcover types was expected. Also, newly established rubber plantations, characterised by land cleared of most trees and woody vegetation but possibly retaining some trees for shade or an herbaceous understory, could be easily mistaken for sparse vegetation.

Misclassifications observed in the 2020 land cover map between rubber plantations and natural forests (Table A1.3 in Appendix) may be attributed to the similarities in spectral reflectance and the spatial configuration of these classes such as the growth of mature rubber plantation crowns and improper planting spacing [101]. Bare surfaces and built-up areas may be misclassified also due to similar spectral appearance of dry soil and light impervious surfaces such as concrete [102,103].

4.3. Landcover Changes

GAP underwent significant land cover changes (Figure 6) — either through reduction or expansion of land cover types — from 2010 to 2020. Between 2010 and 2015, land cover alterations were observed across 114,028 ha, constituting 21% of the entire landscape. An equivalent area, amounting to another 21% of the total landscape, underwent changes from 2015 to 2020. Excluding water, the proportions of land cover change increased to 44% and 43% for the 2010-2015 and 2015-2020 periods, respectively. This recalibration serves to highlight the extent of terrestrial changes more clearly and offers valuable insights into the dynamics of human impacts on the landscape.

The spatial distribution of land cover and patterns of change across the three time periods revealed several notable trends and dynamics. For instance, the peatland classes, particularly mixed swamp and bog plains, exhibited expansion trends, whereas palm swamp forests showed a significant decline (Table 5). These changes may be related to various factors such as climate variability, hydrological alterations, natural processes, and anthropogenic pressures, including land-use changes [104,105,106]. Mangroves expanded from 2015 to 2020 and could be attributed to the successful implementation of conservation and mangrove restoration efforts by NGOs and local communities in the region since 2015 (Figure 9 ; [12,13]).

Bog plain expansions may be linked to natural processes such as peat formation and accumulation [81,107]. Conversely, reduction in sparse vegetation areas (encompass a diverse range of land cover types such as young plantation trees, rainfed croplands, and small-scale agriculture,Table 1) may be attributed to the corresponding expansion of rubber plantations, and urbanisation (Figure 7;Table 5). Rubber plantations experienced 85% growth, increasing from 30,530 ± 110.96 ha in 2010 to 56,617 ± 220.90 ha in 2020 (Table 5). This can be linked to the increasing global demand for natural rubber and favourable market conditions, which have incentivised land conversion to rubber plantations [108,109,110]. Declines in coconut and oil palm plantations may be related to factors such as disease outbreaks, changes in agricultural policies, or shifts in market preferences [111,112,113,114,115].

Built-up and bare areas exhibited steady growth from 2010 to 2020 (Figure 7;Table 5), reflecting the urbanisation and development in the GAP. This urban expansion may be driven by increasing population and reported infrastructural developments, particularly due to the expanding oil and gas activities in the region [12,23,24,116]. Conversion of lands for urban purposes can lead to habitat loss, fragmentation, and alterations in ecosystem functions, which in turn can have significant implications for biodiversity, ecosystem services, and human well-being [6,117].

The relative stability of water areas throughout the study period (Table 5) suggests that hydrological processes, such as surface runoff and evapotranspiration, have remained relatively constant over time [118]. However, further research is needed to investigate potential changes in water quality, groundwater resources, and aquatic ecosystems in response to land cover changes and human activities in the region.

These land cover changes highlight the complex interplay between natural processes, socioeconomic drivers, and policy interventions that shape land use and land cover patterns over time. Understanding these dynamics is important for informing land-use planning, resource management, and conservation efforts aimed at promoting sustainable development and enhancing ecosystem resilience in the face of global environmental change [105,119,120].

4.4. Relative Landcover Transitions

The assessment of land cover transitions within the GAP indicates significant shifts, with both losses and gains among various land cover types (Figure 8; Tables A2.1, A2.2 in Appendix). Establishment of rubber plantations was predominantly at the expense of natural forest and sparse vegetations (including smallholder food crop farms). This trend reflects broader patterns of deforestation for monoculture expansion in tropical regions [108,109,110,121,122], raising concerns about the long-term ecological sustainability of GAP and its capacity to support diverse ecosystems. Monoculture plantations can lead to soil nutrient depletion, increased pest pressures and greater vulnerability to diseases [108]. It could also have significant implications for local livelihoods and food security by reducing the availability of land for food production and potentially disrupt traditional land use systems [109].

Expansion of built-up and bare areas, collectively referred to as artificial surfaces, also occurred predominantly in areas previously occupied by sparse vegetation (Tables A2.1 and 2A.2 in Appendix). This confirms the ongoing urbanisation and infrastructure development within the GAP. Such changes can lead to reduced agricultural land availability, displacement of natural habitats, and potential disruption of local biodiversity and water cycles [109,123].

The results further suggest that palm swamps act as pioneer habitats in the ecological succession towards mixed swamps, initially colonising the landscape before other tree species become established to form a mixed-species swamp (Tables A2.1 and A2.2 in Appendix). This transition is likely influenced by factors such as hydrology and natural succession processes [29], with the establishment of mixed swamp forests potentially offering increased habitat diversity and ecosystem functionality [1]. Such transitions underscore the dynamic nature of peatland ecosystems in the GAP.

In general, these land cover transitions have important implications for land management and conservation efforts. The conversion of natural forests to rubber plantations highlights the need for policies that balance agricultural expansion with the conservation of natural ecosystems. This may include the promotion of sustainable land management practices, such as agroforestry, which can enhance biodiversity and maintain ecosystem services while also supporting rural livelihoods [81,124,125]. Protecting and restoring palm swamp and mixed swamp ecosystems can also provide valuable habitat for many species and contribute to the overall ecological integrity of the GAP. Conservation efforts should focus on preventing further degradation and fragmentation of these ecosystems, as well as promoting the recovery of degraded areas through reforestation or assisted natural regeneration [6,126].

Expansion of built-up and bare surface areas indicates that urban planning and infrastructure development should be integrated with biodiversity conservation objectives. This may involve implementing land-use zoning and green infrastructure strategies that minimise habitat fragmentation and enhance ecological connectivity [11]. It is also important to continuously monitor and analyse land cover changes to inform land management and conservation decision-making. Remote sensing and GIS can be a valuable tool in this process as demonstrated in this study, enabling stakeholders to better understand the drivers and consequences of these changes and to develop more effective strategies for sustainable land use and conservation [46,127].

4.5. Projected Land Cover Changes

Projections of land cover changes within the GAP for 2030 and 2040 were based on land cover transitions observed between 2010 and 2020, terrain characteristics, infrastructure proximity, and population density metrics. The results estimate minor changes in GAP (Table 6 ), implying that the ecological and developmental dynamics within the area are likely to maintain a steady state, with conservation efforts and development pressures achieving a delicate balance. However, these projections do not consider certain important factors capable of influencing the region’s land cover trajectory. Notably absent are considerations of climate change, the implications of proposed mining concessions, the potential impacts of a proposed $60 billion petroleum hub project, population growth and the absence of a formal management for the GAP [12,13,116,128].

The omission of these factors underscores the need for an updated projection model to integrate these variables for a more comprehensive and accurate forecast of land cover changes. Incorporating climate change scenarios for instance is essential, given their profound effects on ecosystems, agriculture, and water resources. Similarly, understanding the potential land cover impact of large-scale projects like the petroleum hub and mining concessions is crucial for effective land management and conservation planning. Also, the region’s higher-than-average national population growth rate [116] adds a layer of complexity, potentially accelerating land use changes and increasing pressure on peatlands.

This initial projection should therefore be viewed as an important baseline from which to gauge the potential impacts of the unconsidered factors, identify areas of both stability and vulnerability within the GAP and inform targeted conservation and development strategies. Integrating the overlooked aspects will allow for a better understanding of the dynamic interplay between natural processes and human activities and facilitate the development of nuanced strategies that ensure sustainability and resilience of GAP’s land cover. By doing so, policymakers, conservationists, and stakeholders can undertake informed actions that balance ecological preservation with socio-economic development, ensuring the long-term well-being of the region and its inhabitants.

4.6. Contrasting Projections with Existing Literature

This study projects a net increase in mangrove areas, a decrease in natural forest cover, and a modest expansion of rubber plantations and built-up area by 2040. Such projections contrast with global trends of mangrove degradation and diverge from literature suggesting widespread forest recovery [45,129,130,131,132,133,134]. The anticipated modest growth in rubber plantations and urban areas also differs from reports of their rapid expansion in certain regions [68,135].

This divergence highlights the unique ecological dynamics and land use pressures within GAP, emphasising the need for context-specific conservation strategies and land management policies. The specificity of the findings underscores the importance of localised environmental assessments in informing sustainable development practices, particularly in regions facing the dual challenges of preserving biodiversity while accommodating economic growth.

5. Conclusion

This study utilised multispectral and radar remote sensing data integrated with machine learning techniques to analyse and predict land cover changes in the GAP from 2010 to 2020, with projections extending to 2040. The findings indicate significant land cover changes, notably, an 85% increase in rubber plantations and a 6% decrease in natural forest cover, underscoring the impact of anthropogenic activities on the region. Future projections suggest minor changes in land cover by 2040. However, these projections do not account for potential impacts from climate change, large-scale development projects, and demographic shifts, which could significantly alter land cover dynamics. Continuous monitoring and adaptive management are therefore essential to ensure the resilience and sustainability of the GAP ecosystem. Overall, the research contributes to the broader understanding of land cover dynamics in tropical peatlands and emphasises the need for integrated conservation and development strategies. The methodological approach demonstrated here provides a robust framework for accurately monitoring and predicting land cover dynamics in tropical peatlands, which can be adapted for similar ecosystems globally. Policymakers and stakeholders are encouraged to utilise these insights to develop targeted actions that balance ecological preservation with socio-economic development, ensuring long-term ecological health and human well-being in the GAP region.