Estimating the Aboveground Fresh Weight of Sugarcane Using Multispectral Images and Light Detection and Ranging (LiDAR)

Abstract

This study aimed to develop a remote sensing method for estimating the aboveground fresh weight (AGFW) of sugarcane using multispectral images and light detection and ranging (LiDAR). Remotely sensed data were acquired from an unmanned aerial vehicle (drone). Sample plots were harvested and the AGFW of each plot was measured. Sugarcane crown heights and volumes were obtained by isolating individual tree crowns using a LiDAR-derived digital surface model of the area. Multiple linear regression (MLR) and partial least-squares regression (PLSR) models were tested for the field-sampled AGFWs (dependent variable) and individual canopy heights and volumes, and spectral indices were used as independent variables or predictors. The PLSR model showed more promising results than the MLR model when predicting the AGFW over the study area. Although PLSR is well-suited to a large number of collinear predictor variables and a limited number of field samples, this study showed moderate results (R² = 0.5). The visual appearance of the spatial distribution of the AGFW map is satisfactory. The limited no. of field samples overfitted the MLR prediction results. Overall, this research highlights the potential of integrating remote sensing technologies in the sugarcane industry, thereby improving yield estimation and effective crop management.

1. Introduction

Sugarcane (Saccharum spp.) is one of the largest broadacre crops in the world. It is grown in 106 countries across the tropics and subtropics and accounts more than 70% of global sugar consumption [1]. Additionally, sugarcane plays a fundamental role in renewable energy production, ranking as the second major feedstock for ethanol production [2]. Moreover, the crop is identified as the most efficient bioenergy crop in tropical regions, and biotechnological tools have been emphasized in maximizing its potential for sustainable energy crop development [3]. These factors underscore the economic and scientific significance of this crop. As sugarcane remains an economic mainstay for many regions of the world, the indirect estimation of sugarcane yields is of crucial importance for the induction of this crop [4].

Sugarcane cultivation is one of Costa Rica’s main agricultural products and plays a key role in boosting the national economy; for example, the sugarcane sector, through production and agroindustry, contributes to 0.33% of the gross national product and 3.83% of the gross agricultural domestic product [5]. The sugarcane industry generates 25,000 direct jobs and more than 100,000 indirect jobs [6]. Given the significant importance of sugarcane cultivation in the agricultural sector in Costa Rica, crop yield assessment plays a vital role. The accurate estimation of factors such as the aboveground fresh weight of crops before harvesting enables effective crop management while increasing the capacity for informed decision-making [7]. For instance, knowledge of the aboveground fresh weight of sugarcane allows producers to predict the amount of raw material to be harvested and, accordingly, to plan production and planning processes more efficiently in terms of managing resources and optimizing production and decision-making [8].

Traditionally, the crop yield has been determined by directly weighing the harvest from a sample of the farm and averaging the whole field. However, this is costly and laborious at the production level for large areas; therefore, there is a demand for developing an alternative method that incorporates the use of new technologies [9,10]. Predicting the sugarcane yield is not an easy task, although some studies have explored alternatives to improve this process. A study carried out in Brazil revealed significant comparative results between data obtained from multispectral images (red, green, and blue bands (RGB images)) acquired by unmanned aerial vehicles (UAVs) and field observations of the aboveground yield [11]. The study aimed to assess the degree of canopy closure in order to predict the potential yield. The crop growth was evaluated by considering the Leaf Area Index (LAI) and Green–Red Vegetation Index (GRVI). The results revealed that the GRVI (R² = 0.69) provided better results compared to the LAI (R² = 0.34) when predicting the yield [11]. Sumesh et al. [8] estimated the sugarcane yield using a consumer-grade red–green–blue (RGB) camera mounted on a UAV with a minimal field dataset. The yield was calculated by integrating the spatial variability of plant height, stalk density, and weight using ordinary least-squares regression. The study outlined a very-high-spatial-resolution RGB image and object-based image analysis (OBIA) approach for mapping the spatial variability of plant height and stalk density and estimating the sugarcane yield [8]. Another study conducted in Thailand demonstrated higher accuracies of OBIA (92% and 96%) in contrast to pixel-based methods (84% and 88%) for estimating sugarcane coverage [12]. After isolating the sugarcane coverage, the number of stems and the average height and weight were calculated to assess the yield [12]. The results indicated the effectiveness of applying OBIA techniques for estimating the sugarcane aboveground fresh weight [12]. Aleman et al. [13] estimated the sugarcane yield by integrating aerial photographs, photogrammetric methods, and field determinations. Digital surface models (DSMs) and digital terrain models (DTMs) were obtained, and based on those, the height and volume of the sugarcane were calculated. These estimates were correlated with different yield variables, such as metric tons of sugarcane per hectare. There was a high linear correlation between this variable and the height of the crop (R² = 0.83). The absolute error of the yield estimation from plant height was 3.7 ton/ha. However, although this study provided promising results, sugarcane spectral information was not considered in the analysis. Therefore, demand still exists for analyzing sugarcane’s spectral information and thus estimating the aboveground fresh weight.

