Assessing the Spectral Information of Sentinel-1 and Sentinel-2 Satellites for Above-Ground Biomass Retrieval of a Tropical Forest

Abstract

Earth Observation (EO) spectral indices have been an important tool for quantifying and monitoring forest biomass. Nevertheless, the selection of the bands and their combination is often realized based on preceding studies or generic assumptions. The current study investigates the relationship between satellite spectral information and the Above Ground Biomass (AGB) of a major private forest on the island of Java, Indonesia. Biomass-related traits from a total of 1517 trees were sampled in situ and their AGB were estimated from species-specific allometric models. In parallel, the exhaustive band combinations of the Ratio Spectral Index (RSI) were derived from near-concurrently acquired Sentinel-1 and Sentinel-2 images. By applying scenarios based on the entire dataset, the prevalence and monodominance of acacia, mahogany, and teak tree species were investigated. The best-performing index for the entire dataset yielded R² = 0.70 (R² = 0.78 when considering only monodominant plots). An application of eight traditional vegetation indices provided, at best, R² = 0.65 for EVI, which is considerably lower compared to the RSI best combination. We suggest that an investigation of the complete band combinations as a proxy of retrieving biophysical parameters may provide more accurate results than the blind application of popular spectral indices and that this would take advantage of the amplified information obtained from modern satellite systems.

1. Introduction

1.1. Tropical Forest Biomass and Importance of Mapping

Approximately 31% of the Earth’s surface is covered by forests [1], which provide essential ecosystem services and host enormous biodiversity. Forests hold an estimated 861 ± 66 Pg of carbon stock, of which 42% resides in biomass [2], a fact which is particularly noticeable in tropical forests [3]. Therefore, tropical forests play a paramount role in global carbon sequestration and climate change; this has been a major reason for countries and organizations to report periodically on their forest status through schemes such as the Reducing Emissions from Deforestation and Forest Degradation (REDD).

Despite the fact that the rate of forest loss has decreased in recent years, it continues to take place, with an estimated 3% reduction in global forested area from 1990 to 2015 [4], mainly due to clearings undertaken to prepare for agriculture expansion. Therefore, monitoring human-induced fluctuation in forest biomass is of paramount importance. Nevertheless, quantifying Above Ground Biomass (AGB) or below ground biomass in the tropics is especially challenging due to the large AGB stock, arduous accessibility to dense forest structures and the associated lack of inventory data. It is noteworthy that AGB is indirectly related to below-ground biomass, and species-specific models have been developed in this framework based on related field campaigns. For example, Singnar et al. [5] showcased the relationship between these two traits in a bamboo forest in northeast India; Liu et al. [6] demonstrated an accurate estimation of the above-ground to below-ground biomass ratio in L. leucocephala; and He et al. [7] derived a coefficient of determination, (R²) = 0.89, and a corresponding significance level, (p) < 0.001, between below- and above-ground biomass cumulatively for 12 species in a natural coniferous and broadleaf mixed forest in Northeastern China.

The quantification of forest biomass is traditionally conducted with destructive sampling, which involves harvesting tree samples, requiring the felling, dissecting, oven-drying, and weighting of trees. While several studies have been based on this methodology (e.g., [8,9,10]), it is unavoidably a resource- and labor-intensive procedure [11], especially in tropical areas where carbon accumulation occurs in the biomass instead of the soil [2]. Tree growth, on the other hand, is traditionally assessed either by repeated measurements of marked trees in permanent sample plots or tree ring analysis (dendrochronology). The latter technique is applicable primarily in temperate environments and species, as it requires tree samples with well-defined annual growth rings [12]. Alternatively, allometric equations have been suggested to estimate AGB (e.g., [7,13,14,15]) and related to height, diameter at breast height (DBH), species, age, crown density, and bioclimatic variables [16]. However, allometric equations are dependent on species composition and landscape characteristics; nonetheless, a few studies have proposed robust equations, such as that of Chave et al. [17], which introduced a pan-tropical model incorporating wood-specific gravity, trunk diameter, and tree height. Consequently, interest in estimating AGB on large spatial scales has brought remote sensing to the forefront of current research, and numerous attempts have been made to approximate AGB based on Earth Observation (EO) data.

1.2. Earth Observation and the Use of Sentinels

EO has played a vital role in biophysical parameter retrieval, especially in recent years, during which the abundance of satellite data and the tools accessible to ecologists has surged [18], with the former becoming increasingly available freely and in near-real time. In the context of biomass retrieval, several sensor types have been used such as optical (e.g., [19,20,21]), Synthetic Aperture Radar (SAR) (e.g., [22,23,24,25,26]), Light Detection And Ranging (LiDAR) (e.g., [27,28,29]), and a fusion of the above (e.g., [30,31,32,33,34,35,36,37]). Other experimental datasets have also been used such as Google Earth images. On this basis, Ploton et al. [38] claimed to have undertaken reliable mapping of tropical forest AGB in the Western Ghats of India. In tropical areas, studies have concentrated on the valorization of SAR imagery as a means to circumvent the frequent cloud presence encountered in this climate (e.g., [23,39,40,41]).

