Improving Remote Monitoring of Carbon Stock in Tropical Forests Using Machine Learning: A Case Study in Indonesian Borneo

1. Introduction

1.1. Biomonitoring of Carbon in Tropical Forests

Tropical countries primarily contribute to carbon emissions through forest degradation and deforestation, and this accounts for approximately 10 percent of total global annual carbon emissions [1]. Local, national, and international initiatives (e.g., REDD+2) are underway to reduce carbon emissions. To assess whether these initiatives are delivering the expected results, it is crucial to develop precise, high-resolution, and cost-effective methods to track aboveground changes in carbon stocks.

Monitoring carbon stock in tropical rainforests can be expensive and challenging [2]. Traditional fieldwork to estimate aboveground carbon is a labor-intensive and time-consuming process, requiring at least 5-10 technicians working for weeks to measure a wide range of structural traits of trees (e.g., wood density of different tree species, density of individual trees, trunk diameters, canopy height) over a relatively small area (few hectares) to provide a localized estimate of carbon stock [2]. By replicating the same measures over a network of field plots, it is possible to derive an allometric relationship that estimates aboveground carbon density (ACD) based on these variables [3]. Whereas field methods provide accurate localized estimates, ACD extrapolation beyond the inventory plot networks remains challenging, as it would require deriving information on the same set of forest/tree structural traits at a fine scale over wide geographical areas [4]. Remote sensing methods are cheaper and are feasible at larger scales.

One useful remote sensing platform for monitoring and managing forest cover is Hansen Global Forest Watch (GFW) [5]. GFW data are available for free on the GFW Open Data Portal, which provides global coverage at several spatial and temporal resolutions. However, while useful for estimating forest area and changes on a yearly basis, Hansen GWF currently does not offer estimates of aboveground carbon density, nor can it be used for quantifying weekly or monthly changes in forest cover [5]. Therefore, novel datasets and pipelines are required in order to produce high-resolution, accurate time series data on ACD changes [6].

By providing high-definition, accurate three-dimensional data on forest canopy height and structure, airborne Light Detection and Ranging (LiDAR) optical sensing technology can be used to estimate ACD beyond the field plot inventory used to parametrize the allometric relationships between canopy height and ACD [6,7]. LiDAR-derived canopy height has been used to map biomass and estimate ACD in tropical forests including Central America [7], South America [8], Asia [9], Africa [10,11,12], and the oceanic islands [13]. However, high-definition airborne LiDAR data is also geographically limited in extent, a result of costs and logistical hindrances associated with the use of aircraft to carry the optical sensing technology [13,14]. Therefore, LiDAR is often paired with other data, such as optical remote sensing data of a range of spatial and spectral properties, to model LiDAR-derived properties across the larger geographic scale afforded by satellites [1,8,15,16].

Coupling airborne LiDAR along with satellite images and other geospatial datasets has become a more frequent approach to map tropical forest ACD, more so across geographic regions that do not have enough LiDAR measurements [17,18]. For instance, Asner et al. [19] coupled Landsat-derived parameters along with Shuttle Radar Topography Mission elevation metrics to regulate an ACD model that used airborne LiDAR for Peru. Baccini and Asner [20] used MODIS, Moderate Resolution Imaging Spectroradiometer, image data along with airborne LiDAR, to create pantropical high-resolution ACD maps. Notably, Mascaro et al. [6], for the first time, not only incorporated LiDAR data, but also used machine learning algorithms to estimate ACD. Today, evaluating forest properties, such as tree canopy heights (TCH) and ACD, from optical images can be accomplished rapidly via machine learning (ML) regression modeling [21,22,23,24,25].

1.4. Machine Learning (ML) to Estimate Forest Properties

Random Forest (RF) is the ML approach of choice, as it has proven to be superior to conventional ML techniques and tools for tropical forest carbon mapping applications [6]. This is in part due to RF being non-parametric and robust even with skewed data and a high number of input variables [26]. RF accuracy can be greatly enhanced by computing two widely used textural measures, namely Fourier transform textural ordination (FOTO) [27,28] and gray-level co-occurrence matrices (GLCM) [29,30,31]. Although textural features have been applied to imagery acquired with a variety of sensors, including WorldView-2 and -4 [32], Cartosat-1a [33], SPOT-5 [34], and IKONOS-2 [35], these texture features have yet to be tested on other gradient boosting decision tree models for large scale ACD mapping, such as XGBoost [36] and LightGBM [37], or with deep learning (DL) methods (e.g., neural network models).

