Delineation of Cocoa Agroforests Using Multiseason Sentinel-1 SAR Images: A Low Grey Level Range Reduces Uncertainties in GLCM Texture-Based Mapping

Abstract

Delineating the cropping area of cocoa agroforests is a major challenge in quantifying the contribution of land use expansion to tropical deforestation. Discriminating cocoa agroforests from tropical transition forests using multispectral optical images is difficult due to the similarity of the spectral characteristics of their canopies. Moreover, the frequent cloud cover in the tropics greatly impedes optical sensors. This study evaluated the potential of multiseason Sentinel-1 C-band synthetic aperture radar (SAR) imagery to discriminate cocoa agroforests from transition forests in a heterogeneous landscape in central Cameroon. We used an ensemble classifier, Random Forest (RF), to average the SAR image texture features of a grey level co-occurrence matrix (GLCM) across seasons. We then compared the classification performance with results from RapidEye optical data. Moreover, we assessed the performance of GLCM texture feature extraction at four different grey levels of quantization: 32 bits, 8 bits, 6 bits, and 4 bits. The classification’s overall accuracy (OA) from texture-based maps outperformed that from an optical image. The highest OA (88.8%) was recorded at the 6 bits grey level. This quantization level, in comparison to the initial 32 bits in the SAR images, reduced the class prediction error by 2.9%. The texture-based classification achieved an acceptable accuracy and revealed that cocoa agroforests have considerably fragmented the remnant transition forest patches. The Shannon entropy (H) or uncertainty provided a reliable validation of the class predictions and enabled inferences about discriminating inherently heterogeneous vegetation categories.

1. Introduction

The mapping of cocoa commodity cropland is essential to the quantification of its ecosystem services and the disservices related to tropical forest cover loss. Agricultural land expansion, predominantly for oil palm, rubber, and cocoa plantations, contributes significantly to tropical deforestation [1,2,3]. Moreover, these commodity cropping lands provide different ecological services in terms of carbon sequestration, habitat provision, and the conservation of biodiversity [4,5]. Thus, a reliable and recurrent mapping of such cropping areas is crucial for customizing forest landscape management to land use expansion.

Agroforestry is an agricultural option for sustainable cocoa production. Cocoa agroforestry refers to a system in which a cocoa tree crop is grown in the understory of multistrata canopy trees [6], which provide various timbers, fruits, and non-timber forest products (NTFPs) [5,7,8]. Cocoa is a perennial crop of high economic importance. In tropical Sub-Saharan Africa [2,9], it provides about 70% of global dry cocoa bean export [10]. Regrettably, the expansion of cocoa production lands significantly contributes to the loss of forest cover [11,12]. Such expansions are somewhat specific to countries and production landscapes [13,14,15,16]; therefore, some are more destructive to forests than others. On a global scale, cocoa production was responsible for 57% of the global agricultural land expansion rate (132,000 ha year⁻¹) in the period 2000–2013 [2]. However, such figures need validation, at the national level, through mapping of actual cropping lands.

From an ecological standpoint, compared to intensive monostratum cocoa plantations and other high-canopy commodity crops, such as oil palm and rubber, cocoa agroforests sustain ecosystem services on a scale that is considered to be second to transition forests [17,18,19]. Regarding management in most cocoa-producing nations, the available literature rarely addresses the spatial mapping of cocoa production areas. Management projections for a production area are based on the Food and Agricultural Organization’s (FAO’s) database on crops, FAOSTAT. This database depends on sporadic annual country reports of a harvest area. Projections from these reports may not represent the actual ground reality [20]. In the case of Cameroon, cocoa is predominantly grown in small-scale agroforests of 1–3 ha [5,7,8]. Thus, the cocoa area of 123,120 ha in Cameroon’s center region production hotspot [21] is in effect the harvested area, which is based on seasonal records from local farmers and cocoa cooperatives. The National Cocoa and Coffee Board (NCCB) monitors and publishes the national statistics. Depending on the variety and propagation technique, established cocoa farms require approximately three or more years of crop tending before the harvest of first produce [22]. Consequently, the FAO records of the harvested area may be, at minimum, three years short of possible expansions in cocoa farms. On the assumption of continuous expansion of cocoa production land, the associated impact on forest cover is, therefore, far greater than management extrapolations made solely from published data on harvested areas. For the sustainable management of cocoa production landscapes, national government programs that stimulate the export of dry cocoa beans [2,9,21] need reliable and updated estimations of both harvest and expansion areas for cocoa agroforests.

The application of Earth observation data provides large-scale mapping of commodity cropping areas. To date, the discrimination of cocoa agroforest areas with a multistrata canopy, using optical reflectance and vegetation indices, has not proven to be successful [2,23]. Cocoa agroforests have a similar canopy structure to transition forests [16]. In moist tropical zones, a high frequency of clouds and atmospheric aerosols hampers the application of optical satellite data. Synthetic aperture radar (SAR) images, in contrast, provide cloud- and season-independent information about land surface features. Based on texture information extraction, the analysis of SAR images has been used for the discrimination of cropland [24,25] and forest biomass estimation [26]. Unlike optical imagery, which captures the reflectance of trees and forest canopies, SAR data capture the water content (a dielectric property) and structure (a geometric property) of target features. Information on the structure of target features is only provided if the target’s size is lower than or close to the wavelength of the SAR sensor. Thus, use of SAR imagery is determined by the sensor’s wavelength, and necessitates image-processing procedures that vary with the vegetation type and scale of assessment.

Long-wavelength SAR, such as L-band (λ ≈ 25 cm) SAR, provides details on volume scattering from branches and stems, which are essential for aboveground biomass estimation [27]. A SAR-based index, the Radar Vegetation Index (RVI), was developed for biomass monitoring using L-band data. However, the application of such an index requires the removal of the contribution from soil surface backscatter [28]. Although such bands may improve the discrimination of vegetation with a high tree canopy, they are less reliable for mapping features of low or sparse vegetation. When mapping features of sparse vegetation, other SAR wavelengths may be equally reliable.