Some studies have evaluated sugarcane crop growth stages using light detection and ranging (LiDAR) technology. Villareal and Tongco [14] achieved 98.4% accuracy for classifying sugarcane crop growth stages and 91.7% for detecting the maturation stage. These results support the usefulness of developing sets of rules, suggesting their applicability to map and verify other sugarcane plantation areas. Furthermore, as LiDAR data easily capture information on the vegetation structure, they are suitable for predicting the aboveground fresh weight and generating spatial distribution maps accordingly [9,15]. LiDAR data provide more consistent and significant correlations of sugarcane biophysical parameters, especially volume and height, than multispectral images [15]. Aerial images provide quantitative information about the vegetation and its health conditions through spectral data analysis, reflecting the expected amount of yield [16]. Hence, an integrated approach combining these two data types might provide accurate statistics on the sugarcane aboveground fresh weight. Although existing studies have demonstrated the effectiveness of individual remote sensing technologies (multispectral imaging or LiDAR) for predicting the aboveground fresh weight (AGFW), according to the authors’ knowledge, there are very limited or no studies integrating these technologies, especially in Costa Rica, where we conducted this research.

The aim of this research was to develop a remote sensing method for estimating sugarcane AGFW to streamline production at the farm level. The proposed data included multispectral images and LiDAR captured from a UAV as well as field observations (AGFW of sample plots) for calibration and validation. The specific objectives were to (1) delineate individual sugarcane tree crowns using LiDAR; (2) develop a mathematical model that estimates the AGFW of sugarcane using spectral values from remotely sensed imagery and field observations; and (3) apply the developed model to predict the AGFW over the entire study area.

2. Materials and Methods

2.1. Study Area

This study was conducted over a sugarcane field located in Guanacaste Province, Costa Rica (10°20′9″ N, 85°9′31″ W) (Figure 1). The farm is called “Agrícola El Cántaro”, and the cultivated sugarcane variety was “CP722086”. The average temperature in this area is consistently high, averaging 28.3 °C with a maximum of 40.1 °C and a minimum of 15.9 °C, which makes it a suitable region for sugarcane cultivation [5]. The summer months span from December to April, typically the sugarcane harvesting months, and coincide with the resowing of sugarcane during the first week of January [17].

The field was in its fourth-year reproduction cycle and was eight months old at data acquisition (in the maturation stage). According to the farm management, this plot was homogeneous and productive compared to other plots within the neighborhood. The farm utilizes a gravity irrigation system and implements scheduled amendments to the crops [18].

2.2. Data Acquisition and Preprocessing

The aerial images and LiDAR and field observations were collected on 24 August 2023. A UAV: DJI Matrice 300 [19] equipped with a Real-Time Kinematic (RTK) differential correction system, Zenmuse L1 laser scanner [20], and Micasense Red Edge-P camera [21] were used for data acquisition. The flying height was selected as 60 m aboveground for LiDAR data acquisition. To capture multispectral images with a high spatial resolution (5 cm), the flying height was 80 m above ground level. The LiDAR point cloud was converted to the LAS format for further processing, as a standard binary format is recommended for industry usage [22]. LiDAR data were cleaned by removing the data noise, and a semi-automatic classification of low points was performed, considering variations in iteration angle, iteration distance, terrain angle, and terrain slope [23]. This process was carried out using Spatix software (Version 023.002) [24].