In this study, we present an open-source, rapid, automated, fine-scaled, cost-efficient, and generalizable framework to accurately estimate ACD. By stacking multiple sources of coarse-resolution satellite imagery and remotely-sensed earth data atop the high-resolution canopy tree height LiDAR data, we have trained multiple ML algorithms, including tree-based models (e.g., RF, XGBoost, LightGBM) and DL models (i.e., neural networks) to estimate TCH and ACD values over areas that lack high-quality LiDAR data. Our work demonstrates a solution for increased accessibility of carbon estimation of forests from freely-available and open-source remote sensing data.

2. Materials and Methods

2.1. Study Area and Data Pipeline

Borneo, a 750,000 km² island, contains some of the world’s most biodiverse and carbon-dense tropical forests. However, Borneo has lost about 30% of its forests within the last 40 years [38]. In 2014, the NASA Carbon Monitoring System (CMS) surveyed 85 sites in Indonesian Borneo (Kalimantan) using plane-based LiDAR [39], creating canopy height digital surface models for 100,000 hectares of tropical forest TCH (Figure 1). We included this TCH data into our data pipeline and utilized it to calculate ACD at the pixel-level, as our ML regression target (see Section 2.2).

Like previous studies [6,40], we integrated several sources of common remote sensing data in our data pipeline, including Landsat-8 (8 spectral bands), Sentinel-1 (4 band Synthetic Aperture Radar), land cover fractional coverages from the Copernicus satellite mission, digital elevation model, and NASA CMS airborne LiDAR. Because of data loss due to cloud cover and other atmospheric effects, Landsat data was temporally sampled from August 2014 – January 2015 (the time period of the LiDAR data plus 2 months before and after). In addition to these data, we also integrated other data sources, such as indices commonly used to assess vegetation (e.g., NDVI, NDWI, NDII, EVI, for Landsat 8), Gray-Level Co-occurrence Matrix (GLCM) texture metrics [29] for all spectral data, and climatic data from WorldClim [41]. All these data sources (listed in Table 1) were compiled into a data processing pipeline which produced our dataset of 336 features for each LiDAR value (examples in Figure 2), resulting in about 2 million observations. The definition of each data feature is listed in Supplementary Material Table S1.

All data was reprojected to the coordinate reference system of the corresponding Landsat imagery (30 × 30 m pixel); LiDAR data was resampled using bicubic interpolation (Figure 3). Land Cover classification values were transformed via label encoding. This diverse remote sensing dataset was critical to providing the foundation on which we produced a robust, predictive ML framework. Here we utilized Google Earth Engine (GEE) [42] programmatically for data acquisition, processing, analysis, and modeling. The data processing pipelines we built in GEE can be further used to facilitate new data acquisition and processing for any future work in different forest regions and scope expansion.

2.2. ACD Estimation

In this study, we used the empirical relationship specific to Borneo, built by [9], to derive ACD based on TCH from airborne LiDAR. The empirical equations are the follows:

ACD = 0.567 × TCH^0.554 × BA^1.081 × WD^0.186

BA = 1.112 × TCH

WD = 0.385 × TCH^0.097

where ACD is in Mg C ha^-1, TCH is in m, basal area (BA) is in m² ha^-1, and wood density (WD) is in g cm^-3. Jucker et al. [9] combined field data of tree height measurements with airborne laser scanning and hyperspectral imaging to characterize how topography shapes the vertical structure, wood density, diversity and ACD of nearly 15 km² of old-growth forest in Borneo. This general form of empirical relation was first developed by Asner et al. [8] and the allometric scaling coefficients in the equations are region-specific and are based on the field work that took place in that region. The error range of the ACD values calculated from the empirical equations is 39.3 Mg C ha^-1 [9].