Optical satellites have long been used to retrieve biophysical parameters such as chlorophyll content [42] and AGB, with a wide range of reporting accuracies. For instance, in a study estimating grass AGB across different fertilizer treatments, Sibanda et al. [43] reported correlations of R² = 0.81 and R² = 0.76 for Sentinel-2 and Landsat 8 imagery, respectively. Pandit et al. [19] reported a correlation of R² = 0.81 in an estimation of the AGB of a subtropical Nepalese forest based on Sentinel-2 imagery. In contrast, Majasalmi and Rautiainen [44] reported a low R² = 0.24 for AGB in a study on three boreal forests in Finland, while for other biophysical variables, namely Fraction of Absorbed Photosynthetically Active Radiation (fPAR) and Leaf Area Index (LAI), the correlation was considerably higher. More sophisticated methods have examined the use of Artificial Neural Networks (ANN) and multiple input layers such as day of year, canopy height, and NDVI [45].

Nevertheless, the usability of optical images is hindered by two main constraints. First, cloud coverage can impede regular image acquisition and, hence, reduce the temporal availability of satellite observations. Secondly, certain biophysical parameters, including biomass, saturate at low values of remote sensing data [46]. As an example, LAI can be retrieved based on a linear regression of an optimized vegetation index (e.g., [47]), the accuracy of which relates to the vegetation type. For instance, Chen et al. [48], in a study conducted in a boreal environment, showcase that, for conifer stands, there does not seem to be a saturation point (up to LAI = 10, which is the highest reported in this case study); however, in a deciduous forest, LAI saturates at values of 2–3, which is similar to the case of crops. This saturation effect becomes even more prominent in tropical environments [49] and restricts the usability of optical satellite images. The retrieval of biophysical parameters also depends on factors such as the soil type and plant architecture; hence, information about the landscape is crucial before the selection of the spectral information for the estimation of biophysical variables [50]. SAR sensors, on the other hand, are not considerably hampered by clouds and are independent from solar angle and radiation intensity. However, saturation of the SAR signal often exists in both the broadleaf and conifer tree-types [51]. Lastly, a few studies have attempted to use synergistically optical and SAR data for biomass estimation (e.g., [35,52]), and they report improvement of accuracy when fusing datasets of differing natures. A review of remote sensing for AGB can be found in Tsitsi [53].

1.3. Biomass Retrieval—Optimal Spectral Bands and Indices Selection

Traditionally, biomass is estimated based on allometric equations or empirical relationships between the satellite data and the biomass measured in situ. However, these two approaches are species- and location-specific and therefore cannot be extrapolated to other species and locations. More complex approaches have been introduced such as Artificial Neural Networks [54], physically-based reflectance models (e.g., PARAS model by [44]), Support Vector Machines [55], Support Vector Regression [35], Random Forests [19,37], and Sparse Partial Least Squares Regression [21]. Nevertheless, it is important to note that Hill et al. [56] claimed that methodological complexity does not lead to improved precision; as a matter of fact, they suggested that a way to improve precision is by integrating uncertainty estimates.

1.4. Band Selection and Complete Band Combinations

Most often, the retrieval of a biophysical variable based on remotely sensed data is derived from a preferred individual spectral wavelength range, such as the red-edge band which highlights vegetation vigor. Alternatively, a spectral index is employed as a proxy for a vegetation trait [57]; for instance, the Normalized Difference Vegetation Index (NDVI) has been a popular and frequently used spectral index. However, several studies have indicated that this is not always the solution yielding the most accurate results for certain vegetation traits [58,59]; it is often the case that information inherited in bands is overlooked, leading to the glorification of mainstream spectral indices, such as the NDVI in the case of AGB quantification.

The current study scrutinizes the selection of the spectral information inherited in the optical and SAR satellite systems for the retrieval of AGB. We calculate the exhaustive spectral combinations from the Sentinel-1 and Sentinel-2 satellites to investigate the degree to which the full spectral information is utilized when retrieving AGB. Based on in situ ground truth information collected in a private forest in Indonesia, we compare the performance of the satellite-based retrieval approach proposed and benchmark our findings against the AGB estimation based on traditional vegetation indices. Finally, we showcase the superiority of examining the full potential of the exhaustive band combinations and suggest it as a methodology for remote sensing-based biophysical retrieval.

2. Materials and Methods

2.1. Study Area

