Accurate Distribution Assessment of Small Inland Wetlands in South-Kivu Province through Integration of Optical and SAR Data with Statistical Models


An overview of the study area

The South-Kivu province is situated in eastern DRC and shares borders with Rwanda, Burundi, and Tanzania. Covering an area of approximately 64,791 km2, it is one of the 26 provinces of the Democratic Republic of the Congo. South-Kivu accounts for around 2.73% of the total land area of the country42,43. The province is home to approximately 6.2 million people, with 47% of the population residing in rural areas44. Geographically, the province is located between latitudes 1.5836° and 5.0103° South and longitudes 26.8106° and 29.3890° East (Fig. 1).

Figure 1
figure 1

Territories of the South-Kivu province in Eastern DR Congo (map created using ArcGIS 10.7 Esri-TM: http://www.esri.com).

The province is divided into eight territories, namely (Shabunda, Kalehe, Uvira, Walungu, Kabare, Idjwi, Mwenga, and Fizi) including the city of Bukavu (Fig. 1), and each territory in municipalities referred to as “Groupement”. The province has a variety of landscapes ranging from low to very high altitudes. The elevation fluctuates from ~ 512 to 3464 m above the sea level (m.a.s.l) and decreases from east to west. The Congolese “Cuvette Centrale” begins in the Shabunda and Mwenga territories, whereas reliefs and valleys, including mountains and the Mitumba Chain Mountains, characterize the eastern and central territories. On the left, the Ruzizi plain is a broad plain extending into the Fizi territory, in the territory of Walungu and Uvira, as well as the highland favorable to the growth of a diverse type of wetlands and especially swamps, marshes, and peats44. Various wetlands occur due to these physical conditions (Figs. 2 and 4).

Figure 2
figure 2

Diversity of wetlands in South-Kivu, Eastern DR Congo (a) very high altitude lake Lubwe in Itombwe, (b) peatland in Lugana, (c) coastal of the Ruzizi river, and (d) flooded plain for rice production in Ruzizi plain. map created using ArcGIS 10.7 Esri-TM: http://www.esri.com.

South-Kivu province is characterized by a humid tropical and equatorial climate. According to Balasha et al.6, the mean annual rainfall of ~ 1500–1800 mm each year and an average annual temperature varies between 11° and 25 °C (Supplementary material S1) are observed in the province. Seven primary soil categories predominate in the South-Kivu mainly Haplic Acrisols, Humic Cambisols, Humic Ferrasol, Luvic Phaeozems, Mollic Fluvisols, Vertisols, and Gleyic Solonchaks45,46. Histisols (organic soil) are observed in some territories in small areas. The peat is used as energy fire material for food cooking in small quantities using small bricks. These bricks can replace a sizable quantity of charcoal or wood. In the South Kivu, peat is hand-harvested in the Kakonda peat (in Kabare), Hogola (at Nyangezi) in the Chiherano, and Kachandja peats (in Walungu). The hydrography is extensive and thick, with several small sources, large rivers (such as the Ruzizi, Elila, Ulindi, Itombwe, and Lwama), as well as lakes. Additionally, the province is a part of the “African Great Lakes” region. It has Lake Kivu, Tanganyika, and a few more small lakes and ponds (such as Lubwe in Itombwe). The vegetation comprises highland forests, herbaceous savanna, wooded swamps, and dense forest.

The land cover map was taken from the 2018 version of the European Space Agency’s (ESA) Climate Change Initiative Land Cover Project (CCI-LC) (cds.climate.copernicus.eu) product. The CCI-LC identified 21 land cover classes in South-Kivu (Supplementary material S2), which can initially be categorized using the FAO’s land cover classification system (LCCS). The province has protected areas, mainly the Kahuzi-Biega National Park (KBNP), the Natural Reserve of Itombwe (NRI), Maniema, South-Masisi, the Mont Kabobo (National Park of Ngamikka), and the Luama-Kivu Hunting Reserve47.

Materials and methods

A complete overview of the methodological flowchart can be found in Fig. 3. The method used is discussed in detail in the following subsections comprising data pre-processing, classification, training data, and statistical models used for prediction and validation. In this study, we adapted the methodology suggested by Mwita et al.38 for small inland wetlands of East Africa; and revised it according to the one proposed by Adeli et al.39, LaRocque et al.30, and Garba et al.32. We combined Sentinel-1 and -2 and ALOS-2 PALSAR data; four statistical classifiers were used. The first step was downloading satellite images. These images were obtained from the ESA and JAXA websites (https://www.eorc.jaxa.jp/ALOS/en/index_e.htm). Once obtained, the processing follows as the second step. It consisted of optical image extraction, correction, and band merge. The corrected bands were clipped to the study area (South-Kivu shapefile obtained from RGC: http://rgc.cd/), and indices were calculated. The indices calculation is presented in Supplementary material S3. The indices formulas were extracted from (https://custom-scripts.sentinel-hub.com/custom-scripts/sentinel-2/indexdb/). In total, 27 indices were used.

Figure 3
figure 3

Flowchart of the methods used to map wetlands in South-Kivu by using optical, SAR and field data and four machine-learning (ML) algorithms. Topographic, vegetation, hydrologic indices were used as variables subsequently calibrated and validated using field samples.

Distribution data

The term “wetland” used in this study is based on the definition suggested by Steinbach et al.48, Chuma et al.5, and Amler et al.7, and which is currently used in eastern Africa and adapted from the Ramsar definition49: ‘wetlands are areas of marsh, swamp, inundated valleys, peatland or water, whether natural or artificial, permanent or temporary, with static or flowing water, the depth of which at low tide does not exceed 6 m’. Conceptually, the “potential distribution, existing wetlands” method has served as the foundation for this study. It concentrated on prospective wetlands, representing wetlands’ size at the most significant level before human interference. In other words, regions where water-related ecosystems are most likely to develop, are highlighted through prospective wetland mapping. Archive points and polygon files constitute the presence data, which were collected by scientists, environmental local non-government organizations (NGOs), and mostly during our field works from September 2020 to August 2022. Other visual fieldworks of well-known sites were made using drone images and delineated in ArcGIS and Google Earth following the digitalization process. In total, 550 samples shapefiles were considered (Supplementary data S4). For inaccessible wetlands, the location’s geographic coordinates were taken from the edge and placed in the center once the image was obtained. More points were taken into consideration for the wetlands that had a large area.

Check new:   Shortlist for Dezeen Awards 2023: The 30 Most Innovative Sustainability Projects in the World

The total sample (shapefiles) was split into two datasets; 70% of plots (385 shapefiles) were used as training samples, and the remaining (30%: 165) as validation samples. Those points were considered “presence samples”. The “Absence” data comprised 2300 shapefiles extracted from the RGC comprising cities and main towns, schools, hospitals, villages, airports, farms, tree plantations, woodland, etc. These ground-truth surveys were conducted during the same period as the image acquisition.

Data on environmental variables

Since wetlands are characterized by specific topography, vegetation, and hydrology, topographic, vegetation, and hydrological indices were used during the mapping and delineation process. Topographic variables were derived from the digital elevation model (DEM) from ALOS-PALSAR (12.5 m resolution). This comprises the elevation (m), the slope (%), land aspect, curvature, Topographic Witness index (TWI), also known as the compound topographic index (CTI), estimated following the formula: TWI = (αtanb) where α is the local upslope area draining through a certain point per unit contour length and b: the local slope in radians. TWI is related to soil moisture influencing rapid runoff and flash floods50. It was calculated in ArcGIS 10.7 Esrti-TM. The ALOSPALSAR DEM was first projected, and the flow direction and flow accumulation were calculated. TWI was calculated as the ratio between the Flow accumulation and the slope. For Sentinel-1, the horizontal transmit and horizontal receive polarization (HH), horizontal transmit and vertical polarization (HV), and vertical transmit and vertical receive polarization (VV); the ratio γ =|HH|2|VH|2 and η =|VV|2|VH|2 were also calculated and integrated. For ALOS-2 PALSAR data, three features were integrated: HH, HV, and HH/HV (ρ =|HH|2|VV|2). They were corrected using Freeman–Durden (following surface scattering and double-bounce scattering), Cloude–Pottier (polarimetric decompositions and the compact polarimetric simulations) comprising alpha, beta, gamma, and lambda, multi-polarizations, dual polarizations, and polarimetric decompositions51. The SAR data preprocessing comprised speckle reduction, terrain correction, and geocoding following steps developed by Veci52, Foumelis et al.53, and Braun54. Since the province typically has two seasons, one dry season lasting from May to August and another dry season for the remaining months, two images were used, comprised of those taken in July (at the height of the dry season) and one in November (during the high rainy period).

This study considered hydrogeomorphic (HGM) wetland types: riverine, depression, slope, flat, and lacustrine fringe. Wetlands were mapped by modeling groundwater using environmental variables. The standard algorithms implemented steps in S1 Toolbox software include applying orbit files, removing low-intensity noise and invalid data on scene edges, removing thermal noise, radiometric calibration, and orthorectification51,54. The correction was made in the Sentinel Application Platform (SNAP) for optical data, starting with Rayleigh correction (computing the bottom of Rayleigh reflectance bands). The image resolution (in m) was resampled at 10 m. The sea level pressure and ozone were maintained at 1013.25 (in hPa) and 300 (in DU). Finally, radiance-to-reflectance conversion was used before indices calculation55. The dataset of exploratory variables included 27 environmental variables divided into 7 topographic and Hydrogeomorphology (HGM) variables, 8 SAR, and 12 vegetation indices. The description, formula, and source of these indices are presented in the Supplementary material S3. The radar backscattering was made following Sigma-nought (σo). It is referred to as the radar backscatter per unit area (m2/m2), expressed in decibels (dB). The standard formula used to calculate σo: σo = 10 * log10(DN2) + K, where DN is the image pixel digital number measured in the SAR amplitude image, and K is a calibration factor that varies depending on the SAR sensor and processor system used. As we use both ALOSPALSAR and Sentinel-1, the factor was fixed to − 83.0 dB and 0 dB, respectively.

Other ancillary information was used, including GPS points, on-site images obtained from drone missions (Fig. 4), field notes on dominant vegetation, accessibility from roads, villages, or markets, and land use. The GPS points were inserted into ArcMap and google earth, and then the boundary delineation was conducted using high very resolution images. The Universal Transverse Mercator (UTM), zone 35S, was used for projection, while the Geographic Coordinate System (latitude–longitude) WGS1984 was maintained for points. Three non-wetland classes: deep water, urban, and upland, were obtained and merged into one class called ‘no-wetland’.

Figure 4
figure 4

Model performances using AUC, Kappa and TSS, correlation coefficients. The four models were compared with the F-test of Welch and DeLong test. For each model 10 experiments were executed and helped for statistical comparison.

Model selection and construction

Four spatial statistic models were used comprising Artificial Neural Network “ANN”21, Boosted Regression Tree “BRT”56, Random Forest “RF”11,56 and Maximum Entropy ‘Maxent’35. It is a computing system inspired by the biological neural networks that make up animal brains, known as neural networks (NNs) or neural nets.

Check new:   IOC Aims to Support Paris 2024 and Future Hosts with Organising Olympic Dossier

(a) ANN is a widely used machine learning (ML) algorithm that can work for big data analysis. The multilayer feed-forward feature is the primary form of neural networks. ANN includes several neurons or nodes that function in parallel to convert the input data into output types. Typically, ANN consists of three-layer types, namely (i) the input, (ii) the hidden, and (iii) the output layers. Depending on the specific application in a network, each layer has some neurons. Each neuron is connected to other neurons in the next consecutive layer by direct links. These links have a weight that represents the strength of an outgoing signal57.

(b) The Random Forest (RF) model has been widely used to map this ecosystem58,59. According to Rapinel et al.11, due to its ability to consider a large number of variables from many sources and low sensitivity to outliers and over-learning, RF models have been more effective for mapping wetlands than other types of models like support vector machines (SVM), maximum likelihood, or decision trees (DT).

It is a non-parametric supervised ML algorithm inside RStudio and R 4.2.160. Because it can manage the significant difference variables and be used to neutralize noisy data, RF has demonstrated its use for classifications with the enormous volumes of data in satellite images. The input variables for the Random Forests were indices generated from the satellite images and sample points (presence and absence). The number of trees used in each Random Forest classification was set at 250, which has been found to be the ideal amount when accuracy and processing speed are considered. To reduce the tree depth, the minimum node size was set at five. The RF algorithm builds numerous bootstrapped, de-correlated random decision trees to categorize a dataset according to the mode. Indeed, when implementing the RF classifier must specify the number of decision trees and randomly chosen variables for dividing the trees. The RF technique was developed with 250 decision trees utilizing the stratified random sample points after multiple iterations and fine-tuning. An ensemble of 10 models was generated for each combination. The number of training samples was set to 5000 and maintained 10 RF models.

(c) Boosted Regression Tree (BRT): To maximize prediction accuracy while revealing information about pertinent variables and their interactions, BRTs adaptively construct several; basic regression-tree models and merge them into a multi-tree model. This differs from conventional regression methods, which result in a single prediction model. With the added benefits of boosting, which makes it possible to model nonlinear functions and improves robustness to data concerns like outliers, BRTs combine the flexibility of regression trees to accept all data types inside a model, including missing or non-independent data. BRTs are described in depth and detail by Berhane et al.59. A BRT model needs two crucial inputs: learning rate and tree complexity. While the latter identifies the number of nodes inside each tree and sets the number of variable interactions that are fitted, the former establishes each tree’s contribution to the final model60. As each tree contributes less to the overall model with a slower learning rate, more trees are in a BRT model. However, building the model requires more computing effort and observations the more trees there are. Averaging 550 trees per model, we used a learning rate of 0.005 and a tree complexity of 5. In order to reduce the inherent stochasticity in each model because of the subsampling and bagging that go into the construction of each tree, an ensemble of 10 models was created for each combination of inputs and then averaged. Boosted regression trees (BRT) combine the strengths of two algorithms: regression trees (models that relate a response to their predictors by recursive binary splits) and boosting (an adaptive method for combining many simple models to give improved predictive performance). The final BRT model can be understood as an additive regression model in which individual terms are simple trees fitted in a forward, stage-wise fashion.

(d) Maximum of Entropy (Maxent): Currently used for species distribution, Maxent assumes all presence locations of the same type (usually species), and there is little variability in the species’ niche preference for locations. It assumes that environmental conditions at extant wetland locations represent the fundamental niche for a particular type of wetland. Maxent cannot adjust for differences in sampling effort type of wetland. In our case, the remotely sensed wetland presence data ensured an unbiased sample, while modeling limitations were collectively addressed using extensive wetland presence coverages to perform validation studies. In this study, we have adapted the methodology developed by Rebelo et al.61. We set the number of iterations to 5000, allowing the model enough time to converge. All other Maxent settings were left at their default values. Ten models were also maintained. This technique uses environmental data (here variables) for several “background” sites and known presence locations. The raster files contain variables that were extracted. The result, which was displayed as a probability raster, was reclassified into two classes above-mentioned.

Check new:   Why the price of GTA 6 could be the highest in the entire industry?Latest

Data were first processed in ArcGIS 10.7 Esri-TM and R 4.2.1 to run the models. The R packages “raster”, “rgdal”, “virtual species”, RStoolbox”, “xlsx”, and “data.table” help to process the images and the training samples. These packages also helped with data stacking, data frame creation, and raster processing to the same extent, resolution, and projection system. All the data were brought to WGS84/UTM zone 35S. To run the models, the package “MASS”, “gbm”, for BRT, “dismo”, “SDMtune”, “NMeval”, “virtualspecies”, “maxnet” for Maxent, “neuralnet” and “caret” for ANN, “randomForest” and “caret” for RF. For the result presentation, “ggplot2” and “gpairs” helped design graphs. The output files are probability images (ranging from 0 to 1) that were split into two classes, “wetland” and “non-wetland”1,35.

Post classification

The post-classification stage of this study aimed to refine the wetland inventory by removing erroneously classified pixels. Indeed, before conversion, post-classification was made. Three spatial analyses were made; we started first with the “majority filter tool” followed by the “boundary clean tool”. The wetland shapefiles were smoothed using the boundary between areas by expanding and shrinking. The small clusters were processed using the “region group tool”62. The “polynomial approximation with exponential kernel smoothing algorithm” was finally used to smooth wetland shapefiles.

Accuracy assessment

To determine the level of accuracy of the small inland wetland map, it is necessary to validate the predicted wetted landscape using optical and SAR imageries. The second dataset (30%: 165) of wetlands was used for validation. These wetlands were physically surveyed to explore their nature and water seasonal availability to validate the predicted wetland areas. 300 sites were collected from Google Earth Pro images and selected from a field survey using GPS. They were used to assess the overall and local accuracy. The accuracy was assessed using the Area under the curve (AUC), Kappa coefficient, correlation, and True Skill Statistics (TSS). First, we started assessing the “Sensitivity” and “Specificity”. In fact, due to their complementarity and the fact that they comprise the core of the Receiver Operator Characteristic “ROC” curve, Specificity and sensitivity are two fundamental metrics in classification63,64. Sensitivity and Specificity in this study relate to the accuracy of the “wetland” and “non-wetland” classes producer, respectively. Both concepts have their roots in the literature on ecological presence-absence distribution models, such as species and habitat distribution models65, which are comparable to the probability of wetland-occurrence models developed. Here, sensitivity is a term used to define the percentage of presences (wetlands) identified; it measures the degree of omission errors. Specificity is the fraction of accurately anticipated absences (i.e., non-wetland) that occurred. The formula (i) and (ii) are used to ass the two parameters and combined in the confusion matrix as presented below. The two parameters are True Positive Rate (TPR) and False Positive Rate (FPR). Sensitivity is generally used to describe the proportion of observed presences (i.e., wetlands) that are correctly predicted as such and is a measure of the level of omission errors. Specificity, on the other hand, is a measure of errors of commission, representing the proportion of correctly predicted absences (i.e., non-wetland).

Sensitivity = \(\frac{\alpha }{\alpha + \delta }\) (i) and Specificity = \(\frac{\beta }{\beta + \gamma }\) (ii), with n = α + β + γ + δ. α (TP) is the true presence, β (TN): true negative, γ (FP): the false positive and δ (FN): the false negative.

AUC: is a comprehensive assessment of the model’s overall performance. It provides a conceivable categorization threshold and is interpreted differently; one is the probability that the model classifies a positive example more carefully than a negative one. AUC is a robust and commonly used measure of predictive model performance. The AUC ranges from 0 to 1. An entirely false prediction model will have an AUC of 0, whereas an accurate prediction model will have an AUC approaching1,63,64. The value 1 specifies that the diagnostic test is perfect, 0.5 stands for a worthless test, and 0 signifies the result is entirely wrong. According to Kanti et al.65, the AUC value of 0.90–1.00 indicates excellent, 0.80–0.90 means good, 0.70–0.80 means fair, and 0.60–0.70 means poor accuracy level.

Correlation coefficient (r2): The correlation coefficient (r2) is a standardized measure of the predictive accuracy of a model. The formula (iii) was used to assess the r2.

$${\mathbf{r}} = \frac{{\sum {\left( {{\text{Y}}_{{{\text{sim}}}} – \overline{{{\text{Y}}_{{{\text{sim}}}} }} } \right)\left( {{\text{Y}}_{{{\text{obs}}}} – \overline{{{\text{Y}}_{{{\text{obs}}}} }} } \right)} }}{{\sqrt {\sum {\left( {{\text{Y}}_{{{\text{sim}}}} – \overline{{{\text{Y}}_{{{\text{sim}}}} }} } \right)^{2} } } \sqrt {\sum {\left( {{\text{Y}}_{{{\text{obs}}}} – \overline{{{\text{Y}}_{{{\text{obs}}}} }} } \right)^{2} } } }}\quad\hfill \left( {{\text{iii}}} \right)$$

Yobs is the wetland, and Ysim is the value from the model. Ȳobs, Ȳsim are the average values of Yobs and Ysim. As for AUC, r2 varies from 0 to 1. 0.90–1.00 indicates excellent, 0.80–0.90 means good, 0.70–0.80 means fair, and 0.60–0.70 means poor accuracy level66.

Kappa statistic: The Kappa coefficient is more valuable than the overall accuracy as it indicates how the classification rate compares to the likelihood of correctly classifying pixels by chance. The formula (iv) was used to assess the Kappa coefficient.

$${\text{Kappa}} = \frac{{2 \times \left( {{\text{TP}} \times {\text{TN}} – {\text{FN}} \times {\text{FP}}} \right)}}{{\left( {{\text{TP}} + {\text{FP}}} \right) \times \left( {{\text{FP}} + {\text{TN}}} \right) + \left( {{\text{TP}} + {\text{FN}}} \right) \times \left( {{\text{FN}} + {\text{TN}}} \right)}}\hfill\quad {\text{(iv)}}$$

True Statistic Skill (TSS): The TSS, often called the Hanssen-Kuipers discriminant, compares the number of correctly classified samples to that of a hypothetically perfect classification while removing those correctly classified samples that could be attributed to a chance agreement. Equation (v) below is used to calculate the TSS, which provides a more objective measurement of classification accuracy comparable to the Kappa coefficient often employed in the literature on remote sensing.

$${\text{TSS}} = {\text{Sensitivity}} + {\text{Specificity}}{-}1\quad {\text{or}}\quad TSS = (\alpha \times \beta ) – (\beta \times \gamma )(\alpha + \gamma ) \times (\beta + \delta )\hfill\quad {\text{(v)}}$$

Source