The C-band sensor has a short wavelength (λ ≈ 5 cm), and the satellite systems operating in C-band, so far, have provided large temporal series of SAR images. Notable spaceborne C-band SAR image sources are the European Remote Sensing Satellite (ERS) 1 and 2 series, Radarsat constellation, and the Sentinel-1 constellation. Its utility for mapping tropical land cover is little-explored, especially in commodity croplands under the predominantly heterogeneous farming conditions in Sub-Saharan Africa. The C-Band SAR penetrates the vegetation canopy only to a limited extent. However, as in the case of the settlement and grassland land cover classes, SAR backscatter signals from forests have high temporal stability. Thiel et al. [29] showed that the contrast between these land cover classes and agricultural land is high in cross-polarized (horizontal transmit, vertical receive (HV) or vertical transmit, horizontal receive (VH)) SAR image bands. Stimulus et al. [30] reported that texture measurements are needed to discriminate settlement areas from forests. Thus, considering the seasonal changes in the structure and water content of vegetation canopy elements, a temporal metric of texture from C-band SAR images may be able to discriminate perennial agroforestry land cover.

Texture measurements from the grey level co-occurrence matrix (GLCM) provide reliable information on the spatial relationship of an image’s pixels [31]. The GLCM provides a joint probability distribution or co-occurrence frequency of the grey levels (or intensity tones) in an image based on three parameters: pixel(s) distance, angular displacement, and image sub-region analysis window size. Several second-order, i.e., between two pixels, statistics from the GLCM are proposed in [32] to describe the texture in an image. The use of GLCM texture measurements depends on the geometry of target features and their characteristic spatial structure in the landscape [31]. For land cover classification in a heterogeneous landscape, Mishra et al. [33] observed that texture information was more valuable for improving the classification accuracy in an SAR image than in an optical image. However, the authors of [33] suggested that an optimum combination of texture features is needed for the specific type of landscape heterogeneity under consideration. Land cover classification using GLCM texture extraction has focused on the scale or window size [31,33,34,35,36]. However, the importance of grey level quantization in a GLCM texture analysis has been emphasized [37,38,39]. Moreover, for such texture feature extraction, the use of grey levels beyond the depth of pixels (the range of values) may increase the uncertainty in the results [40]. Thus, although applications of the GLCM matrix for land use and land cover (LULC) classification [30,32,37] do not take grey level quantization into account, this aspect may be vital in mapping heterogeneous agricultural landscapes; in particular, the inherently heterogeneous LULC categories.

The purpose of this study was to assess the temporal contribution of SAR volume scattering, essentially by vegetation canopy, to the discrimination of perennial cocoa agroforest land use. We used multi-seasonal and multipolarization Sentinel-1 C-band SAR images: (1) to evaluate the performance of GLCM texture-based discrimination of cocoa agroforest land use from transition forest cover in comparison to a typical classification from a multispectral optical image using a RapidEye image; (2) to assess the contribution of grey level quantization to improving the texture-based classification’s performance by comparing four different grey levels of quantization or dynamic pixel ranges (32 bits, 8 bits, 6 bits, and 4 bits); and (3) to assess the information gained from Shannon Entropy (H) or uncertainty as a classifier performance estimator.

We applied the GLCM to estimate four texture statistics (contrast, entropy, variance, and correlation) based on [31]. The texture statistics provided texture information on the structure and water content of vegetation (volume scattering). Then, we averaged the SAR volume scattering values across seasons by using a machine learning classification algorithm (Random Forest). We included other land cover classes in the classification analysis to derive a thematic land cover map of the heterogeneous landscape.

This paper extends the analysis and results in a conference paper [41]. We provide literature on texture analysis of SAR images and specify the research gap in the context of mapping agroforests land use. The data and image processing steps are described meticulously in Section 3—with added illustrations of radar backscatter intensity and the flowchart of image processing and analysis. We conducted an additional experiment to assess the contribution of grey level quantization on classification accuracy and uncertainties in the predicted land cover/use. We extensively illustrate the major results, and thereafter provide a discussion in the context of cocoa agroforests land use mapping using C-band SAR from Sentinel-1. Finally, in Section 5 and Section 6, we suggest application and management implication of the methods and results.

2. Study Site

This study was conducted in the landscape of Bakoa (in the UTM zone 32N, and bounded by the UTM coordinates 734280 m E 510975 m N and 747435 m N 501480 m E—a surface area of 123.28 km²), which is located in the Bokito District of the Mbam and Inoubou Department in the center region of Cameroon (Figure 1). This area is classified as a savannah–forest transition zone. The topography features a rolling terrain, and the altitude ranges between 500 and 900 m above sea level (a.s.l). The vegetation is a mosaic of bush-savannah, subsistence farming, and perennial cocoa agroforests. These perennial agroforests were established mainly within or along patches of transition and gallery forests. The study area is situated in a bimodal humid forest agro-ecological zone, which is characterized by two dry and wet seasons. The total annual rainfall ranges between 1300 and 1500 mm, and there is a long rainy season from August to November. The main dry season lasts for approximately 5 months (from November to April). The mean annual temperature is 25 °C.

3. Materials and Methods

3.1. Satellite Data: Optical and Radar Imagery

We acquired a multispectral optical image of 5 m spatial resolution from RapidEye, which was recorded in the dry season of 2015. The image comprises five spectral bands in the Blue (400–510 nm), Green (520–590 nm), Red (630–685 nm), RedEdge (690–730 nm), and Near Infrared (760–850 nm) range of the electromagnetic spectrum. Four different image tiles, acquired on the same date, were needed to cover the study landscape.

We accessed Sentinel-1A C-band (λ = 5.5 cm) SAR images for the study area from the Sentinel Scientific Data Hub of the European Space Agency (ESA). The SAR data were acquired with dual (VV (vertical transmitted, vertical received) and VH (vertical transmitted, horizontal received)) polarization in the Interferometric Wide (IW) swath imaging mode. We used Level-1 preprocessed Ground Range Detection (GRD) images of 10 m spatial resolution. We selected a temporal series of 50 images, acquired between March 2015 and April 2017, which covered both the dry and wet seasons. Using the image-processing tools of the Sentinel Application Platform (SNAP) version 5.0, we prepared image subsets and preprocessed them sequentially from radar backscatter intensity values to sigma naught (sigma0). The SAR image preprocessing chain, in SNAP toolbox, was as follows: thermal noise removal, apply orbit files, radiometric calibration, and geocoding. The digital elevation model (DEM) of the Shuttle radar topographic mission (SRTM) was applied in SNAP with the Sentinel-1 toolbox (S-1TBX) for terrain correction and geometric rectification of the SAR images. We used both the co- (VV) and cross-polarized (VH) bands of all images. We then projected the preprocessed 10-m resolution images in World Geodetic System (WGS) 1984 Universal Transverse Mercator (UTM) Zone 32 N. The main remote-sensing data that were used for analysis are presented in