Figure 3. (Top) An original full extent, full resolution (1 × 1 m) LiDAR-derived Canopy Height Model (CHM) with a rectangle (pink) showing inset. (Left) Zoomed-in inset of CHM at original 1 × 1 m resolution, showing individual trees as light-colored circular shapes and a river/low ground in dark. (Right) CHM image after resampling to Landsat-8 scale (30 × 30 m).

Figure 3. (Top) An original full extent, full resolution (1 × 1 m) LiDAR-derived Canopy Height Model (CHM) with a rectangle (pink) showing inset. (Left) Zoomed-in inset of CHM at original 1 × 1 m resolution, showing individual trees as light-colored circular shapes and a river/low ground in dark. (Right) CHM image after resampling to Landsat-8 scale (30 × 30 m).

2.3. ML Framework

In our framework of ACD estimation from ML modeling, we experimented with several ML and DL algorithms, including artificial neural network [43], random forest (RF) [44], XGBoost [36], and LightGBM [37], as well as the most recent AutoMLs [45,46].

From our data pipeline, we ran each ML algorithm with our dataset containing 336 features and 2 million TCH and ACD observations. The regression target for each ML model is the ACD value in pixel level. The goal was to develop a generalized ML model that is capable of predicting ACD solely based on open-source remote sensing data in the regions that LiDAR TCH is generally not available. The framework workflow is presented in Figure 4.

The workflow for ML model development consisted of two stages: (1) training and validation and (2) testing. The selected ML model takes the input data (GEE compiled remote sensing data) and ACD as the predicted output values. During the training and validation phase, the cost (or loss) function, the mean squared error (MSE), was used to evaluate the performance of the ML model by comparing the ACD estimation to the output values predicted by the ML model. Since this is a regression task, the evaluation metrics we used here is Coefficient of determination (R²). We split the dataset into a 90% training subset, a 5% validation subset, and a 5% test set. The subsets were kept consistent for each round of experiment of the selected ML algorithms. The validation set was used to evaluate if the model is over-fitting to the training data, while the test set was treated as out-of-bag data that was not used in in the training phase.

Note that in general, with large datasets (e.g., > 1 million data points), the artificial neural network (ANN) model is more generalizable and likely has higher predictive performance than the tree-based algorithms, such as RF. As RF has been tested rigorously in the remote sensing literatures in recent years, we treat the tree-based algorithms tested in this study as the baseline.

We used Python environment to establish our code base and ML framework. The ML computation, model training, and inference were done in Google Cloud Platform’s Vertex AI.

Figure 4. ML framework with GEE remote sensing data pipeline and GCP infrastructure with ML model development and deployment for ACD prediction.

Figure 4. ML framework with GEE remote sensing data pipeline and GCP infrastructure with ML model development and deployment for ACD prediction.

2.3.1. Tree-Based Models

In addition to commonly used RF in the remote sensing community, we experimented with gradient boosted tree algorithms, such as XGBoost and LightGBM, as they have previously been effective in remote sensing image analysis and classifications of forests and land cover [47,48,49].

RF is a supervised ensemble algorithm that generates unpruned classification trees to predict a response [44]. RF implements bootstrapping samples of the data and randomized subsets of predictors. Final predictions are ensembled across a forest of classification trees (n = 100) based on an averaging of the probabilistic prediction of each classifier. Tree depth, number of leaf nodes, and number of features are tested with the dataset in each iteration in the training process.

XGBoost (eXtreme Gradient Boosting) is designed for both linear and tree-based supervised models [36]. XGBoost builds an ensemble of weak learners, where each weak learner is a decision tree. These decision trees are trained in a sequential manner, with each tree trying to correct the mistakes of the previous tree. The final output is a combination of all the individual trees, which results in a highly accurate model. XGBoost also includes a variety of regularization techniques, such as L1 and L2 regularization, and the ability to handle missing values and categorical variables.