The camera provided multispectral images with five bands and a panchromatic image.

Table 1. Characteristics of the Micasense RedEdge-P images (adapted from https://ageagle.com/drone-sensors/rededge-p-high-res-multispectral-camera/; accessed on 10 September 2023).

Band	Wavelength (nm)	Bandwidth (nm)
Blue	475	32
Green	560	27
Red	668	14
Red edge	717	12
Near infrared (NIR)	842	57
Panchromatic	634.5	463

shows the spectral sensitivity (center wavelength) of each band. Multispectral images were processed to obtain the orthomosaic (with spectral reflectance values) for further analysis using Pix4DMapper software (v 4.9) [25].

Five sampling plots (3 m² each) were identified based on the green appearance (homogeneity) of the plants, the spatial distribution, and the access to the plots. The location of each plant on each sampling plot was captured using a Global Navigation Satellite System with an RTK Fix solution. The field measurements showed that the distance between sugarcane rows on the farm was not uniform; however, the average distance between rows was 0.9 m, and the average distance between plants was 0.65 m. The total weight of millable stalks (aboveground fresh weight (AGFW)) in each harvest plot was measured; 70% of sampling was used to create the models, and 30% was used for validation. Figure 2 shows the complete workflow for this study.

2.3. Individual Tree Crown Delineation

A digital terrain model (DTM) and the digital surface model (DSM) were extracted from the LiDAR data (Figure 2). A DTM contains only information about the terrain surface, excluding objects and surface features, and a DSM consists of an elevation dataset that includes both terrain information and information about the objects and surface features present [26]. The canopy height model (CHM) was created for a spatial resolution of 0.2 m using Equation (1).

$$CHM=DSM-DTM$$

(1)

The individual tree crowns were isolated using the “lidR” package in R software [27]. First, the CHM’s highest points (local maxima) were extracted as sugarcane tree top locations. A threshold value of 0.5 m was applied to avoid small plants and stalks in this calculation. Then, the algorithm developed by Dalponte and Coomes [28] was used to extract individual tree crowns. The approximate crown diameter was selected as 1.7 m by manually measuring several crowns on-screen within the study area. Heights of individual trees and the canopy coverage were extracted (Figure 2).

2.4. Estimating the Aboveground Fresh Weight

As shown in Figure 2, vegetation indices were derived from the multispectral image. A vegetation index transforms the spectral values in a meaningful way to enhance the contribution of vegetation properties and thus ease the identification of vegetation vigor [29]. Nine vegetation indices (

Table 2. Details of vegetation indices selected for this study.

Name	Equation	References
Green Normalized Difference Vegetation Index (GNDVI)	$\frac{NIR - Green}{NIR + Green}$	[34]
Normalized Difference Red Edge Index (NDRE)	$\frac{NIR - Red Edge}{NIR + Red Edge}$	[35]
Normalized Difference Vegetation Index (NDVI)	$\frac{NIR - Red}{NIR + Red}$	[36]
Soil-Adjusted Vegetation Index (SAVI)	$\left(\frac{\left(NIR - red\right)}{\left(NIR + red+L\right)}\right)\times \left(1+L\right);where L = 0.5$	[31]
Modified Soil-Adjusted Vegetation Index (MSAVI)	$\frac{2 x NIR + 1-\sqrt{{\left(2 x NIR + 1\right)}^2-8 x\left(NIR - Red\right)}}{2}$	[37]
Simple Ratio (SR)	$\frac{NIR}{Red}$	[38]
Normalized Green–Red Difference Index (NGRDI)	$\frac{Green - Red}{Green + Red}$	[33,37]
Normalized Green–Red Edge Difference Index (NDGRed)	$\frac{Green - Red Edge}{Green + Red Edge}$	[37]
Triangular Greenness Index (TGI)	$-0.5\left[\left({\lambda }_r-{\lambda }_b\right)\left(Red - Green\right)-\left({\lambda }_r - {\lambda }_g\right)\left(Red - Blue\right)\right]$*	[33,37]

* λ_r, λ_b, and λ_g represent the center wavelengths for the red, blue, and green bands, respectively.

) were selected for this study based on their ability to detect the greenness of vegetation and biomass [30]. For instance, the Normalized Difference Vegetation Index (NDVI) is well-known for analyzing sugarcane biomass [10,12]. Additionally, there were some yield gaps in the study area due to replanting (small plants—early stages and bare soil), so the SAVI and MSAVI were selected considering their effectiveness in this situation, especially when there is abundant soil and small plants [31]. Meanwhile, the positive values of the Triangular Greenness Index (TDI) correspond to green vegetation and thus the leaf chlorophyll content [32], while the negative values may represent soil. Since the TDI is calculated using blue, green, and red wavelengths, it is suitable for RGB cameras in agricultural applications such as this one [33].

Once vegetation indices were created, two algorithms, multiple linear regression (MLR) and partial least-squares regression (PLSR), were tested to calculate the AGFW. For the MLR method, the correlations between the AGFW of each plot and independent variables were tested to select the best predictor variables. Each plot’s canopy cover, mean canopy height, volume, and mean vegetation indices (NDVI, NDRE, GNDVI, SAVI, MSAVI) provided the highest correlations with field samples. Equation (2) shows an example of the MLR model.

$$Y={\beta }_0+{\beta }_1{X}_1+{\beta }_2{X}_2+\cdots +{\beta }_n{X}_n+ϵ$$

(2)

where Y is the dependent variable to be predicted, $X_1, X_2, \mathrm{and} X_n$ are the independent variables, ${\beta }_0$ is the intercept, ${\beta }_1, {\beta }_2, \mathrm{and} {\beta }_n$ are the coefficients for each independent variable, and ϵ is the error.

The accuracy of the model was assessed using the coefficient of determination R-squared ($R^2$). $R^2$ is a measure of the goodness of a fit of a model and indicates how the AGFW varies with predictor variables.

Since we had limited AGFW samples, the PLSR algorithm was also selected to predict the AGFW. PLSR is an efficient and optimal regression model based on covariances. The PLSR model generalizes and combines principal component analysis and multiple linear regression, producing a linear regression model [39]. It reduces the independent variables to a smaller set of predictors internally. The algorithm adds complexity and utility to the model, as it is designed to deal with situations where there are many independent variables, possibly correlated, and relatively few samples as dependent variables [40].

Random points were generated inside individual tree crowns of each sample plot to ensure the homogeneity of spectral values. Corresponding spectral band values (Table 1 except panchromatic band) and different vegetation indices (Table 2), canopy height, and volume were extracted for each point. The initial PLSR model was created using 15 predictor variables, assuming that more variables would improve the prediction power. The “pls” package in R software was used [41]. The model performance and the best number of components were evaluated using cross-validation (root mean square error of prediction (RMSEP)). Once the best number of components was identified, the final model was created using four components only.

Finally, the PLSR model was selected to predict the spatial distribution of AGFW over the study area. A low-pass filter was applied to the resultant raster to remove data noise. The accuracy of the predicted AGFW raster was assessed using R-squared ($R^2$).

3. Results

3.1. Individual Tree Crown Delineation

Individual tree crowns were detected successfully. Their accuracy was assessed visually compared to the multispectral image. For example, Figure 3a shows a part of the study area with delineated tree tops and crowns. However, some areas do not show tree tops or crowns (Figure 3b). These areas have trees shorter than 0.5 m (threshold applied in the algorithm) or soil.

Also, the number of trees in each sample plot was compared with the number of trees delineated in each lot plot (Figure 4). There was a 100% agreement between these data for three plots. In plots 2 and 5, an extra crown was detected, mainly due to the foliage of neighboring plants.

3.2. Estimating Aboveground Fresh Weight

3.2.1. Multiple Linear Regression (MLR)

Table 3. The correlation between aboveground fresh weight and the predictor variable of each plot.

Variable	Correlation Coefficient (r)
Canopy area (m2)	0.86
Canopy height (m)	0.63
NDVI	0.72
NDRE	0.67
GNDVI	0.63
SAVI	0.91
MSAVI	0.95
Volume	0.16

shows the correlation between the AGFW of each sample plot and the mean value of each predictor variable. All predictor variables where we obtained a correlation (r) greater than 0.7 were selected for the final regression model—the MSAVI, SAVI, NDVI, and crown area (Table 3).

Figure 5 shows the R² values of each chosen predictor variable. In this study, the NDVI had a lower R² than the MSAVI and SAVI. Although we expected a high correlation between the crown area, canopy height, and aboveground fresh weight, they showed a moderate relation. The final MLR model showed a high correlation between sample data and selected predictors (R² = 0.7).

3.2.2. Partial Least-Squares Regression (PLSR) Model

As shown in Figure 6, the lowest RMSEP was obtained with only four components. Although three to seven components were within the lowest RMSEP range, increasing the number of components did not improve the prediction performance. The cross-validation results also confirmed the variability of RMSEP with a changing number of components.

3.3. Spatial Distribution of Aboveground Fresh Weight

The prediction results of the MLR model were unreliable and exhibited over prediction. The resultant AGFW over the study area ranged from 200 kgm⁻² to 300 kgm⁻², indicating model overfitting.

The PLSR model predicted the AGFW over the study area (Figure 7a). The results ranged from 33.1 kgm⁻² to 70.0 kgm⁻². The R² compared to field samples was 0.5, indicating moderate prediction results. The AGFW showed higher values closer to areas with high precipitation, especially south of the plot, adjacent to a large river. Lower values were exhibited closer to the access road, where the ground elevation was higher than that of other areas. As such, the visual analysis of the map indicated a reasonable pattern related to the geography and availability of water in the study area. Figure 7b illustrates the distribution of the mean AGFW of each crown in the study area, and the red density curve fitted to the data reflects normality in the distribution of these and the concentration of values around the mean, which has a value of 47 kgm⁻². Overall, the study area has 11.48 ha and an aboveground fresh weight of 184 ton/ha.

4. Discussion

4.1. Individual Tree Crown Delineation

This study tested the tree crown delineation method developed by Dalponte and Coomes [28] for forest tree crowns. Since it was successful in the sugarcane environment with some customization, we demonstrate the method’s applicability in sugarcane studies and its advantages over image segmentation. For instance, when identifying individual tree crowns using aerial image classification algorithms (pixel-based or object-based) without incorporating CHM, tree heights are not considered. Hence, small, juvenile plants are identified as mature trees in the calculation of the potential yield. Pixel-based image analysis has long been the mainstay approach for classifying remotely sensed imagery due to its simplicity, effectiveness in numerous studies, and wide range of applications [42]. However, the pixel-based classification experienced difficulties in segmenting sugarcane and non-sugarcane (noise, soil, and grass). It cannot classify the leaf structure and canopy due to limitations when considering texture and spatial structure only, leading to erroneous stem counting [12]. Object-based approaches offer a significant improvement in this regard, integrating a wide range of features that go beyond spectral information, such as shape and texture, resulting in a more comprehensive representation of the study areas [43]; however, the latter require more detailed parameterization and can be more susceptible to errors in heterogeneous environments [12]. Nonetheless, this study has demonstrated the appropriateness of the method developed by Dalponte and Coomes [28] to overcome these challenges in sugarcane studies.

4.2. Estimating Aboveground Fresh Weight

The results indicated a significant correlation between some predictor variables and aboveground fresh weight; however, some showed high multicollinearity, especially between the SAVI and MSAVI. This occurred because these indices use reflectance from the same set of bands but are based on different mathematical combinations. Canopy volume and height did not correlate well with the AGFW samples, suggesting they may not be as descriptive as vegetation indices based on spectral reflectance. Furthermore, it is well-known that the vegetation indices indirectly represent plants’ biophysical variation [44]. However, this could be influenced by the variability of the sample plots and the procedure to measure each plot’s aboveground fresh weight. For instance, once each plot was harvested, the top part of each plant (leaves) was removed before weighing. The mean canopy volume and height calculated from LiDAR included that information. Hence, we recommend removing the corresponding slice of LiDAR data points from the top before calculating canopy volumes and heights.

Although the MLR model performed well (R² = 0.737) relative to sample data, the prediction results were unreliable. The R² was relatively lower than those in existing studies that mapped several crop yields using MLR models. For example, Niedbała [45] used the MLR model to predict the winter wheat yield, achieving between 0.518 and 0.5457, indicating moderate results. The study tested three MLR models with large number of predictor variables (13, 17, and 19) and field samples. Abrougui et al. [46] used a large number of field samples and achieved high performance with an MLR model (R² = 0.89) for agricultural yield prediction. Aleman et al. [13] achieved an R² of 0.83 using images acquired from a drone and MLR algorithm when estimating the sugarcane yield in metric ton/ha; their observation was that the model could present a better fit if the training data were more varied and larger in quantity. However, the study’s limited number of field samples made it difficult to conclude on the performance. The model was overfitted with a limited number of field samples, as MLR models are known to be sensitive to the number of sample points [47].

Partial least-squares regression reduces the given predictors to a set of uncorrelated components, which will be used for least-squares regression instead of original data [48]. Hence, PLSR is identified as a fast, efficient, and optimal regression method that provides a high performance based on covariance [40]. With four components, this study’s PLSR model provided moderate prediction results. This was mainly due to the multicollinearity of predictor variables and the limited number of samples. However, PLSR is primarily recommended when there are more predictor variables than the number of samples [40]. Although PLSR is commonly used in the chemical, drug, food, and plastic industries [48], increasing applications can be seen in natural resources management and agriculture. Lopez et al. [49] obtained an R² of 0.66 and RMSE of 10% when estimating yield losses using the PLSR model, with results similar to those in this study. Jeong et al. [50] indicated that overfitting can be a problem when the training data are not sufficiently varied or large enough, which highlights the importance of a well-designed training dataset that covers a wide range of predictor variability to minimize this problem. Xu et al. [9] showed a significant positive correlation between the height and weight of sugarcane field samples. Their prediction results were also promising.

To improve the accuracy of these models, we recommend increasing the number of samples and considering the spatial variability of sugarcane when selecting plants. Sanches et al. [11] suggested the inflection point for biomass accumulation by sugarcane as a good indicator for use to estimate the yield. Hence, we recommend combining images from different growing stages (multitemporal) to capture the spectral variations over time and decide on the inflection point, thus potentially enhancing the robustness of the results.

5. Conclusions

Sugarcane is one of Costa Rica’s main agricultural products and is key in boosting the national economy. Traditionally, the crop yield has been determined by directly weighing the harvest from a sample of the farm and averaging the whole field. However, this is costly and laborious at the production level for large areas; therefore, there is a demand for developing an alternative method that incorporates the use of new technologies such as remote sensing. Hence, this research aimed to develop a remote sensing method for estimating the sugarcane aboveground fresh weight to streamline the production at the farm level.

This study used multispectral images acquired from a Micasense Rededge-P camera and LiDAR data. The aboveground fresh weight measured in the field for sample plots were used for model calibration and validation. The digital terrain and digital surface models extracted from LiDAR data were used to create a canopy height model (CHM). Individual sugarcane tree crowns were delineated using the method developed by Dalponte and Coomes [28]. The extent of each crown area was calculated. The average height of each crown was extracted from the CHM. Two regression models (multiple linear regression (MLR) and partial least-squares regression (PLSR)) were developed to check the correlation between aboveground fresh weight samples and vegetation indices derived from images. The MLR model obtained higher accuracy (R² = 0.7) compared to the field samples of aboveground fresh weight. However, the PLSR model showed more promising results than the MLR model when predicting aboveground fresh weight over the study area. The accuracy of the prediction raster was calculated as R² = 0.5 compared to field samples.

This research has highlighted the potential of integrating remote sensing technologies into precision agriculture, thereby improving yield estimation and effective crop management. However, this study faced challenges when assessing the model accuracy and the prediction performance. Despite the successful adaptation of the individual tree crown delineation method, there were limitations mainly due to the limited number of samples collected in the field. Also, there was low correlation between canopy volume, height, and the aboveground fresh weight of samples. This could have been influenced by the variability of the sample plot and the procedure for measuring each plot’s aboveground fresh weight. For instance, once each plot was harvested, the top part of each plant (leaves) was removed before weighing. We recommend removing the corresponding slice of LiDAR data points from the top before calculating canopy volumes and heights. Furthermore, we suggest combining images from different growing stages (multitemporal) to capture the spectral variations over time and potentially help decide on the inflection point of sugarcane biomass production, to enhance the robustness of the results. Lastly, diversifying training samples will help address challenges related to multicollinearity and prediction performance of model accuracy.