Table 1. The single-date RapidEye image (5 m resolution) and 10 subsample (of the 50 selected initially) multi-date and multi-seasonal Sentinel-1 synthetic aperture radar (SAR) images (10 m resolution).

Satellite Mission	Scene ID(s)	Acquisition Date (DD/MM/YYYY/)	Sensing Stop Time (HH:MM:SS UTC)	Acquisition Mode (Polarization)	Data Level
RapidEye: RE-3	3241224_ 3241225_ 3241124_ 3241124_	9 January 2015	10:35:41.00	MSI, Optical	L3A
Sentinel-1A	_006256_008304_78DE	6 June 2015	17:28:11.147769	IW Ascending (Dual: VV,VH)	Level1 GRD
	_007306_00A05D_2111	17 August 2015	17:28:14.323283
	_007831_00AE86_4926	22 September 2015	17:28:15.577539
	_008706_00C641_B612	21 November 2015	17:28:15.454239
	_010456_00F838_64CF	20 March 2016	17:28:13.302867
	_011156_010D64_7E35	7 May 2016	17:28:15.219784
	_012031_012962_8F08	6 July 2016	17:28:18.678294
	_012906_01465C_878C	4 September 2016	17:28:21.557129
	_014831_0182BC_16C4	14 January 2017	17:28:19.166952
	_015706_019D94_BC50	15 March 2017	17:28:18.681981

MSI, Multispectral Imager; IW, Interferometric Wide Swath; GRD, Ground Range Detected.

.

3.2. Field Campaigns

During the field campaigns, which were conducted in 2015, 2016, and 2017, we collected ground information on land use and land cover. The field data comprised ground Global Positioning System (GPS) information and an inventory of representative areas that characterize the different land cover types and land uses in the landscape (see Figure 2 and

Table 2. A description of the thematic land cover types used for the classification of land cover (Figure 2).

Class Acronym	Class Name	Description
Bu	Built up	Residential, commercial/market, industrial, and administrative settings
Es	Earth road/bare soil	Land areas of exposed soil and bare rocks
Sv	Shrub/grassland Savannah	Imperata sp. savannah land: shrubby and grassland areas that have not been converted to farmland
W	Water	Rivers, ponds, and seasonal and permanent swamps
Af	Perennial cocoa agroforests	Land areas used for cocoa production with various degrees of canopy stratification. The canopy/shade trees are mainly deciduous
Fa	Subsistence farming	Savannah and forest land areas that have been converted essentially for permanent or seasonal subsistence crop production, including farm fallows
Sf	Transition/Secondary forests	Disturbed and gallery forest patches, secret/cultural forests, and hunting forests. These forests have a rather permanent and less stratified canopy structure

).

3.3. Image-Processing Workflow

The RapidEye images were preprocessed using the following protocol: atmospheric rectification by dark object subtraction (DOS), radiometric calibration to reflectance values, geometric correction, and the computation of different vegetation indices from a mosaicked image.

After subsetting the SAR images, we used the batch-processing mode of SNAP for the following preprocessing steps: radiometric calibration to Sigma0 (decibels) and geocoding with the SRTM 3-s DEM using RangeDoppler Terrain Correction. We used an intensity backscatter profile (Figure 3) and the Random Forest (RF) important variables criterion, the mean decrease in Gini coefficient [42], to select a subsample of 10 (of the 50) important images that represent six wet and four dry seasons between 2015 and 2017 (see Table 1). The seasonal variation in backscatter intensity, for 4 (of the 10 subsample) images, is shown in Figure 4.

The image-processing steps, which are detailed in the following subsections, were comprised of three major categories: feature extraction, image classification, and post-processing. Feature extraction consisted of computing images of vegetation indices and GLCM texture images. During image classification, we coregistered the vegetation index and texture images into the separate stacks or models described in

Table 3. The image stacks used to evaluate the Random Forest (RF) algorithm’s classification accuracy.

Data Categories	Model	Image Stack
Dry season Multispectral RapidEye Image (RE).	RE1	TOA Reflectance of B, G, R, Red Edge, and NIR: 5 Bands
	RE2	TOA Reflectance and Vegetation Indices (VIs): 10 Bands
Multidate and season SAR GLCM Textures (GL).	GL1	Multi-date VV GLCM Textures: 40 Bands
	GL2	Multi-date VH GLCM Textures: 40 Bands
	GL3	Multi-date VV and VH GLCM Textures: 80 Bands
Multidate and season SAR intensity and GLCM Textures (GLI).	GLI1	Multi-date SAR VV Simga0 intensity and VV GLCM Textures: 50 bands
	GLI2	Multi-date SAR VH Sigma0 intensity and VH GLCM Textures: 50 bands
	GLI3	Multi-date SAR VV plus VH Sigma0 intensity and, VV plus VH GLCM Textures: 100 bands

NIR, near-infrared; VIs, vegetation indices; GLCM, grey level co-occurrence matrix; TOA, top of atmosphere.

. Then, we ran eight RF ensemble (machine learning) classification algorithms using the image stacks as input (Table 3). Post-processing consisted of the estimation of uncertainties in the classified maps, in addition to accuracy metrics, as the basis for validating the classifier models. Finally, for the texture-based model with the highest overall accuracy, we evaluated the GLCM texture images at four different grey levels of quantization to improve the classification uncertainties.

3.3.1. Feature Extraction: Vegetation Indices (VIs) and GLCM Texture Features

The monitoring of vegetation status and extent is often based on the normalization ratios of spectral bands in the Visible and Near-Infrared (NIR) spectrum [43] in spaceborne imagery. These ratios are based on contrasting spectral responses of vegetation to the Red and NIR wavelengths.