Like XGBoost, LightGBM also builds an ensemble of weak learners, where each weak learner is a decision tree [37]. The key difference between LightGBM and XGBoost is the way they handle the decision tree learning process. LightGBM uses a technique called gradient-based one-side sampling to sample the data, which improves the speed of the training process. It also uses an algorithm called exclusive feature bundling to handle categorical variables, which reduces the number of splits required in the decision tree. Another difference between LightGBM and XGBoost is that LightGBM uses a leaf-wise tree growth algorithm, whereas XGBoost uses a level-wise tree growth algorithm. The leaf-wise algorithm tends to grow deeper trees than the level-wise algorithm, which results in better accuracy but can also lead to overfitting if not carefully controlled.

2.3.2. Neural Network Models

Deep learning algorithms (i.e., neural network models) have previously been effective for discovering underlying patterns in large datasets [43]. Therefore, we expected that neural network models could find highly predictive features from our diverse remote sensing data, which contained many data features (336) and data points (2 million). We specifically constructed the simple form of feedforward artificial neural network (ANN).

Feedforward ANNs consist of inter-connected computational units referred to as neurons that are arranged in layers. Input data is passed in one direction through an input layer and subsequently propagated through a defined number of hidden layers whereby the sum of the products of the connection weights from each layer approximate a function to estimate the observed output values [43]. Under repeated iteration and adjustment of the connection weights, a function between inputs and output is as closely estimated as possible given the parameter space available in the network [50,51]. Detailed development of model architecture lies in the selection of node functions (i.e., activation functions), layer sizes (number of hidden layers and numbers of nodes in each layer), learning rate, regularization parameters, and parameter dropout, which is discussed in Section 2.3.3.

As described above, the neural network used the input parameters as the first layer of neurons, and the last layer of neurons represents the predicted output values. Stochastic gradient descent, a common optimization method for ANN, was then used to iteratively adjust the weights and biases for each neuron to allow the ANN to best approximate the training data output. At each iteration, a backpropagation algorithm estimated the partial derivatives of the cost function (MSE) with respect to incremental changes of all the weights and biases, to determine the gradient descent directions for the next iteration. Note that in our model, the neurons of each hidden layer were composed of Rectified Linear Units (i.e., a ReLU activation function), to avoid the vanishing gradients [52]. Validation data were not used in the optimization or backpropagation algorithms. Instead, the cost function was evaluated over the validation data and served as an independent tuning metric of the performance of the ANN.

2.3.3. AutoML (Automated Machine Learning)

AutoML (Automated Machine Learning) is a technique that automates the process of building machine learning models, from data preparation to model selection and hyperparameter tuning [45]. The goal of AutoML is to make the process of building models more accessible and efficient for ML practitioners [46]. In this study, we primarily utilized AutoML for model selection and fine-tuning, as manually constructing the architecture of ANN model and fine tune hyperparameters, such as number of hidden layers, numbers of nodes in each layer, learning rate, and dropout rate, can be a time-consuming and complex process. We implemented a Python library, Auto-Keras [53], to automate the search process for ANN architecture that yields the best model prediction R². To evaluate the robustness of the ANN model, we compared it to the best performing tree-based model as our baseline. We implemented a Python library, TPOT, a tree-based pipeline optimization tool for automating machine learning [54], to automate the search process for tree-based models that yields the best model prediction R².

3. Results

In our experiments, with 336 data features and 2 million data points on ACD compiled from our GEE processing, the tree-based ML models achieved good R² performance on the test set. Our tree-based baseline ML models all have similar and consistent predictive performance: R² = 0.572 for the RF model, 0.583 for the XGBoost model, and 0.581 for the LightGBM model. The 3 most important features of each model are presented in Table 2. TPOT, the AutoML approach, found the best tree-based model to be the decision tree with fine-tuned parameters (maximum depth = 9, minimum sample leaf = 2, minimum samples split = 15) and the best R² was 0.593. These results are better than the comparable studies [6,55], here used as benchmark (R² of 0.3-0.5), due to the large number of features used in our study. For instance, the RF model in Mascaro et al. [6] was trained with data that contained only multispectral satellite channels and land cover features, which had a total of 11 data features and 80,000 data points.