The study area is located in the island of Java, within the Girisekar and Jetis private Forest Management Units (FMUs) in Gunung Kidul Region, Yogyakarta Province, Indonesia. The Gunung Kidul Region covers an area of 3133 km², has a highest elevation of 2955 m, and receives precipitation of approximately 2681 mm. The two adjacent FMUs cover a cumulative area of 2666 ha and are situated at the geographical coordinates 8°02′20.4″ S and 110°28′42.6″ E (Figure 1). The greater Gunung Kidul Region has become a major area of private forest activity and has attracted interest surrounding the forestry’s community-based certification, management, and contribution to economic development (e.g., [60,61,62,63,64]). Private forest plantations in Gunung Kidul are comprised predominantly of teak (Tectona grandis), although efforts have been undertaken in planting acacia (Acacia mangium), mahogany (Swietenia mahagoni), munggur (Sasmania saman), and sono (Dalbergia latifolia) [61]. Deforestation in the island of Java dates back to 300–400 years ago [65] and began to be perceived as a problem from the middle of the 19th century onwards [66]. Gunung Kidul Region specifically has observed a transition to re-vegetation since the 1970s [67].

At the country level, Indonesia has suffered intense deforestation and degradation during the last decades primarily due to the conversion of forest into estate crops and agriculture and logging concessions [65,68]. Despite the fact that the Global Forest Resources Assessment [1] reported reduced deforestation between 1990 and 2015 at a global scale with a few tropical countries transitioning to forest gain, specifically for the tropical counties, the forest loss and degradation persists [69].

Indonesia has been a main contributor of deforestation globally. Stibig et al. [70] estimated that, during the period 2000–2010, 0.82 Mha of forest per year were lost, accounting for 56% of the total forest cover loss in Southeast Asia. The consequent greenhouse gas emissions from deforestation, and the concomitant vulnerability of Indonesia to climatic change, sensitized the Indonesian government to recognize the need for action on climate change in 2009 [71]. Subsequently, a reduction in greenhouse gas emissions by 29% by 2030 compared to the business as usual (BAU) scenario was pledged by taking actions such as sustainable forest management and restoring degraded ecosystems [72]. Mitigation efforts included the REDD project, peat land conversion, and forest management and governance [73]. One specific way of supporting ecosystem and landscape resilience has been through social forestry [72,74] and community forest management during the last decade in Indonesia, which have fostered important efforts [75].

2.2. Datasets

Thirty plots were selected based on a stratified random sampling and sampled in September 2017. Each plot was defined as a rectangle sized 20 m × 20 m, which is the nominal pixel size for most bands of Sentinel-1 and Sentinel-2. The latitude and longitude coordinates of the center of each plot were recorded with a Global Navigation Satellite System (GNSS) device. Trees with a DBH of more than 8 cm were considered and the species, circumference of the tree trunk at breast height (1.3 m above the ground) and height of each specimen were recorded. A total of 1003 trees were sampled (

Table 1. Tree species recorded during the fieldwork, and the species-specific allometric equation used to estimate the biomass (adapted from Purwanto et al. [15]).

Common Name	Scientific Name	Allometric Equation
Teak	Tectona grandis	AGB = 0.0149x(D2xH)1.0835
Earleaf acacia	Acacia auriculiformis	AGB = 0.0775x(D2xH)0.9018
Mahogany	Swietenia mahagoni	AGB = 0.9029x(D2xH)0.6840
Others (gembilene, melinjo, pulai, durian, coconut, kayu besi, trembesi)	-	AGB = 0.0219x(D2xH)1.0102

where AGB: total biomass of stems, branches, and leaves; D: diameter of the tree trunk at breast height calculated from the circumference measurement; and H: height of the tree.

). Plots containing more than 90% of a single tree species from all the trees recorded were considered as monodominant. It is worth noting that frequently monodominance is frequently defined as ≥60% of canopy-level trees belonging to the same species with values regularly attaining 80–100% [76]. Nevertheless, in the current study we set a higher threshold (i.e., 90%) as trees with a very small diameter were excluded from field sampling.

The remote sensing dataset used in the current study is comprised of the Sentinel-1 and Sentinel-2 satellites, two twin constellation systems from the Copernicus programme [77]. Sentinel-2 is an optical satellite with 13 spectral bands and 290 km swath width supporting primarily land use and land cover applications [78]. Sentinel-1 is the corresponding C-band SAR satellite mission providing Vertical transmit-Horizontal (VH) and Vertical transmit-Vertical (VV) polarizations. The Sentinel-1 and Sentinel-2 satellite constellations have a revisit frequency at the Equator of 3 and 7 days, respectively, and cover the majority of the Earth, making them unprecedented information sources for frequent monitoring of biophysical parameters at a global scale. Satellite data were selected on the basis of the temporal proximity of the image acquisition and the field data and taking into account the quality assessment of the cloud coverage for the optical data over the area of interest.