The application of indices, such as the Normalized Difference Vegetation Index (NDVI), in vegetation monitoring has faced several challenges [44]. One notable issue is biomass saturation above certain thresholds, which is common in moist tropical vegetation. Although saturation may not be an issue over agricultural landscapes, reflectance from the soil background often perturbs the discrimination of sparse vegetation or cropland from bare soil [45]. In this study, we used VIs whose values indicate the status and abundance of vegetation and biomass and that minimize the effect of soil background on vegetation reflectance values [46]. We used the NDVI, the green NDVI (gNVDI), the Enhanced Vegetation Index (EVI2), the Soil Adjusted Vegetation Index (SAVI), and the Modified SAVI (MSAVI) [47,48,49,50,51]. However, to provide additional information on vegetation characteristics and vitality, recent optical sensors have included an additional spectral band: the red-edge band featured in RapidEye and Sentinel-2 [52]. This band is located between the red absorption (by chlorophylls) zone and the NIR waveband. Since radar backscatter signals from a ground resolution cell are pseudorandom, the interaction of microwaves with terrain objects may be difficult to predict. Moreover, SAR images have a speckle effect because the response signal of a resolution cell is a form of coherent interference from multiple scattering elements within the cell. Based on texture information extraction, an analysis of SAR images has been used to discriminate cropland [25] and estimate forest biomass [26]. Often, the GLCM statistical approach is used to estimate SAR textures. The GLCM is a sparse matrix that stores the co-occurrence probabilities of interpixel grey levels in an image [32]. These probabilities provide a second-order measure for texture features in an image. They represent the conditional joint probabilities of all pairwise combinations of grey levels (G) in the spatial window of analysis and depend on both the spatial orientation (θ) and displacement distance (δ). The computation of a GLCM is faster for images with fewer grey levels because the matrix is dimensioned to G. The conditional probabilities are estimated as follows:

$$\mathrm{Pr}\left(x\right)=\left\{C_{ij}|\left(\theta ,\delta \right)\right\}$$

(1)

where C_ij is the co-occurrence probability between grey levels i and j, and is defined by:

$$C_{ij}=P_{ij}/{{\displaystyle \sum }}_{i,j=1}^G{P}_{ij}$$

(2)

where P_ij is the number of occurrences of grey levels i and j within the given window for a pair (θ, δ); and G is the quantized number of grey levels. The denominator sums up to the total number of grey level pairs (i, j) within the analysis window.

Although different second-order statistics are commonly used to classify single images [53], some GLCM texture measures are autocorrelated [32]: a few selected texture measures may be sufficient to achieve specific image analysis objectives [31]. We assessed the accuracy of SAR images using four less-correlated GLCM texture measures: Contrast, Entropy, Correlation, and Variance [31]. We estimated the GLCM texture measures using a 5 × 5 moving window, an aggregate orientation of four directions (0°, 45°, 90°, and 135°), and a one-pixel displacement (interpixel distance).

$$Contrast={\displaystyle \sum }_{i=0}^{G-1}{\displaystyle \sum }_{j=0}^{G-1}{P}_{i,j} {\left(i-j\right)}^2$$

(3)

$$Entropy= {\displaystyle \sum }_{i=0}^{G-1}{\displaystyle \sum }_{j=0}^{G-1}{P}_{i,j}{\log }_2{P}_{i,j}$$

(4)

$$Variance={\sigma }_i{}^2={\displaystyle \sum }_{i,j=0}^{G-1}{P}_{i,j}{\left(i- {\mu }_i\right)}^2$$

(5)

$$Correlation={\displaystyle \sum }_{i,j=0}^{G-1}{P}_{i,j}\left[\frac{\left(i- {\mu }_i\right)\left(j- {\mu }_j\right)}{{\sigma }_i{\sigma }_j}\right]$$

(6)

where P_i,j is the joint probability distribution of the grey levels i and j at two ends of a displacement vector in the assessment window, and G is the number of rows or columns. Since we considered a symmetrical GLCM, ${\mu }_i\equiv {\mu }_{j }$ and ${\sigma }_i{}^2\equiv {\sigma }_j{}^2$. For Entropy, 0 × ln(0) = 0, since ln(0) is undefined.

3.3.2. Classification: The Random Forest Ensemble Algorithm

Mature Machine-Learning (ML) algorithms, such as Random Forest, Support Vector Machines (SVMs), single and boosted Decision Trees (DTs), and Artificial Neural Networks (ANNs), are frequently used for land use/land cover classification. These algorithms are able to model complex class signatures, can accept a variety of input predictor data, and do not require assumptions about the data’s distribution (i.e., they are nonparametric) [54]. Several studies have shown that these methods outperform traditional parametric classifiers, especially for complex data with a high-dimensional feature space. The literature on remote sensing, however, does not provide straightforward advice about the choice of algorithm. However, RF and SVM yield high and comparable classification accuracies [55]. The choice between the two may be case-specific and depend on the characteristics of the classes being mapped, the training data’s quality, the predictor variables, and the optimization parameters [54].

In this study, we decided not to work with a single DT or an ANN because these methods have been shown to be sensitive to the training data’s size and quality [56]. We preferred RF over SVM due to its relative robustness with respect to training set size and parameter settings and, most importantly, because RF works quite fast with complex, high-dimensional datasets, such as the multiseason co- and cross-polarization intensity and texture image bands used in this study. Making SVM fast on large datasets remains the subject of research. In addition, RF is able to model nonlinear concepts, such as backscatter and texture features, and account for their variability, which the SVM subsets of predictor variables (the support vectors) that define the boundary or margin conditions may exclude. Moreover, RF allows us to assess classification uncertainties at the pixel level [57], which was one of the objectives of our study. Finally, the choice of RF was supported by several recent and successful classification experiments [25,53,58,59,60].