For the DL experiments, the AutoML approach via Auto-Keras found the best neural network model to be 3 dense layers that applied multi-category encoding and normalization; each hidden layer had 32 neurons. ReLU activation function layers are implemented twice in this best-fit model, which allows for the neural network to learn hierarchical non-linear representations of the input data. The model architecture summary is presented in Table 3. This neural network had 13,090 total parameters in which 11,296 parameters were trainable. This custom feed-forward ANN achieved R² = 0.821 on the test set, which is higher than our tree-based baseline models and outperformed comparable studies [6,55]. After obtaining the best neural network model, we used it to produce a prediction ACD map over a section of Borneo where high-resolution LiDAR does not have coverage (Figure 5).

4. Discussion

Here we show that large remote sensing datasets produce tree-based ML models with high predictive accuracy. Our three tree-based ML models trained with 336 data features and 2 million data points all had R² values greater than 0.57, which is substantially higher than prior studies with fewer features and data points that had R² values in the 0.3 – 0.5 range [6,55]. This resonates with the recent data-centric AI approach, which emphasizes the importance of improving data quality and quantity, instead of focusing on improving ML algorithms [56,57,58].

We also found that our deep learning neural network model yielded substantially higher R² than the tree-based models (R² > 0.8). This is consistent with previous studies that showed that neural networks are more generalizable and capable of extracting underlying patterns from a large dataset, resulting in a higher predictive performance [59]. This is a large benefit, but there is also a cost, because neural network models are more difficult to interpret than tree-based models. While it is mathematically possible to determine which nodes of a deep neural network were activated, it is difficult to assess how each of the neuron layers of neurons performed collectively. This is less straightforward than interpreting conventional tree-based models, which provide crisp rules on feature selections and encode for feature importance, making it easier to interpret the reasoning behind the model. In this study, we focus more on model predictive performance than interpretability, so we considered the DL based neural network modeling to be the best approach for our ADC estimation application.

As our models were not built for interpretability, but rather to maximize predictive power by including a large number of features, it is not sensical to select a single feature and comment on its significance. However, we did calculate the feature importance from the RF, XGBoost, and LightGBM models. In general, the top features that have the most predictive power for the ACD estimation were Sentinel-1 VH Cluster Shade (a GLCM texture feature), Landsat-8 thermal and infrared GLCM features, and Landsat-8 Red and Blue GLCM features. GLCM features comprised the top 5 features of every model run, which indicates the significance of texture analysis of medium-scale satellite imagery in the analysis of forest canopy structure.

Note that the Landsat imagery used in our data pipeline is 30 × 30 m, so the sub-meter LiDAR data was resampled to derive mean height in each 30 × 30 m plot. The new generation of high temporal and spatial resolution satellite data, e.g., PlanetScope, which began observations over Borneo in 2018, consists of pixel cell sizes at ~3-5 m² and would thus reduce the amount of data resampling, potentially providing a better estimation of TCH and, in turn, of ACD [60]. PlanetScope sensors capture eight spectral bands, the same number as from Landsat, but don’t have a thermal infrared sensor. Nevertheless, in our models, texture features derived from red, blue, and infrared spectra consistently were ranked the most important, suggesting that PlanetScope sensors would still capture the most useful information, especially as the amount of textural information would increase greatly with the finer resolution. The use of PlanetScope imagery is particularly appealing as a future avenue of research because PlanetScope satellites now cover the entire globe daily, allowing high temporal frequency imagery of the area of interest in Borneo. This would greatly help with dealing with cloud cover, which is a common challenge in rainforests like the ones we considered in Borneo.

5. Conclusions

In this study, we developed an improved method for evaluating fine-scale carbon stock by integrating LiDAR with satellite data and using a range of ML and DL algorithms. We found that Neural Network models trained on a large dataset of observations and features outperformed tree-based ML algorithms and allowed to achieve better prediction accuracy than any previous studies. Future efforts will employ Google Cloud to host and deploy this trained model on GEE. Our pipeline can be modified to estimate ACD in regions other than Borneo where location-specific equations exist to relate TCH to ACD and will be especially accurate in regions lacking large LiDAR coverage. This will provide insights for the necessary actions required to mitigate the effects of climate change.