Sentinel-2 includes a ShortWave InfraRed (SWIR) band at 1.375 μm dedicated to Cirrus clouds and a respective cloud mask is provided with each delivery. However, cloud coverage is spatially inconsistent in tropical areas; for this reason, the overall cloud coverage statistics provided were disregarded during the selection of the appropriate image. Instead, the visual assessment of the cloud coverage over the area of interest was preferred. A cloud-free Sentinel-2B image acquired on 12 August 2017 was selected and downloaded from the Copernicus Open Access Hub. The image was provided at Top Of Atmosphere (TOA) 1-C level. In addition, a Sentinel-1 image acquired on 4 September 2017 was used; the standard high resolution Ground Range Detected (GRD) product selected was sensed in Interferometric Wide swath (IW) acquisition mode in dual VV and VH polarization. A suitable Sentinel-1 image with dual HH and HV polarization for the area of interest was not found in the Copernicus Open Access Hub. It is worth noting that while the Sentinel-1 image fell within the field data collection time window, the Sentinel-2 image was acquired approximately 3 weeks earlier as a consequence of a worse candidate due to local cloud coverage conditions. Nevertheless, considerable AGB growth is not expectable within 23 days, and, moreover, the climatic conditions remained relatively stable during this time period; the climate is characterized as a tropical monsoon climate (Am in the Köppen–Geiger system), with a short dry season from July to September [79] when the rainfall remains relatively low.

2.3. Methods

The AGB was estimated from allometric equations based on the tree height and circumference measured in the field. A collection of allometric equations (Table 1) was used to estimate AGB for each individual tree based on D and H as the dependent variables. The allometric equations were adapted from Purwanto et al. [15], who studied the biomass distribution of 25 species of perennial woods, including mahagony, acacia, sonokeling, teak, sengon, and others, in a community forest in Nglanggeran Village, Indonesia. In the current study, a separate allometric equation was selected for each species as it has been proven that biomass predicting power is higher for species-specific models than for all species combined [7]. Thereafter, the AGB estimates for each individual tree were added to extract the representative AGB of each plot (in ton/ha).

The Sentinel-2 image was first atmospherically corrected with the Sen2Cor plugin version 2.8 [80] to construct the level 2A surface reflectance product. Sentinel-2 bands inherit various spatial resolutions (at 10, 20, and 60 m), and the spatial resolution of each band was kept unchanged during the atmospheric correction. The bands 1, 9, and 10, acquired at 60 m spatial resolution, were dismissed from further processing as they relate mainly to cloud screening and atmospheric correction. Finally, the bands acquired at 20 m spatial resolution were resampled to 10 m with the nearest neighbor algorithm to retain the native spectral information of the image.

In regard to the Sentinel-1 imagery pre-processing, several studies have presented efficient processing workflows for a variety of applications (e.g., [81,82,83,84]). In the current paper, we follow a baseline pre-processing of sequential steps undertaken with the open-source software Sentinel Application Platform (SNAP) version 8.0.0. First, the orbit file was applied based on 3 polynomial degree followed by the thermal noise removal. Then, the sigma nought (σ⁰) was calculated and the Lee Sigma speckle filter with a window size 7 × 7 and sigma = 0.9 was applied on both polarizations to eliminate noise. Thereafter, a terrain correction operation was implemented using the Shuttle Radar Topography Mission (SRTM) 3 sec Digital Elevation Model (DEM) and bilinear interpolation at 10 m spatial resolution, and, finally, the pixel values were converted to dB.

The two Sentinel images were then subset to the Girisekar and Jetis FMU extent (Figure 1). The nearest neighbor resampling algorithm was used to preserve the raw pixel values. Thereafter, all Sentinel-1 and Sentinel-2 bands were stacked and reprojected to the WGS84 coordinate system. Correspondingly, the GNSS coordinates recorded at the center of each field plot were converted from Degrees-Minutes-Seconds (DMS) to decimal degrees. Subsequently, the Sentinel-1 and Sentinel-2 band values of the pixels containing the aforementioned coordinates were extracted. The correlation between the spectral indices and the spatially corresponding in situ-collected data was examined through the exhaustive combination of the Ratio Spectral Index (RSI) and the Normalised Difference Spectral Index (NDSI) spectral indices, which are defined in Equations (1) and (2).

RSI_i,j = R_i/R_j

(1)

NDSI_i,j = (R_i − R_j)/(R_i + R_j)

(2)

where R_i and R_j are the reflectance values of bands i and j, respectively.

