E-Agriculture Planning Tool for Supporting Smallholder Cocoa Intensification Using Remotely Sensed Data

Abstract

Remote sensing approaches are often used to monitor land cover change. However, the small physical size (about 1–2 hectare area) of smallholder orchards and the cultivation of cocoa (Theobroma cocoa L.) under shade trees make the use of many popular satellite sensors inefficient to distinguish cocoa orchards from forest areas. Nevertheless, high-resolution satellite imagery combined with novel signal extraction methods facilitates the differentiation of coconut palms (Cocos nucifera L.) from forests. Cocoa grows well under established coconut shade, and underplanting provides a viable opportunity to intensify production and meet demand and government targets. In this study, we combined grey-level co-occurrence matrix (GLCM) textural features and vegetation indices from Sentinel datasets to evaluate the sustainability of cocoa expansion given land suitability for agriculture and soil capability classes. Additionally, it sheds light on underexploited areas with agricultural potential. The mapping of areas where cocoa smallholder orchards already exist or can be grown involved three main components. Firstly, the use of the fine-resolution C-band synthetic aperture radar and multispectral instruments from Sentinel-1 and Sentinel-2 satellites, respectively. Secondly, the processing of imagery (Sentinel-1 and Sentinel-2) for feature extraction using 22 variables. Lastly, fitting a random forest (RF) model to detect and distinguish potential cocoa orchards from non-cocoa areas. The RF classification scheme differentiated cocoa (for consistency, the coconut–cocoa areas in this manuscript will be referred to as cocoa regions or orchards) and non-cocoa regions with 97 percent overall accuracy and over 90 percent producer’s and user’s accuracies for the cocoa regions when trained on a combination of spectral indices and GLCM textural feature sets. The top five variables that contributed the most to the model were the red band (B4), red edge curve index (RECI), blue band (B2), near-infrared (NIR) entropy, and enhanced vegetation index (EVI), indicating the importance of vegetation indices and entropy values. By comparing the classified map created in this study with the soil and land capability legacy information of Bougainville, we observed that potential cocoa regions are already rated as highly suitable. This implies that cocoa expansion has reached one of many intersecting limits, including land suitability, political, social, economic, educational, health, labour, and infrastructure. Understanding how these interactions limit cocoa productivity at present will inform further sustainable growth. The tool provides inexpensive and rapid monitoring of land use, suitable for a sustainable planning framework that supports responsible agricultural land use management. The study developed a heuristic tool for monitoring land cover changes for cocoa production, informing sustainable development that balances the needs and aspirations of the government and farming communities with the protection of the environment.

1. Introduction