The RF ensemble classifier algorithm builds multiple decision trees for the same dataset based on random bootstrapping of sample training data [61]. The classifier is less-influenced by the common issue of overfitting and is able to handle a large number of variables [42]. Firstly, each tree is built from a random subset (n) of two thirds of the original samples (N) (the “in-bag” data). Secondly, from a subset (m) randomly selected from the total (M) number of variables in the dataset (mtry), in each decision tree nodes are split using the variable that yields the highest decrease in impurity (the “best split” variable) [42]. The algorithm is a soft classifier based on the probability of pixels belonging to the considered classes (Table 2). Compared to other nonparametric classification algorithms, it is less constrained by the need for extensive training and test data samples. This is due to an integrated out-of-bag (OOB) error estimation and accuracy test that follows bootstrap subsampling on the input data. Additional details on the random forest algorithm can be found in [42,62].

We ran eight RF models for the different image stacks as classifier inputs (Table 3). For each model, we evaluated the OOB error curve and mtry to prune the decision trees down to an optimal number. To yield a spatially explicit and unbiased representation of each land cover class in the RF models, we divided the extracted pixel information for each class into stratified random samples of 70% and 30% of the pixels for training and testing the models, respectively. The image classification was performed using the random forest package [63] of the R programming software, version 3.4.3.

Although poor predictions have been reported for RF [64], the algorithm has performed better in cropland [25] and mangrove vegetation [53] classification. Thus, the performance of the RF algorithm may vary with different landscapes and cropping systems. For example, Loosvelt et al. [57] observed high classification uncertainty for mixed pixels at the heterogeneous boundaries of internally homogeneous cropping fields. Similarly, Van Tricht et al. [60] reported low classification accuracies at such field boundaries. Mixed cropping systems are very common in moist tropical landscapes. However, reports on the processing and use of SAR images for the mapping of heterogeneous tropical cropping land, such as perennial agroforests, are scarce.

3.3.3. Post-Processing and Classification Uncertainty Assessment

In remote-sensing-based mapping, the validity and reliability of classified maps are often determined on the basis of the estimated overall accuracy and the kappa coefficient [25]. Such values such as the user’s and the producer’s accuracy are prone to errors and uncertainties [65]. As a soft classifier, however, the RF algorithm provides us with the possibility to assess data- and computation-related uncertainties [57]. In our analysis, we used the user’s accuracy (the omission error), the producer’s accuracy (the commission error), the overall accuracy, and Kappa statistics, which compare the results of a chance classification versus our RF model’s accuracy. However, the pixel-based classification methods are prone to uncertainties from the use of unreliable data [65]. The RF algorithm, as a soft classifier, provides a vector (P_u) of classification probabilities for each image pixel: P_u = P₁, P₂, P₃,……, P_n for a classification with n categories, where P_i denotes the probability of belonging to class i (Table 2).

In this study, in addition to OOB error estimation, we evaluated the classification uncertainties of RF models using the maximum classifier probability (U) and a weighted uncertainty measure: the Shannon entropy (H) [66,67]. These uncertainties were calculated as:

$$U=1 - P_{max}$$

(7)

$$H=-{{\displaystyle \sum }}_{i=1}^N{P}_i·\log P_i$$

(8)

where P_i is the probability of belonging to class i, P_max is the maximum probability for a pixel’s class, and N is the total number of classes considered in the analysis.

The maximum probability class assignment by the soft classifier for a pixel does not always result in the true class label being assigned to the pixel. Thus, considering the entire range of values in a pixel’s probability vector H, compared to U, which only makes use of P_max, provides a more robust measure of uncertainty. It has a maximum value at the highest entropy: equal probability for all considered classes.

Loosvelt et al. [57] showed that H is reliable for evaluating the uncertainties in mapping croplands from SAR images. However, our study area is characterized by heterogeneous cropping systems and is located in a tropical landscape (Figure 2). For the best-performing RF models, based on the kappa accuracy, we computed and analyzed the U and H uncertainties for the classified maps and the considered land cover classes in the study area (Table 2). The estimations of uncertainty and the analysis were conducted in the Spyder Integrated Development Environment (IDE) of the Anaconda distribution for the Python programming software version 3.0 (Anaconda 3). The overall image-processing workflow is shown in Figure 5.

4. Results

All of the RF models yielded classification accuracies above 70%. The classification error and the sensitivity in discriminating land cover classes were different for each model. The models with the highest classification reliability, in increasing order of importance, were RE1, GLI3, and GL3.

4.1. Classification Accuracy

Table 4. The classification accuracies of different feature models based on the Random Forest (RF) classifier algorithm.

Model	Overall Accuracy (OA)% (95% CI)	Kappa	Out-of-bag (OOB) Error %
GLI1	78.80 (75.85, 81.53)	0.738	19.66
RE2	80.15 (78.76, 81.48)	0.757	19.46
RE1	81.04 (79.68, 82.35)	0.769	19.18
GL2	81.65 (78.74, 84.32)	0.773	18.47
GL1	82.74 (80.02, 85.23)	0.787	17.12
GLI2	82.97 (80.21, 85.48)	0.789	18.71
GLI3	85.07 (82.42, 87.47)	0.817	13.69
GL3	88.07 (85.52, 90.31)	0.853	12.85

summarizes the classification results for all eight RF models. All models had an overall accuracy (OA) of above 70%. However, compared to using VV or VH bands separately, the use of both co- and cross-polarization bands (GL3) resulted in the highest classification accuracy. The GL3 model had the highest overall accuracy (88.1%) and kappa (0.85); and, compared to other models, the OOB error estimate was the smallest (12.8%). The classification from the multispectral optical image (the RE1 model) had an overall accuracy of 81.1%; however, it featured a lower kappa (0.76769). Compared to the GL3 model, an OOB error difference of +7% was observed for RE1. Thus, the GLCM textures can be considered reliable for discriminating land cover/land uses. Considering the heterogeneous and dynamic vegetation in the landscape, improved feature selection using the GLCM approach is necessary to reduce class uncertainties.

The thematic land cover maps from the RE1 and GL3 models are shown in Figure 6. Separately, both the VV and VH GLCM-derived texture measures were poor in the prediction of non-vegetated land cover, and more so when both bands were included in the same model (Figure 6b). When included as input layers, the SAR backscatter intensity did not improve the classification accuracy. Similarly, the inclusion of vegetation indices from the multispectral optical image, taken during a dry season, did not improve the classification accuracy (Figure 6a). The texture measures from both VV and VH backscatter provide comparable, and potentially complementary, LULC mapping accuracy to the commonly used vegetation indices.