The coefficient of determination (R²) and the p-value (p) at the critical thresholds (i.e., 0.050, 0.010, and 0.001) were used as the metrics to evaluate the performance of the statistical relationship. Cumulative R² and p graphs were constructed to aid the visualization and comparison between each pair of bands as in similar studies [59,85,86,87]. Furthermore, several scenarios were investigated based on the structural components of the private forest; the biomass relationship for the main individual tree species sub-group (i.e., acacia, mahogany and teak) was investigated separately, as well as the case of monodominant and non-monodominant forest stands. Subsequently, the indices of the spectral pairs with the highest correlation with in situ measurements were identified and the linear regression was plotted in an attempt to investigate the relationship of the most robust spectral combination resulting from the methodology of the current study.

In order to quantitatively evaluate the importance of our results, 8 popular vegetation indices frequently encountered in the literature were extracted from the same satellite dataset and the corresponding R² and p were calculated. Several popular vegetation indices are a specific case of RSI or NDSI (e.g., NDVI [88], Green NDVI [89]) and, hence, are already included in the complete combination of spectral indices of this experiment. Therefore, we primarily considered indices with a more mathematically complex definition to avoid duplication in the results. The list of selected indices is provided in

Table 2. Popular vegetation indices selected for the current study, their respective equation, and bibliographic references. Bands wavelengths used in the equations are translated to Sentinel-2 central bands according to the following correspondence: B (Blue): band 2 at 490 nm, R (Red): band 4 at 665 nm, NIR (Near InfraRed): band 8 at 842 nm, REDGE (Red Edge): band 8a 2 865 nm, REDGE2 (Red Edge 2): band 5 at 705 nm, and REDGE3 (Red Edge 3): band 6 at 740 nm. L is the canopy background adjustment factor = 0.5.

Index	Reference	Equation
EVI (Enhanced Vegetation Index)	[92]	$\frac{2.5\ast R_{NIR}-R_R}{R_{NIR}+6\ast R_R-7.5\ast R_B+1}$
SAVI (Soil Adjusted Vegetation Index)	[93]	$\frac{\left(1+L\right)\left(R_{NIR}-R_R\right)}{\left(R_{NIR}+R_R+L\right)}$
MSAVI2 (Improved Soil Adjusted Vegetation Index)	[94]	$\frac{2R_{NIR}+1-\sqrt{{\left(2R_{NIR}+1\right)}^2-8\left(R_{NIR}-R_R\right)}}{2}$
ARVI (Atmospherically Resistant Vegetation Index)	[95]	$\frac{R_{REDGE}-2\ast R_R+R_B}{R_{REDGE}+2\ast R_R-R_B}$
mNDVI (modified Normalized Difference Vegetation Index)	[96]	$\frac{R_{REDGE3}-R_{REDGE2}}{R_{REDGE3}+R_{REDGE2}-2\ast R_B}$
CRI (Carotenoid Reflectance Index)	[97]	$\frac{1}{R_B}-\frac{1}{R_R}$
NDVI (Normalized Difference Vegetation Index)	[88]	$\frac{R_{NIR}-R_R}{R_{NIR}+R_R}$
VOG1 (Vogelmann)	[98]	$\frac{R_{REDGE2}}{R_{REDGE1}}$

. The results of the popular vegetation indices were then judged against the complete combination rations. The Geospatial Data Abstraction Library (GDAL) [90] and the R programming language [91] were used for the analysis.

Finally, to validate the algorithm and prove the robustness of the proposed approach, we applied the same methodology with an extended field dataset collected between September and November 2017 consisting of 45 plots (a total of 1517 trees). While this dataset includes an increased number of observations, a longer period is considered extending beyond the date of the satellite image acquisition; thus, the concurrency of field data and satellite imagery is not consistent as in the development phase of the current study. However, this experimental set up allowed us to reproduce the results and investigate if the R² graph resembles the one produced in the controlled experiment and whether the proposed band combinations are still superior to the traditional vegetation indices. The analysis was rerun with the same satellite imagery.

3. Results and Discussion

The range of the structural characteristics of the 1003 trees that were recorded is provided in

Table 3. Cumulative number and species of trees sampled in the 30 plots and their structural characteristics.

Species	Number of Trees	Circumference Range (cm)	Height Range (m)	Biomass Range (kg)
Teak	550	20–98	1.23–23.98	4.63–486.27
Earleaf acacia	165	19–100	4–60–25.66	9.00–699.18
Mahogany	268	25–94	4.00–23.98	41.98–829.65
Others	20	27–100	5–19	12.08–129.41

. The majority of the samples are teak (54%), while mahogany and earleaf acacia account for the rest (43%). The range of the circumference and height for all samples is 19–100 cm and 1.23 m–25.66 m, respectively, while intra-variability for these two attributes seem to be similar across all species. However, variability in biomass is noticeable especially at the upper limit of the range: teak caps at 486 kg, earleaf acacia at 699 kg, mahogany at 829 kg, and other species at 129 kg. This is a consequence of the tailored allometric equations used for each species and the different coefficients applied as presented in Table 1. Despite the fact that generic allometric equations have been developed and employed in the literature (e.g., [17]), it has been shown that biomass predicting power is higher for species-specific models than for all species combined [7].