The Autonomous Region of Bougainville (Bougainville) is the easternmost island of Papua New Guinea (PNG). The economy depends on subsistence agriculture, cash crops, including coconut and cocoa, and small-scale mining (https://bougainvillenews.files.wordpress.com/2018/06/bougainville-strategic-plan-2018-2022.pdf (accessed on 19 July 2021)). To meet the ambitious targets set by PNG’s government [1], opportunities for cocoa production need a kickstart by establishing an E-agriculture strategy (E-Ag strategy) that promotes balanced land use for sustainable economic, social, and environmental development [2]. This renewed focus on cocoa production should prioritise intensification, aiming for increased yield per hectare rather than expansion through deforestation [3]. Cocoa thrives under established coconut palm shade, making cocoa expansion within existing coconut plantations a viable intensification opportunity. Further, it has been reported that coconut plantations can be distinguished [4] from forested areas with ease, unlike cocoa smallholder orchards that are difficult to distinguish from other vegetation covers, such as natural forest [5]. Therefore, we hypothesise that remote sensing can help distinguish areas where cocoa is established or where it can potentially be grown by classifying areas under coconut palm. In this manuscript, the coconut–cocoa areas will be referred to as potential cocoa regions or orchards for consistency.

Numerous studies have investigated land cover classification at local, regional, and global scales using remote sensing (RS) platforms to extract valuable information and capture land cover changes. For large-scale crop type mapping, studies have focused on using the Moderate Resolution Imaging Spectroradiometer (MODIS) [6,7,8,9,10], which provides high temporal resolution and potential daily monitoring of land use changes. However, the spatial resolution of MODIS (250–1000 m), having the potential to provide the areal coverage required [11,12], may not be sufficient to resolve the expected land use changes in smallholder areas (<1 ha) due to potential errors caused by pixel mixing between different classes and land cover heterogeneity.

To address the aforementioned issues, researchers have used data retrieved from satellites with higher spatial resolution, such as Landsat (30 m) [13,14,15], or Sentinel-1 and Sentinel-2 (10 and 20 m) [16,17,18,19,20]. In addition to improved resolution, Landsat and Sentinel-2 provide visible and near-infrared (NIR) data that better interpret differences within an area. Sentinel-2 offers three bands in the red-edge region, which are crucial for crop type mapping [21].

Forkuor et al. [21] showed that Sentinel-2 improved the classification accuracy in rural Burkina Faso by 4 percent compared to Landsat bands. Persson et al. [22] found that the red-edge and short-wave bands of Sentinel-2 were particularly effective in discriminating tree classification. Therefore, land cover classification can be enhanced by utilising the spectral information in the red-edge region of the spectrum provided by Sentinel-2, along with its high spatial resolution compared to other platforms like MODIS and Landsat. Combining Sentinel-1 and Sentinel-2 may further improve these enhancements [23,24].

Moreover, synthetic aperture radar data can help overcome data noise caused by clouds and shadows in surface reflectance data [25]. More significantly, the application of grey-level co-occurrence matrix (GLCM) textural features and indices to very high-resolution imagery can address the issue of canopy heterogeneity [26,27,28,29,30]. This approach has proven successful in differentiating between coconut plantations and native forests [29].

Additionally, this study offers an opportunity to assess the relationship between regions where cocoa can be grown and the main soils of Bougainville, as well as agricultural land capabilities, using legacy information. The current study focuses on monitoring changes in the land area used for coconut palm plantations and with the potential to sustain cocoa production, by combining remote sensing information, machine learning algorithms and agricultural capability legacy data. Generally, cocoa grows well on deep, well-drained, and reasonably fertile soil [31]. The maps created in this study also have the potential to contribute to report on the Sustainable Development Goals of the United Nations (UN-SDGs) and the E-agriculture initiative implemented by the PNG government.

Specifically, this paper aims to:

(a)

Investigate cocoa intensification opportunities by classifying areas under coconut palm as proxies for locations where cocoa is already grown or can be cultivated.

(b)

Explore cocoa expansion opportunities based on highly suitable and capable areas for agricultural production that are not currently utilised.

(c)

Examine the sustainability of cocoa production by comparing the resulting potential cocoa maps with legacy information, that is, soil/land capability/cocoa suitability maps.

(d)

Explore the potential use of this tool as a land cover change monitoring tool for measuring deforestation in the future, particularly as evidence to address the PNG E-Ag strategy and UN-SDGs.

2. Materials and Methods

2.1. Study Area

The study was conducted in the Autonomous Region of Bougainville, Papua New Guinea, which covers an area of 9318 km². This region was chosen based on its potential for Theobroma cocoa L. (Cocoa) production. A significant area in Bougainville is used for agriculture, with cocoa and coconut being the primary cash crops [2,32]. Therefore, we anticipate observing a relation between land capability classes and the cocoa orchard potential in the classified map. Figure 1 shows the elevation, mean annual temperature, and annual rainfall for Bougainville obtained from WorldClim version 1 [33].

2.2. Datasets Used and Preprocessing

This study used two main datasets for the detection of potential areas for cocoa orchards. These datasets included the tiles of synthetic aperture radar (SAR) from Sentinel-1 and surface reflectance data from Sentinel-2 (Sentinel-2 MSI, Level-2A) satellites [34], which were processed through Google Earth Engine [35]. Ancillary data, such as a terrain slope obtained from the processing Shuttle Radar Topography Mission (SRTM) at a spatial resolution of 30 m, was also incorporated into the study [36]. The Sentinel datasets were filtered to limit their extent to the study region and the date range from 2019 to 2020.

Surface reflectance images that contained more than 20% of clouds, stored in the “cloud pixel percentage” property of the collection, were filtered out. An additional cloud mask was applied to the collection by appending cloud probability scenes of Sentinel-2, processed through the s2cloudless library [37], and setting a threshold of cloud probability to 30% [38]. Due to the high occurrence of clouds in the study region (Figure 2), the masked surface reflectance collection was averaged. After processing the Sentinel-2 data, raw bands (B2-B8A and B11) were extracted and used as covariates to calculate vegetation indices, which were combined with the raw bands. Six vegetation indices were calculated for different calendar years using the visible, near-infrared (NIR), and red edge bands. Table 1 shows the vegetation indices used, which can capture specific vegetative features, such as leaf greenness, tree height, and systematic cocoa block planting, compared to a forest that may contain vegetation of various shapes and sizes.

The vegetation indices include the enhanced vegetation index (EVI), the normalised difference moisture index (NDMI), and the chlorophyll index green (Chlgreen), which are visible NIR-based indices. Additionally, The Vogelmann Red Edge Index (VREI), the chlorophyll red edge index (CIr), and the red edge curve index (RECI) were used. The RECI was derived to describe the vegetation curve in the region between the red edge band 1 (RE1; 698 nm) and the NIR (899 nm). The SAR data were filtered to a single vertical transmit and vertical receive (VV) co-polarisation to maximise the number of images in the region. The SAR data were also filtered to the Interferometric Wide Swath (IW) mode of operation, which is the primary acquisition mode that leads to images of a 250 km swath at 5 m by 20 m resolution. Furthermore, to remove edge effects, a mask was applied to remove all pixels with less than −30 dB from the collection. As coconut plantations (potential cocoa orchards) exhibit similar spectral patterns to forests but display a strong homogeneity in terms of vegetation growth, the grey-level co-occurrence matrix (GLCM) texture features analysis was employed [39]. The GLCM with a square kernel of three pixels was applied to distinguish between forests and potential cocoa orchards. From this analysis, angular second moment (VVasm), variance (VVvar), inverse different moments (VVidm), and entropy (VVent) bands were selected and used in the detection.

To quantitatively and qualitatively evaluate the results and assess the sustainability of potential cocoa expansion, a comparison was conducted using independent datasets of cocoa orchards [40], soil series in Bougainville [41], the cocoa suitability map from Papua New Guinea [41], and land use capability maps [42]. These were compared with the map of potential cocoa orchards developed in this work (refer to Section 2.3).

The cocoa suitability map [41] identified cocoa suitability classes in PNG using a multicriteria decision analysis, considering three biophysical and two market access criteria (Figure 3). This map was used to establish a proxy for sustainability evaluation by comparing the areas where cocoa can be grown under coconut palm plantations with the identified suitability classes.

The land use capability map for Bougainville in Figure 4 shows the areas suitable for cultivation based on various aspects describing soil, physiography, access limitations, and hazards. The land suitability in this map ranges from class I to class V, indicating lands unsuitable for cultivation to lands exceptionally suitable for any commercial use, respectively.

The soil map of Bougainville (Figure 5) presents the major soil types (http://worldmap.harvard.edu/data/geonode:DSMW_RdY (accessed on 23 April 2021)) found within the province. We ranked the soil types based on their capability for cocoa production using the following principles:

• Acrisols with a clay-rich subsoil (present in southeast Bougainville), were ranked as having low suitability for cocoa due to acidic conditions, low fertility, and the need for management.
• Luvisols with a clay-rich subsoil (present in north Bougainville—Lonahan), corresponding to fertile soils, were ranked as having a moderate suitability.
• Other soil types categorised are:

  (a)

  Relatively young soils or soils with very little or no profile development, or very homogenous sands with moderate fertility, predominately Cambisols (present in central Bougainville), were ranked as having moderate suitability.

  (b)

  Soils strongly influenced by water, such as Fluvisols (present in south and central west Bougainville around Siwai), were ranked as being moderately suitable for cocoa as they require drainage management.

  (c)

  Soils of volcanic origin, that is, Andosols (main soil type of Bougainville spread across central and south), were ranked as exceptional because they are deep and present adequate physical soil properties for cocoa.

2.3. Detection of Potential Cocoa Orchards and Uncertainty Estimation

All the different rasters derived from Sentinel collections for each calendar year were stacked together with the slope of the terrain and the distance to roads to create a single multi-band image used for potential cocoa orchard detection (Figure 6). To achieve this, a random forest model with 250 decision trees was used. This algorithm was selected because it outperformed other classifiers, including Support Vector Machines (SVM), Classification and Regression Trees (CART), Gradient Tree Boost, and Naïve Bayes [43] (Supplementary Materials Table S1). The supervised classification was conducted by sampling anchor pixels from polygons drawn in two land use classes, including coconut orchards and others (such as water, roads, forest, and bare land). A total of 268 polygons were drawn from known potential cocoa orchards (coconut plantations) and randomly sampled based on pixels. The sampling order was randomised and limited to 20,000 samples. Additionally, 444 polygons from other land uses were drawn, sampled, randomised, and from these, another 60,000 samples were selected. These two categories were combined, randomised again, and split into training, validation, and testing subsets in proportions of 10%, 80%, and 10%, respectively.

Table 1. Vegetation indices used in this study.

Index	Equation	Source
Enhanced vegetation index (EVI)	$\mathrm{EVI}=\frac{2.5\left(\mathrm{NIR}-\mathrm{red}\right)}{\left(\mathrm{NIR}+6\times \mathrm{red}-7.5\times \mathrm{blue}+1\right)}$	[44]
Normalised difference moisture index (NDMI)	$\mathrm{NDMI}=\frac{\left(\mathrm{NIR}-\mathrm{SWIR}1\right)}{\left(\mathrm{NIR}+\mathrm{SWIR}1\right)}$	[45]
Chlorophyll index green (ChI green)	$\mathrm{ChI}\mathrm{green}=\frac{\mathrm{NIR}}{\mathrm{green}-1}$	[46]
Vogelmann Red Edge Index (VREI)	$\mathrm{VREI}=\frac{\mathrm{NIR}}{\mathrm{Red}\mathrm{Edge}}$	[47]
Chlorophyll index red edge (CIr)	$\mathrm{CIr}=\frac{\mathrm{Red}\mathrm{Edge}3}{\mathrm{Red}\mathrm{Edge}1}-1$	[48]
Red edge curve index (RECI)	$\mathrm{RECI}=\left(\mathrm{Red}\mathrm{Edge}2-\mathrm{Red}\mathrm{Edge}1\right)\times \frac{\mathrm{NIR}}{\mathrm{Red}\mathrm{Edge}1}$	Proposed in this study

Table 1. Vegetation indices used in this study.

Index	Equation	Source
Enhanced vegetation index (EVI)	$\mathrm{EVI}=\frac{2.5\left(\mathrm{NIR}-\mathrm{red}\right)}{\left(\mathrm{NIR}+6\times \mathrm{red}-7.5\times \mathrm{blue}+1\right)}$	[44]
Normalised difference moisture index (NDMI)	$\mathrm{NDMI}=\frac{\left(\mathrm{NIR}-\mathrm{SWIR}1\right)}{\left(\mathrm{NIR}+\mathrm{SWIR}1\right)}$	[45]
Chlorophyll index green (ChI green)	$\mathrm{ChI}\mathrm{green}=\frac{\mathrm{NIR}}{\mathrm{green}-1}$	[46]
Vogelmann Red Edge Index (VREI)	$\mathrm{VREI}=\frac{\mathrm{NIR}}{\mathrm{Red}\mathrm{Edge}}$	[47]
Chlorophyll index red edge (CIr)	$\mathrm{CIr}=\frac{\mathrm{Red}\mathrm{Edge}3}{\mathrm{Red}\mathrm{Edge}1}-1$	[48]
Red edge curve index (RECI)	$\mathrm{RECI}=\left(\mathrm{Red}\mathrm{Edge}2-\mathrm{Red}\mathrm{Edge}1\right)\times \frac{\mathrm{NIR}}{\mathrm{Red}\mathrm{Edge}1}$	Proposed in this study

To assess the uncertainty of the classification, the output of the random forest model was set to the probability mode to obtain probabilities for correct pixel-based classifications. These probabilities were used to evaluate the uncertainty using the normalised Shannon entropy [49]:

$$H\left[P\right]=\frac{-{{\displaystyle \sum }}_{i=1}^n{p}_{i \ln p_i}}{S_{max}}$$

(1)

where p_i is the probability of class i, n is the number of classes, and S_max corresponds to ln n. The normalised Shannon entropy ranges from 0 to 1, where 0 implies total certainty about the classification and 1 for a total uncertainty [50].

2.4. Validation of Classification

For this analysis, a holdout cross-validation analysis was subsequently applied to the sample subsets for the data corresponding to the calendar year 2019. This allowed the development of a confusion matrix, from which the performance was evaluated using metrics such as kappa, overall accuracy, f1-score, consumer’s accuracy, and producer’s accuracy [51]. The kappa metric compares the agreement obtained to the agreement obtained by chance [52]. The overall accuracy is the ratio between correctly classified samples and the total samples. The f1-score corresponds to the harmonic mean between precision and recall. The producer’s accuracy for each class is calculated by dividing the number of correctly classified samples by the total number of samples per column in the confusion matrix, while the consumer’s accuracy is the correctly classified samples per class divided by the total samples per row.

In simple terms:

(a)

The producer’s accuracy (relates to omission errors) represents how well reference pixels of the land-use type are classified.

(b)

The omission error refers to excluding a pixel that should have been included in the class (i.e., omission error = 1 − producer’s accuracy).

(c)

The consumer’s accuracy (relates to commission errors) represents the probability that a pixel classified into a given category represents that category on the ground.

(d)

The commission error refers to including a pixel in a class when it should have been excluded (i.e., commission error = 1 − consumer’s accuracy).

Since cocoa orchards are expected to have a multidecadal lifespan, all samples previously selected for detection through the random forest algorithm (training and validation from 2019) were combined and used as a validation source for the year 2020.

3. Results

The results showcase the classification conducted in this study involving three main components. First, the use of the finer resolution Sentinel-1 and Sentinel-2 imagery (Figure 6); second, the processing of imagery (Sentinel-1 and Sentinel-2) for feature extraction using all the 22 variables together; third, fitting the RF model to detect potential cocoa orchards. An attempt has also been made to find relationships between this classified map and soil and land capability legacy information. Therefore, the following subsections of the results report on the main components.

3.1. Classification Performance

The classification accuracy was above 90% for both the calibration and validation datasets (

Table 2. Classification accuracy for calibration and validation datasets.

Dataset	Kappa	Overall Accuracy	Producer’s Accuracy	Consumer’s Accuracy	F1-Score
Calibration	0.999	0.998	Cocoa region: 0.999	Cocoa region: 0.999	Cocoa region: 0.999
			Non-cocoa: 0.998	Non-cocoa: 0.999	Non-cocoa: 0.999
Validation	0.929	0.974	Cocoa region: 0.992	Cocoa region: 0.992	Cocoa region: 0.983
			Non-cocoa: 0.918	Non-cocoa: 0.977	Non-cocoa: 0.948

). The model provided a high overall calibration accuracy (0.99) as well as a high validation accuracy (0.97). This high overall accuracy needs to consider the imbalance between classes, as machine learning classification can be greatly influenced by imbalance in classes [53]. A decline in the more conservative Kappa accuracy for the validation model (0.929) still shows a very good performance in the classification. However, Kappa results should be taken with caution as they refer to an agreement by chance instead of the actual agreement, which may lead to a wide range of values, even if the overall accuracy is high, and difficulties in comparing Kappa coefficients [54].

The source of errors when evaluating the producer’s and consumer’s accuracy is mainly in the non-cocoa classified regions, which implies larger omission errors (false negatives). Furthermore, the commission error is higher for non-cocoa regions compared to potential cocoa regions. This is reflected in the lower consumer’s accuracy for non-cocoa regions (0.977) compared to cocoa regions (0.992), which also represents within-class errors in a map (false positives). These results, which imply the combination of different metrics (i.e., Kappa coefficient, overall accuracy, f1-score, producer’s accuracy, and consumer’s accuracy) show a reliable performance in the classification, even considering the particular limitations indicated, providing a robust calibration dataset of areas with cocoa potential. This can be used to assess future land changes, providing a tool utilisable by government policymakers responsible for monitoring and reporting on land cover changes.

Additionally, the resulting map was compared against the independent dataset of known cocoa orchards [40] (Figure 5). All the cocoa farms from the independent dataset are within the classified areas with cocoa potential, which is consistent with the general land capability map for Bougainville shown in Figure 4.

Following a robust optimisation process, the RF model created in this study shows classification reliability. In this study, the top five variables contributing the most to the RF model were B4 (red), RECI, B2 (blue), NIR entropy, and EVI (Figure 7). The B4 and B2 represent raw bands, which are strongly absorbed by vegetation. Similarly, EVI relies on red, NIR, and blue bands, and varies depending on changes in the canopy structure [55]. These covariates provide insights into the vegetation types. On the other hand, RECI is based on differences in the red edge range of the wavelength spectrum and the ratio between NIR and red edge bands. The red edge wavelengths exhibit high reflectance variations in vegetation [56] and are known to be good indicators of the vegetation physiological status [57]. Entropy in the NIR and red bands corresponds to a textural analysis of images as it uses the pixel neighbourhood (kernel) for calculation. When applied to these wavelengths, entropy helps distinguish between land uses with different reflectance characteristics. All surface reflectance-derived bands took precedence over SAR and ancillary data.

Although multi-collinearity does not impact the performance of random forest results [58,59,60], it can lead to unreliable feature rankings [61]. Therefore, Figure 8 presents the correlogram of covariates from the samples. Bands in the visible and red edge wavelength spectrum exhibit high correlation, and a similar pattern can be observed among indices calculated from reflectance band combinations. As a result, the importance of covariates must be taken with caution.

3.2. Detection of Potential Cocoa Orchards

The classified area with the potential for cocoa expansion (under coconut plantations, Figure 9a) accounts for 5% of Bougainville’s total area (479.97 km²), as assessed by the degree of uncertainty map, which ranges from 0.20 to >0.80 (Figure 9b). Some observations align clearly with the land capability map (Figure 4). Areas with high uncertainty (Figure 9b) correspond to regions with significant variability in land capabilities (Figure 4).

This relationship is further illustrated in Figure 10, where most lands with cocoa potential are located along the coast in the central region, in line with the identified land capability in Figure 4. Similarly, the limited or unsuitable land classification in the south (Figure 4) aligns with the regions of non-cocoa potential in Figure 10. There is a notable coincidence of cocoa potential in both maps for the western region (Figure 4 and Figure 10). Additional examples of good and poor detection for cocoa potential are provided in Figure 11.

This heuristic tool offers users the opportunity to assess the potential for cocoa expansion, evaluate land cover change, and detect deforestation.

Comparing the soil map with the classified map of cocoa potential (

Table 3. Cocoa grown in Bougainville on different soil types.

Soil Type	Soil Area (ha)	Area of Cocoa Grown by Soil Type (%)
Af-Ferric Acrisol	50,177.7	1.46
Bh-Humic Cambisol	68,292.1	0.14
Je-Eutric Fluvisol	98,743.7	8.34
LC-Chromic Luvisol	43,727.0	22.32
TM-Andosol	633,407.4	75.24

) and quantifying the extent of lands with cocoa potential in each soil type, it becomes evident that local farmers have successfully chosen where cocoa may be produced. The majority of land with potential for cocoa expansion have Luvisols and Andosols and a few with Acrisols. The distinction between the use of the Cambisols in the central part of Bougainville and the moderately suitable Fluvisols and Luvisols is likely influenced by limited road accessibility in the central region. From this, we understand that suitability to produce cocoa depends mainly on soil capability and market access constraints amongst other biophysical variables (e.g., climate).

4. Discussion

Native forests throughout the tropics face an increasing threat from human-induced land use changes and climate change [62,63] and there is the need to support the livelihoods of smallholders, including those growing cocoa [2,32]. Previously, it has been nearly impossible to accurately map cocoa orchard distribution, as cocoa trees are identified as forests in traditional satellite images [64]. However, regularly planted coconut-dominated canopies can be easily detected [29]. Coconut is a common shade tree in the cocoa production systems of PNG [65]. Therefore, this study has identified an opportunity for a systematic and broad-scale assessment of coconut to evaluate the areas that can be used for cocoa intensification.

We used freely available Sentinel-1 and Sentinel-2 satellite imageries and successfully identified the processing required to discriminate the coconut classification from the surrounding forest [29] and overcame the high frequency of cloud cover associated with the region. Our results demonstrate that combining remote sensing techniques using very high-resolution imagery, GLCM textural analysis, vegetation indices, and a RF classifier gave an overall validation accuracy of 97 percent for potential cocoa and non-cocoa classification. These results were not entirely a surprise given the huge body of successful work along similar lines [26,27,28,29,30]. Additionally, the confusion matrix helped reveal the sources of errors; based on this, it will be possible to improve the classification accuracy in the future, perhaps by adding more polygons/pixels from non-cocoa regions.

Further, we expect the top five covariates used in the model (B4, RECI, B2, NIR entropy, and EVI) to differentiate between systematically planted cocoa (4 m spacing) and coconut (12 m spacing) trees and their associated cropping cycle from the forest. For example, the red edge and NIR bands are probably able to capture information on planting patterns, the colour of leaves, and the height of cocoa (average 3–5 m) and coconut (average 15–25 m) trees from the forest [66]. We cannot undertake this analysis at this stage, and this will require some ground-based data to create a calibration dataset to distinguish cocoa from coconut trees. The land capability was mapped (Figure 4) in 1967 in Bougainville without the use of finer scale information currently available [42]. The map illustrates the general land suitability for cultivation across the island, but not to the scale needed for the current monitoring of intensification or expansion potential. While the 1967 map and our remote sensing effort cannot be compared directly due to differences in methodology, the general trend shows that the area under cocoa complements the regions that are suitable for cultivation, as shown by the coincidence of uncertainties (Figure 9 and Figure 10) with land capability (Figure 4). Therefore, we propose consulting the legacy data to assess the importance of land capability on present cocoa orchard dynamics [67].

Studies on land suitability [65,68] show that much is known about the capability of a given area to be able to grow cocoa trees in PNG. However, there are still uncertainties around many aspects of interactions between soil capability and conditions [30]. For example, the highly suitable Andisols of Bougainville may have low magnesium and calcium with high but unavailable phosphorous [69]. These are of the type 1 class (soils with no physical limitation to root development within 1.5 m of the surface) [68]. This shows that these soils can produce cocoa with additional nutrient management required to maintain or enhance soil conditions based on deficiencies [70] using fertilisers [30].

We found that most farmers already grow plantations where the soil is suitable, and these results can be used to identify local limitations for more targeted and effective soil management through training opportunities. The map can help understand the spatial extent for future intensification of cocoa farming and can be used as a tool for preemptive decisions of where extension services and demonstration sites for soil/land management can be directed to improve cocoa and coconut production and their value chain (https://bougainvillenews.files.wordpress.com/2018/06/bougainville-strategic-plan-2018-2022.pdf accessed on 19 July 2021).

For sustainable intensification purposes, an attempt has been made to investigate the suitability of areas where cocoa can potentially grow under coconut plantations. Specifically, this addresses the challenges outlined in the PNG E-Ag Strategy (https://png-data.sprep.org/dataset/papua-new-guinea-e-agriculture-strategy-20172023/resource/379a5244-5b98-4b72-95e6 accessed on 19 July 2021) developed in accordance with the framework proposed by the FAO ITU E-agriculture Strategy Guide (http://www.fao.org/3/a-i5564e.pdf accessed on 19 July 2021). This approach and the resulting data and interpretation presented here provide a monitoring and evaluation framework for land cover change and suitability for intensification. Further, we outline two main impacts of the study:

(a) As a land cover change monitoring tool: Comparing the 2019 and 2020 data, a preliminary analysis of the spatial distribution of potential cocoa regions indicates an increase in area under plantations within 1 year. For the 2019 classified map, the total area under crop production is 479.97 km² and an increase in the plantations by 7.31 km² is seen (

Table 4. A potential land cover monitoring tool. Area under cocoa in years 2019 and 2020 in Bougainville.

Year	Region	Area (km2)
2019	Non-cocoa	8893.16
	Cocoa	479.97
2020	Non-cocoa	8885.85
	Cocoa	487.28
Change 2019–2020	Non-cocoa	−7.31
	Cocoa	7.31

). This result shows an increase in plantations in Bougainville from the year 2019 to 2020. This land use change is not insignificant, since it represents a growth of 1.52% of plantation extent in one year. Changes of this magnitude should be monitored and managed to avoid damage to natural ecosystems. Therefore, the calibration dataset created in this study can be used as a cocoa expansion and monitoring tool to quantify short- and long-term changes in land cover.

(b) Sustainable cocoa intensification: Comparing the cocoa suitability map (Figure 3) to the classified map, we see that (

Table 5. Suitability class areas for cocoa orchards in Bougainville and the corresponding planted areas.

Cocoa Suitability	Area (ha)	Area of Cocoa in Bougainville (%)	Area of Cocoa within Class (%)
Suitability 1	0	0	0
Suitability 2	8738.7	0.24	1.16
Suitability 3	437,677.1	24.68	2.35
Suitability 4	413,485.0	69.21	6.97
Suitability 5	12,467.1	5.86	19.58

) plantations (coconut and coconut–cocoa) are currently grown in highly suitable areas. This is important because unsustainable development implies using land resources beyond their limits, and therefore, coconut–cocoa orchard expansion towards unsuitable lands can be used as an indicator of unsustainable development. Further, planned roads would also create new deforestation hotspots via the rapid expansion of logging, mining, and oil-palm plantations [71].

Since our results show that coconut and cocoa are already produced in highly and exceptionally suitable areas, this sets a limit to their potential expansion in a sustainable manner. Yet, these values should not be used as absolute. For instance, extremely suitable areas for cocoa production are under plantations in about 20% of their extent. This does not mean that the remaining 80% of these areas could be used for cocoa plantations. These regions may have excellent potential for cocoa, but they might be prioritised for other uses (e.g., residential, environmental protection, horticulture, etc.). Income diversification has been shown as a significant factor in smallholder wealth in Bougainville [2]. In this regard, economically sustainable agriculture must consider expansion limits to foster ecosystem functioning and biodiversity, which is commonly affected by agriculture, while intensification has proven insufficient on its own to cope with this issue [72].

Being an understory tree, cocoa’s intensification can be considered a sustainable cash crop in a multi-layered forest system [73,74], adopting agroforestry principles [75]. The mapping shows great potential for intensification as an understory crop with coconut. It is expected that the cocoa grown in this type of farming system holds enormous potential for environmental [75] and cultural [76] conservation in regions under pressure from climate change, deforestation, and mono-cropping [63,76,77,78].

5. Conclusions

This study combined grey-level co-occurrence matrix (GLCM) textural features and vegetation indices from Sentinel datasets to evaluate the sustainability of cocoa expansion in given land suitability classes for agriculture and soil capability classes. We mapped areas where cocoa smallholder orchards either already exist or can be grown by using fine-resolution C-band synthetic aperture radar and multispectral instruments from Sentinel-1 and Sentinel-2 satellites, respectively. Further, the processing of imagery (Sentinel-1 and Sentinel-2) for feature extraction using 22 variables was conducted as was fitting an RF model to detect and distinguish cocoa orchards from non-cocoa areas. We were able to differentiate potential cocoa and non-cocoa regions with 97 percent overall accuracy and over 90 percent producer’s and user’s accuracies.

Our planning tool provides a foundation for sustainable cocoa production in other coconut-growing regions by identifying the opportunity to intensify cocoa and coconut intercropping. Secondly, this work addresses the regional needs aspired to by the Bougainville government by identifying the land cover change in relation to land capability to meet the sustainability goals set in the Bougainville Strategic Plan 2018–2022. Demand for sustainably produced cocoa must be supported with technology, extension services, and consultation with governmental and non-governmental organisations. Further assessment on sustainability needs to be evaluated by defining a cocoa production extent threshold for sustainable development and the difficulty in its formulation or quantification. We need empirical and quantitative ways to assess it, therefore, this is a huge field for future work, along with region-specific research to identify sustainable production plans.