A visual analysis of the RE1 map revealed a relatively intact and continuous expanse of transition forest patches (Figure 6c). In contrast, the classified map from GL3 revealed that cocoa agroforests have fragmented the transition forest cover into smaller patches (Figure 6d). In addition, from the classification reliability estimates (100: commission error) in Figure 6b, the RE1 model was found to be more reliable in delineating non-vegetation land features. The SAR-based texture images had a high degree of reliability in delineating vegetation landscape features (Sv, Savannah; Af, Cocoa Agroforest; Fa, Subsistence farms; and Sf, Transition forests). Thus, although the multispectral optical image had better classification performance for the prediction of land cover classes in general, it was less reliable in discriminating perennial agroforests from transition forest land cover.

4.2. Uncertainty in Discriminating Vegetation Land Cover

The classification results from the RapidEye multispectral optical image RE1 had low overall and class uncertainties. From the cumulative estimates of the class probabilities shown in Figure 7, the classification uncertainty from the RE1 model converges at a probability of around 0.6 for both U and H. The uncertainty from the GL3 map converges at higher probabilities: 0.7 and 0.9 for U and H, respectively (Figure 7a). About 90% of the pixels classified by RE1 had H uncertainties below 0.4, compared to about 50% of the pixels for GL3. This difference is less obvious in the cumulative plot of U. Thus, the difference in uncertainty between RE1 and Gl3 was better revealed by Shannon entropy or H uncertainty.

The individual class uncertainties are compared in Figure 7. Although the classified map from the multispectral image (the RE1 model) had a lower accuracy, the class uncertainty was, compared to other land cover types, high for perennial cocoa agroforests and transition forest cover (Figure 7b). In comparison to RE1, the multiseasonal SAR image textures from the GL3 model had a high overall uncertainty in the prediction of a pixel’s class. However, perennial agroforests and transition forests were discriminated with a relatively lower individual class uncertainty (Figure 7c). The median of the class uncertainties was in a range between 0.2 and 0.4, which is comparable to that obtained from the single-date multispectral image (RE1). The uncertainties in land cover/land use discrimination by RE1 may reflect the vegetation status/phenology (canopy greenness). In the GL3 model, the volume scattering of the radar signal reflects changes in the water content and structure of different vegetation canopies.

4.3. The Contribution of Pixel Depth to Texture Feature Extraction

The SAR backscatter intensity images have a pixel quantization of 32 bits. The likelihood of grey level co-occurrences was lower in such a high dynamic pixel range, as observable from the unclassified pixels (the Un class) in Figure 6c. To reduce the classification uncertainty from the SAR image textures, we computed and compared three different image pixel quantizations or grey levels: 32 bits (GL3, as in the original SAR intensity image), 8 bits (GL3_B8), 6 bits (GL3_B6), and 4 bits (GL3_B4). The grey level computation was achieved using the Geospatial Data Abstraction Library (GDAL) in Python program.

The contribution of grey levels to the GLCM feature co-occurrence and classification probability is shown in Figure 8 and

Table 5. The classification results from different grey level GLCM models and land cover/land use surface area estimates. See Appendix A for details of class errors.

Model	Overall Accuracy (OA) % (95% CI)	Kappa	OOB Error %
RE1	81.04 (79.68, 82.35)	0.769	19.18
GL3	88.07 (85.52, 90.31)	0.854	12.85
GL3_B8	88.23 (85.74, 90.41)	0.854	11.84
GL3_B6	88.83 (86.48, 90.90)	0.862	9.92
GL3_B4	88.86 (86.50, 90.94)	0.862	10.38

. Pixel co-occurrence was low at high pixel dynamic ranges (Figure 8a), resulting in a large clustering of features with low and no predictions in Figure 8a,d, respectively. The dynamic pixel range of 64 grey levels significantly reduced the prediction error (Figure 8e). The class prediction probability was improved to a maximum of 0.99 in GL3_B6, with a difference of 7% from the value (0.92) recorded in GL3. The prediction did not improve with a further reduction in grey levels to 4 bits (Figure 8d). Optimizing the image pixel depth resulted in a marginal improvement in classification accuracy to a kappa value of 0.86862 (Table 5). The OOB error of prediction was reduced remarkably, from the initial 12.8% (GL3) to 9.9% in GL3_B6. These results show that the dynamic pixel range was vital in feature selection (the co-occurrence between pixels) and texture feature extraction by the ensemble algorithm.

The improvements in the predicted land cover/use maps, after optimizing the dynamic pixel ranges, are illustrated in Figure 9 and Figure 10. Figure 9 shows the predicted maps for the different texture-based models. Compared to the estimated 647.3 ha of unclassified area (Figure 9a) from GL3, at a pixel depth of 32 bits, the use of the 6-bit dynamic pixel range reduced the unclassified land area to about 7.4 ha (Figure 9c). See

Table A1. The predicted surface area (ha) of land cover/land use for the optical multispectral image and the GLCM textures at four different grey levels of quantization.

LULC Class	Model
	RE1	Gl3	Gl3_B8	Gl3_B6	Gl3_B4
Bu	68.89	547.38	650.35	545.89	719.07
Es	160.73	73.41	235.14	83.14	76.15
Sv	4485.2	3118.63	3035.37	3465.95	3629.72
W	137.08	0.27274	0.15079	0.16927	0.16446
Af	2733.63	2986.84	2904.18	3355.87	3254.77
Fa	2052.11	2787.95	3081.47	3120.90	3002.42
Sf	1706.94	1202.86	1306.70	1210.16	1094.61
Un	0	647.32	562.62	7.42	12.60

in the Appendix for more details.