Figure 2, Figure 3 and Figure 4 present the statistical distribution of the structural characteristics categorized by species (Figure 2), monodominance (Figure 3), and dominance and species (Figure 4). Figure 2 depicts the histogram distribution of the biomass, diameter, and height for the three major species. Regarding the biomass distribution, Askar et al. [20] estimated the average AGB in this private forested area to be 72.54 ton/ha; private forests contain a relatively lower AGB compared to natural tropical forests as it is the case that, in the former, mature trees are harvested, while the latter accumulate large amounts of terrestrial AGB over long periods of time [99]. From the findings of the current study, it is apparent that the diameter distributions of the three tree species (colored curves) have similar characteristics and shape; consequently, the distributions of biomass and height are similar. This conclusion leads us to the observation that, for the current study area and species, the diameter does not associate with biomass, and the height of the tree is, rather, the characteristic that defines the biomass. Such a characteristic can be measured by Lidar, and we can speculate that, with Lidar-derived height, an estimation of the biomass per tree and per plot would be feasible.

Figure 3 presents the statistics for the field plots, depending on the classification of each plot as monodominant or not. Monodominant plots inherit a smaller variation and mean biomass. The smaller variation in monodominant plots is anticipated and explained from the homogeneity of the species therein. It is noticeable that there is one outlier in the monodominant plot with a considerably higher biomass than any other sampled plot. This is attributed to the existence of mahogany, which is an exotic tree species for Indonesia (being native to the Caribbean), growing up to 30 m high and with a large spherical crown and many heavy branches [100], facts which indicate a species with a high biomass. A closer look at the field data reveals that the only four plots with biomasses higher than 150 ton/ha are all comprised of mahogany as the dominant species.

When investigating the statistical distribution of the plot biomass as a function of the dominant species (Figure 4a), it becomes apparent that mahogany-dominated plots contain a higher amount of biomass. Moreover, when separating between monodominant and non-monodominant cases (Figure 4b), the same trend is observed; the single monodominant mahogany plot exhibits a considerably higher value of biomass. The biomass levels for acacia and teak are at similar levels; however, the standard deviation for acacia is considerably lower for both cases of monodominant and non-monodominant plots.

The current study focused on investigating the relationship between the plot biomass and the Sentinel-1 and Sentinel-2 band combinations. The R² map for each experiment was derived and the statistically most important scenarios are presented in Figure 5, Figure 6 and Figure 7. Figure 5 depicts the results when considering the entire dataset of the 30 field plots. It is evident that there are two spectral regions where the correlations seem to be significant. Specifically, the Ri/Rj combinations of bands 2–5 and bands 6–8A correlate with R² > 0.5 (max R² = 0.70 for the band combination 2/8). This is an indication that spectral bands constructed based on the visible and near-infrared bands correlate strongly for this private forest. Interestingly, there is a second spectral region at the combined bands of the near-infrared with the SWIR bands (i.e., 11 and 12) which also yields similar correlation (0.62 ≤ R² ≤ 0.66). Both of these spectral windows indicate a low p-value and therefore are statistically significant. This is in agreement with a study which also utilized the complete band combinations of Sentinel-2 to estimate AGB in a boreal forest [44], which found a low correlation (R² = 0.24) based on the NDSI combination of bands 8a and 10, 6 and 11, and 3 and 12; nevertheless, all band combinations had similarly low values and the aforementioned bands cannot be considered as a suggestion of AGB proxy. Relatedly, band 8a was also the feature of all bands and indices with the highest relative importance in a prediction model that Pandit et al. [19] used to estimate AGB in a sub-tropical forest in Nepal. Moreover, it is worth noting that the R² of the combinations of the visible and SWIR regions overall perform poorly and are, for all cases, below 0.5.

The combination of the Sentinel-2 VH and VV bands does not yield high correlations. Moreover, the combination of the Sentinel-1 bands the Sentinel-2 bands perform adequately for only two combinations, namely that of Sentinel-1 VH with the red-edge and NIR bands of Sentinel-2 (R² = 0.53 and R² = 0.58, respectively). Incorporation of Sentinel-1 and Sentinel-2 bands has been shown to yield higher accuracies than the stand-alone sensors’ bands [37,52]; nevertheless, that was not the case in the current study. In Figure 5, it is actually evident that information from Sentinel-2 is a better predictor of AGB in comparison to Sentinel-1, which concurs with a study mapping AGB in a West African dryland forest from Forkuor et al. [36]. It is perhaps the case that more sophisticated algorithms for incorporating information from SAR and optical satellites can yield higher accuracies and that the spectral indices used in the current study cannot adequately capture the complementary information from these sensors, resulting in lower R² values when using the Sentinel-1 and 2 datasets synergistically. Lastly, the RSI spectral indices did not yield results differing considerably from the NDSI; hence, the latter are not presented.