The model validation results, by the Shannon Entropy estimates, are shown in Figure 10. The difference in class uncertainty between GL3 and GL3_B6 is most evident for the following classes: Built up (Bu), Cocoa agroforests (Af), subsistence farms (Fa), and transition forests (Sf). GL3_B6 had a comparably lower class error for these classes: 16.4%, 4.4%, 5.3%, and 6.1%, respectively (see

Table A2. The classification confusion matrix of the GL3_B6 Model. OA, Overall Accuracy; UA, User Accuracy; PA, Producer Accuracy.

nTrees = 550, mTry = 8, OOB error = 9.9%, OA = 88.8%
		Reference							Class Error	PA
		Bu	Es	Sv	W	Af	Fa	Sf
Predicted	Bu	244	2	5	0	23	17	1	0.164	83.6
	Es	14	89	19	0	6	21	0	0.403	59.7
	Sv	0	0	424	0	0	12	0	0.028	97.2
	W	2	0	3	0	6	7	0	1.000	0
	Af	3	0	0	0	325	8	4	0.044	95.6
	Fa	1	0	18	0	3	415	1	0.053	94.8
	Sf	0	0	0	0	16	0	246	0.061	93.9
UA		92.4	97.8	90.4	0	85.8	86.5	97.6

in the Appendix). Remarkably, the class uncertainty estimate for Sf is low and comparable to that for GL3_B6 and RE1. Unlike the RE1 model, with no difference in class uncertainty between Af and Sf, the significant difference between their class uncertainties in GL3_B6 confirms the reliability of the model to discriminate these two vegetation categories.

The estimated land cover/use against the reference model (RE1) is shown in Figure 11. The total area for each class is summarized in Table A1 (see the Appendix). A signification reduction in unclassified area is shown in Figure 11c,d. The transition forest area for all models was lower than that estimated by RE1. However, the class error and uncertainty for GL3_B6 were comparable to that of RE1 (Figure 10). The land area for the subsistence farming, cocoa agroforest, and built up classes was remarkably larger in the predicted texture-based maps.

5. Discussion

Land use and land cover (LULC) classification using SAR data often entails image filtering or multi-look preprocessing to reduce the speckle noise in SAR images [60]. However, such preprocessing reduces the resolution of the images. Considering the landscape structure and the inherently heterogeneous vegetation categories, we did not consider image speckle to be noise. We computed a temporal average of the measured textures across seasons, which presumably reduced any potential noise from individual image pixels. Meanwhile, the seasonal differences in volume scattering over the vegetation cover provided the texture information for discriminating the vegetation types. Such averaging, however, is less appropriate for the mapping of non-vegetation land cover; notably, water cover. The low classification performance for water surfaces can be explained as follows: variability in water cover is a result of the occurrence of seasonal swampy areas, the seasonal conversion of some swampy areas into subsistence farms of adapted crops, and the low backscatter intensity over water, which culminate in a low number of, or no, co-occurring pixels in the GLCM. Consequently, the GLCM has a low likelihood of extracting the texture features for the water class. Nonetheless, this land cover type was not the main interest of this study.

Compared to a “business as usual” classification using a single-season multispectral optical image (RE1), a combination of textures from both the VV and VH bands and the 6-bit grey level of quantization prior to the GLCM texture classification had the highest OA (88.1%) and kappa (0.85). Moreover, this accuracy resulted in a 3% and 9.3% reduction in prediction error over the GLCM texture at the default 32-bit grey level and the optical image, respectively. Notwithstanding the relatively high prediction error for the RE1 model compared to GL3_B6, the class prediction error was low for the built up, bare soil, savannah, and subsistence farmland land cover classes. These have rather distinctive optical spectral signatures. For vegetation cover with high tree canopies,

Table A3. The classification confusion matrix of the RE1 Model.

nTrees = 250, mTry = 2, OOB error = 19.2%, OA = 81.0%
		Reference							Class Error	PA
		Bu	Es	Sv	W	Af	Fa	Sf
Predicted	Bu	1238	0	0	0	0	1	0	0.001	99.9
	Es	2	529	8	0	0	117	0	0.194	80.6
	Sv	0	0	1469	1	52	172	2	0.134	86.6
	W	0	0	9	52	16	0	2	0.342	65.8
	Af	0	0	64	8	855	31	335	0.339	66.1
	Fa	0	57	233	0	55	1447	3	0.194	80.6
	Sf	0	0	3	0	309	4	665	0.322	67.7
UA		99.8	90.3	82.3	85	66.4	81.7	66.0

shows the high degree of confusion between cocoa agroforests and transition forests, confirming the reported challenges to discriminate them using their spectral signatures [2]. For such vegetation classes, the classification performance from optical data may reflect the vegetation status in the phenological cycle. Thus, spectral information from the dry season was less distinctive for cocoa agroforests versus transition forests.

The confusion between cocoa agroforests and transition forest was low compared to other classes (Table A2). This indicates that optimizing the grey level improved the classification and helped to distinguish the vegetation classes with a highly heterogeneous canopy. In the study landscape, the average range of backscatter intensity for both the VH and VV bands was 34 dB. The 6-bit grey level quantization indeed reflected this range. Therefore, the GL3_B6 model was optimal for discriminating different vegetation types, particularly those featuring a high canopy, and supports recommendations on the same pixel properties [38,39,40]. Other studies on heterogeneous cropland mapping featured an accuracy of 71% using C-band SAR intensity images [25]. In terms of OA, our result is in line with the accuracy observed in different heterogeneous cropping landscapes using a combination of C-band SAR and optical data [25,33]. However, the authors of [25,33] mapped cropping lands with an inherently homogeneous canopy. The landscape in this study is, however, tropical and dominated by vegetation and cropping fields with an internally heterogeneous canopy. Although no speckle filtering was applied to the used SAR images, the reduction of image grey level is a kind of smoothing filter over GLCM, which removes “noise” (between pixel pairs) and ensures pixel co-occurrences for estimate texture features. The “noise” removal reduced classification uncertainty of discriminating cocoa agroforests from other vegetation cover types. Similarly, application of grey level smoothing of medical images, natural and Magnetic Resonance (MR) images, improved classification accuracy while reducing the computational costs [68].

The texture-based land cover map shows spatial fragmentation of forest cover by cocoa agroforest land use. The forest area, estimated to cover 1706.9 ha (mean H = 0.25, class error = 0.32) using the RE1, was 500-ha less (mean H = 0.28, class error = 0.06) than the estimates from GL3_B6. This difference was consistent with field observations. The remnant transition forest patches are mostly owned by families and community groups, are used for hunting and performing traditional rituals, and serve as potential cocoa agroforest parcels. However, the land cover estimates for each model were made under different degrees of uncertainty. As expected, the GL3_B6 model discriminated the two vegetation cover types with significantly different class uncertainties: the mean class uncertainty for transition forests (0.28) was significantly different from the average class uncertainty for cocoa agroforests (0.4). The uncertainty estimates for the optical data, i.e., the RE1 model, confirmed the reported similarity in canopy structure of cocoa agroforests and transition forests [16]. Although the uncertainty for cocoa agroforests was significantly reduced by the GL3_B6, their relatively high class entropy compared to other land cover classes can be explained by the following management characteristics that influence their geometric properties: (1) the high variability in tree density and canopy structural components in cocoa agroforest plantations; and (2) the diversity in farm subsistence crop types and stages. In savannah lands, there is also a frequent change in canopy volume and structure along with phenological cycles. This is not the case in areas that have been converted to subsistence farming or previous farms that have been left to fallow. These biophysical changes cause unavoidable instances of classification confusion to occur [60]. Therefore, pixel-level uncertainties may not be easily reduced in a texture-based classification of SAR data. Following the high overall accuracy and a corresponding low individual class uncertainty, the multidate texture information from SAR images was, therefore, reliable as a classifier input for discriminating cocoa agroforestry land cover from transition forest land cover.

Classification validation based on such accuracy metrics as the overall, the user’s, and the producer’s accuracy is influenced by the sample class distribution in the training data [69]. The differences in pixel resolution between the image types (5 m for RapidEye and 10 m for SAR C-Band) resulted, therefore, in different numbers of reference pixels for the respective land cover types. This can be considered as an explanation for the high overall accuracy for the SAR images. The size of the training dataset was not as influential in classification validation by estimates of entropy or uncertainty. Entropy is a classic metric for biodiversity in ecology; however, in this study, it described the likelihood of a pixel belonging to either one of the considered classes [57]. Unlike the classification uncertainty reported for homogeneous crop fields [57], the U and H uncertainty trends were not comparable. For example, for a cumulative proportion of pixels classified with an uncertainty value of 0.4, the entropy (or H), compared to U, revealed a considerable difference between the RE1 and GLE models (Figure 7a). This can be explained by the heterogeneous vegetation classes and landscape. The slope of the cumulative Shannon entropy (H) was linear for RE1 and sigmoidal for GL3. The former reflects the high likelihood of a single land cover class for each pixel. Conversely, the sigmoid curve represents a complete range of probabilities for each image pixel: a soft boundary between land cover classes. Compared to perennial agroforests, the low class uncertainty for transition forests explained the high temporal stability of radar volume scattering over forest cover [30].

This study contributes to the application of the GLCM in SAR image processing for mapping tropical croplands with heterogeneous and multistrata canopies. We highlighted the important consideration of grey level quantization in the satellite remote-sensing-based mapping of heterogeneous land cover/land use. Depending on the landscape structure and vegetation types, the application of SAR imagery for cropland mapping may, however, warrant a different procedure [25,33,60].

A number of amendments to our used procedures can be proposed. The texture measures were preselected following the recommendation in Reference [31], and were based on inferences from optical image classification results and tailored for vegetation mapping. This preselection of texture features may be biased towards non-vegetation cover. We estimated the GLCM textures using a 5 × 5 window scale considering a 10-m pixel resolution for Sentinel-1 SAR. However, different window sizes may influence the estimates of texture features and the classification’s accuracy. From this perspective, we assume that a fusion of texture information with optical images [25,33] at a temporal scale may decrease the uncertainty in discriminating vegetation classes. Moreover, considering the diversity of vegetation across cocoa production landscapes and the complex structure and canopy of cocoa agroforests, the continuous archiving of, and improvement in, SAR image resolution can help to train deep-learning procedures for land cover/land use classification at the scale of several landscapes. This, however, will inevitably entail a considerably large collation of training data.

6. Conclusions

This study initiated an application of multipolarization and multi-temporal C-band SAR data for the discrimination of cocoa agroforest cropping land in a multi-use tropical landscape. We used seasonal differences in volume scattering from the dielectric (water content) status and structure of the vegetation canopy as a metric to discriminate the vegetation types. We summarize our conclusions as follows:

• For the same window size and an invariant direction, reducing the SAR image’s grey level or quantization marginally improved the texture-based classification accuracy, but significantly reduced the uncertainty in discriminating cocoa agroforests from other vegetation cover types.
• The classification validation using Shannon entropy (H) estimates revealed subtle differences in individual class prediction and provided reliable information for drawing inferences about the vegetation structure in a multi-use and heterogeneous landscape.
• The magnitude of forest fragmentation by cocoa agroforests, which is concealed by vegetation indices from spectral reflectance, can be reliably mapped using texture measures from C-band SAR images.

The current use and evaluation of Machine Learning (ML) algorithms for land cover classification is directed towards achieving higher accuracy [54]. While they may be reliable in object-level classification experiments, accuracy results alone are less informative for understanding the non-thematic characteristics of land cover, such as the variation in vegetation structure. Using an algorithm (RF) that provides classification uncertainty at the pixel level, our measure of uncertainty is especially relevant in the context of mapping agroforestry land cover/land uses that have a heterogeneous canopy. Thus, while ML has increasing potential for remote-sensing-based land cover classification, the choice of algorithm for mapping vegetation may consider other metrics of the land cover structure.

This study addressed the problem of optimizing the mapping of cocoa agroforest land cover. While the research was conducted in Cameroon, it is likely that similar landscapes exist elsewhere in humid tropical Africa, notwithstanding that the structural components might be different. For instance, the crop component might be different from cocoa (e.g., coffee or other tree crops) or the amount and structure of natural forest remnants might be different. Hence, while this study addressed one particular landscape case, we assume that it represents a feasible starting point for validation in other similar tropical climatic zones featuring agroforestry land cover with an inherently heterogeneous vegetation composition. Therefore, our approach can be considered generic in the sense that it can contribute to estimation of the contribution of agroforestry to national and regional REDD+ (Reducing Emissions from Deforestation and Forest Degradation, and the role of the conservation, sustainable management of forests, and the enhancement of forest carbon stocks) strategies. However, there is a need to assess classification uncertainties in different agroforestry-dominant landscapes for operational regional mapping, such as in the Congo Basin sub-region.