When considering only the field samples for which teak is the dominant species (Figure 6), the same spectral windows with the full dataset (Figure 5) seem to appear, only with higher correlations at the band combinations 2, 3 and 6, 7, and 8A (0.70 ≤ R² ≤ 0.74). The highest correlation found for the teak-dominant samples was R² = 0.74 for the band combination 3–7. The fact that we considered plots with a single dominant species seems to increase the correlation. Moreover, the combination of the Sentinel-1 and NIR, which demonstrated high correlation for the entire dataset, is not present in the teak dataset; in this case, the Sentinel-1 yielded high correlation when coupled with the visible (bands 2–5 with R² = 0.62 for band 4) or the SWIR bands (R² = 0.58). For the other 2 dominant species, namely mahogany and acacia, there was no statistically acceptable correlation achieved (p < 0.001); hence, the results are not presented herewith. This is perhaps attributed to the small size of the dataset, as most of the sample plots are primarily planted with teak.

In the third scenario (Figure 7), the role of monodominance was examined. The generic pattern observed is similar to the results when using the entire dataset, especially for the case of monodominance; nevertheless, the abundance of band combinations with R² > 0.5 are comparatively limited. A strong correlation can be observed for the combinations of the band 2 and bands of the NIR (6,7,8,8A) (0.71 ≤ R² ≤ 0.78), with the ratio 2/7 having the highest correlation (i.e., R² = 0.78) and yielding the highest correlation for all the scenarios of this study, which is perhaps justifiable by the fact that only monodominant plots were taken into consideration in this subset of results.

Figure 8 depicts the linear regression for the highest correlation of the monodominant plots scrutinized in this paragraph, as well as the highest correlation for the entire dataset; it is worth noting that there does not seem to be a specific species which performs better than others, which is perhaps attributable to assigning species-tailored allometric equations for the estimation of plot biomass.

The non-monodominant dataset depicts a similar spectral window for the Sentinel-2 band combinations; band 2 correlates with bands 6–8A (0.62 ≤ R² ≤ 0.69). Moreover, the non-monodominant dataset indicates the band candidate of Sentinel-1 VH and Sentinel-2 band 8 as a cross-satellite combination correlating to AGB (R² = 0.66), which is the exact same combination highlighted in the entire dataset scenario (Figure 5); this band combination, and actually any combination involving Sentinel-1 bands, does not appear in the cases of the monodominant samples in Figure 7, and it is a remaining scientific question whether or not the SAR information from Sentinel-1 is actually able to correlate more strongly with heterogeneous plots than homogeneous plots.

In order to evaluate the relevance of our findings in the context of traditional methods, popular vegetation indices were estimated and similar statistics were extracted.

Table 4. Empirical vegetation indices extracted from the Sentinel-2 satellite image and their respective correlations to the plot biomass.

Vegetation Index	Adjusted R2	p-Value	Adjusted R2 (45 Plots)	p-Value (45 Plots)
band8/band2—current paper	0.7	<0.001	0.57	<0.001
EVI	0.65	<0.001	0.51	<0.001
ARVI	0.61	<0.001	0.45	<0.001
mNDVI	0.6	<0.001	0.44	<0.001
VOG1	0.55	<0.001	0.39	<0.001
SAVI	0.53	<0.001	0.42	<0.001
NDVI	0.53	<0.001	0.42	<0.001
MSAVI2	0.51	<0.001	0.40	<0.001
CRI	0.29	0.0014	0.22	0.0014

presents the results of the empirical indices and the correlations found.

The best performing popular index was EVI (R² = 0.65) followed by ARVI (R² = 0.61) and mNDVI (R² = 0.60). Compared to the SI suggested in the current study, all popular vegetation indices underperformed. As a matter of fact, there were 7 band combinations identified in Figure 5 with an R² higher than the one of EVI and 19 band combinations with R² higher than the one of mNDVI. It is worth commenting on the band variation of these three best performing popular indices; EVI is based on the NIR, red, and blue bands; ARVI is based on the red edge, red, and blue bands; and mNDVI is based on two red-edge and blue bands. The only common denominator in these three indices is the use of the blue band, which was also found in the best performing RSI index that the current study suggests.

Regarding the red-edge bands, the mNDVI relies heavily on two red-edge bands, and the VOG1 relies exclusively on a combination of different red-edge bands. These indices provided values of R² = 0.60 and R² = 0.55, respectively. The fact that biomass correlates sufficiently with a vegetation index that uses exclusively red-edge bands, such as the VOG1, is encouraging for the usefulness of the red-edge bands for vegetation trait retrieval. Askar et al. [20], in a study conducted at the same area, estimated AGB based on Sentinel-2 images and reported that the Normalized Difference Index 45 (NDI45), which is constructed from red and red-edge bands, correlated strongly with AGB compared to other indices (R² = 0.79). Other studies have also emphasized the importance of the red-edge band of Sentinel-2; for instance, Delegido et al. [101] reported improvement of accuracy in estimating chlorophyll when including the bands in the red-edge region.

Another important point is that the NDVI, SAVI, and MSAVI2 are indices constructed based on the red and NIR bands, and these three exhibit similar results with a slightly lower accuracy than the group of three indices discussed in the previous paragraph. The reason for this might be the fact that the aforementioned type of indices saturate at high values of biomass which is a fact in tropical environments. Lu et al. [102], in a study comparing Landsat-based NDVI and SAVI, suggested that the vegetation indices that include SWIR spectral bands have higher correlations with AGB in the case of forest stand with complex structure, while the vegetation indices that incorporate NIR associate more strongly with forest stand with simple structure. As an example, Askar et al. [20] suggested that substituting the NIR band with the red edge on the NDI45 index improved the prediction of the biomass, a finding which was seconded by Frampton et al. [103].

Biomass retrieval from remote sensing data from several sensors has been demonstrated frequently in the scientific literature. Spatial resolution is one of the satellite sensor’s specifications that matters when quantifying biomass from space. For instance, Naidoo et al. [37], in a study examining herbaceous AGB of a wetland in South Africa, reported comparatively higher AGB prediction accuracies for WorldView-3 than Sentinel-2, which have spatial resolutions of 1 m and 10 m, respectively. Similarly, the spectral resolution is of importance. For example, Sibanda et al. [43], in a study in which they quantified grass AGB across different fertilizer treatments in South Africa, while they determined that hyperspectral-resampled Sentinel-2 and Landsat-8 data provided a good correlation (R² = 0.81 and R² = 0.76, respectively), hyperspectral data nevertheless exhibited a better correlation (R² = 0.92). A similar conclusion was reached by Sibanda et al. [104], who found that the accuracy of Sentinel-2 was comparatively, but acceptably, lower than that of the HyspIRI hyperspectral sensor (R² = 0.58 and R² = 0.69, respectively). Sentinel-1 and Sentinel-2 satellites have brought new capabilities in this domain, as they offer a combination of very high radiometric resolution at key spectral regions for vegetation-related information. The current study reinforces the usability of Sentinel-1 and Sentinel-2 satellites in retrieving AGB remotely.

Finally, regarding the experiment on the robustness of the methodology, we reran the analysis based on 1617 tree samples collected from 45 plots between September and November 2017. The R² graph (Figure 9) presents a similar pattern with the results from the 30 plots (Figure 5), which indicates that the proposed approach can be applied with a more diversified dataset. However, the correlation for each of the band combinations was reduced, which can be attributed to the temporal dimension of this dataset; the field data collected in November 2017 were collected four months apart from the acquisition of the Sentinel-2 satellite image. Nevertheless, this reduced correlation was homogeneous across both the full combination of spectral indices and the traditional indices calculated. As presented in Table 4, all empirical indices had a decrease in R² ranging from 0.07 for the case of CRI up to 0.16 for ARVI, mNDVI, and VOG1. The best-performing combination in this study, band8/band2, had a reduced R² of 0.13 compared to the 30 plot stable experiment and was still superior to all the traditional vegetation indices. The above results showcase the applicability of our approach, and it is of great importance to conduct further experiments with cloud-free satellite images in the future to thoroughly investigate the robustness of the proposed methodology.

4. Conclusions

Satellite remote sensing has been an indispensable tool in the retrieval of vegetation traits. The application of traditional vegetation indices has long been the norm for translating satellite images into tangible information and subsequently using these indices as proxy to biophysical properties. The advent of contemporary satellites has made available augmented information in key spectral bands, such as ultraviolet, red edge, and thermal, and at increased numbers of bands and radiometric accuracies, opening the way for new methods of approximating biophysical traits. In the current case study conducted in a tropical private forest, we demonstrated that the exhaustive combination of the simple spectral ratio between all the Sentinel-1 and Sentinel-2 bands yields considerably higher correlation compared to nine traditional vegetation indices. The spectral ratio from bands 8 and 2 from Sentinel-2 yielded an R² = 0.70, while the best-performing popular spectral index, namely EVI, yielded R² = 0.65. When restricting the data pool to only monodominant plots, the R² was 0.78 for the band combination 7/2, a fact which can perhaps be attributed to the homogeneity of the field specimens. Overall, we claim that an investigation of the complete band combinations for small scale studies can provide a stronger correlation compared to the random application of popular vegetation indices and circumvent the assumptions on the environmental, structural and ecosystem characteristics that are made for the selection of a vegetation index.