Can Multi-Temporal Vegetation Indices and Machine Learning Algorithms Be Used for Estimation of Groundnut Canopy State Variables?

: The objective of this research was to assess the feasibility of remote sensing (RS) technology, specifically an unmanned aerial system (UAS), to estimate Bambara groundnut canopy state variables including leaf area index (LAI), canopy chlorophyll content (CCC), aboveground biomass (AGB), and fractional vegetation cover (FVC). RS and ground data were acquired during Malaysia’s 2018/2019 Bambara groundnut growing season at six phenological stages; vegetative, flowering, podding, podfilling, maturity, and senescence. Five vegetation indices (VIs) were determined from the RS data, resulting in single-stage VIs and cumulative VIs ( ∑ VIs). Pearson’s correlation was used to investigate the relationship between canopy state variables and single stage VIs and ∑ VIs over several stages. Linear parametric and non-linear non-parametric machine learning (ML) regressions including CatBoost Regressor (CBR), Random Forest Regressor (RFR), AdaBoost Regressor (ABR), Huber Regressor (HR), Multiple Linear Regressor (MLR), Theil-Sen Regressor (TSR), Partial Least Squares Regressor (PLSR), and Ridge Regressor (RR) were used to estimate canopy state variables using VIs/ ∑ VIs as input. The best single-stage correlations between canopy state variables and VIs were observed at flowering (r > 0.50 in most cases). Moreover, ∑ VIs acquired from vegetative to senescence stage had the strongest correlation with all measured canopy state variables (r > 0.70 in most cases). In estimating AGB, MLR achieved the best testing performance (R 2 = 0.77, RMSE = 0.30). For CCC, RFR excelled with R 2 of 0.85 and RMSE of 2.88. Most models performed well in FVC estimation with testing R 2 of 0.98–0.99 and low RMSE. For LAI, MLR stood out in testing with R 2 of 0.74, and RMSE of 0.63. Results demonstrate the UAS-based RS technology potential for estimating Bambara groundnut canopy variables.


Introduction
A key strategy for adapting to changing climatic conditions and meeting the increasing global food demand is the development and promotion of underutilised crops [1,2].These crops possess significant potential to enhance food security, diversify agrifood systems, and reduce environmental impacts [3,4].Underutilised crops, such as Bambara groundnut (Vigna subterranea) have evolved specific traits and physiological responses to tolerate harsh environments, including water scarcity and heat stress [5,6].In many regions of Africa, Bambara groundnut is the third most important legume after peanut and cowpea, with an annual production of 300,000 tones [7,8].Despite its economic value, it remains underutilised due to the lack of information on its phenotypic development and performance in different growing environments.Constraints in field phenotyping capability limit plant breeders' ability to dissect the genetics of quantitative crop traits, especially the traits related to yield and stress tolerance.
It is thus essential to monitor at each phenological stage canopy state variables especially those related to yield, such as leaf area index (LAI), canopy chlorophyll content (CCC), aboveground biomass (AGB), and fractional vegetation cover (FVC).LAI is used to assess leaf cover, monitor growth and photosynthesis and provides information on crop health and nutrient status [9].CCC is the main indicator of photosynthesis, senescence, nutritional status, disease, and stress.Real-time monitoring of CCC serves as a guideline for fertilisation [10].AGB is a key parameter that reflects the growth status and is linked to yield and solar energy utilisation [11].FVC is defined as the ratio of the vertically projected area of vegetation to the total surface extent.It is an important biophysical parameter to observe vegetation cover trends and describe canopy vigour.Moreover, FVC is a controlling factor in photosynthesis and transpiration [12].Traditional methods of phenotypic assessments to estimate LAI, CCC, AGB, and FVC are based on manual measurements or visual scoring, both of which are time-consuming, destructive, expensive, and restricted to point estimates which fail to capture the spatial dynamics of the crop growth [13][14][15].
Remote sensing (RS) technology allows fast, non-destructive, and efficient monitoring of crop growth and development [16,17].Moreover, RS technology enables recurrent and systematic information to be obtained from the local to the global scale, thus allowing characterisation of the spatiotemporal variability [18].Unmanned aerial vehicles (UAV) represent effective and low-cost high throughput phenotyping platforms (HTPPs).Their great flexibility and low operational cost make UAV-based HTPPs a promising tool for precision agriculture.Vegetation indices (VIs) extracted from UAV imagery have been used to estimate several biophysical parameters, including LAI [19], AGB [20], FVC [21], and yield [8].In applied agricultural research, the use of consumer grade red, green, and blue channel (RGB) digital cameras is preferred for their simplicity, affordability, and practicality.However, the RGB-based applications are considered inferior especially due to the lack of the near infrared (NIR) band which is highly effective for crop monitoring.Examples of widely used VIs, that incorporate NIR bands, are the normalised difference vegetation index (NDVI) [22,23], simple ratio (SR) [24], green normalised difference vegetation index (GNDVI) [25], enhanced vegetation index 2 (EVI2) [26], and green chlorophyll index (CI green ) [27].
Recent research demonstrated the potential of utilising VIs derived from UAV-based RGB cameras together with machine learning (ML) algorithms to estimate soybean leaf chlorophyll content (LCC), FVC, and maturity [28].Specifically, the reported accuracy metrics for LCC estimation were R 2 = 0.84 and RMSE = 3.99, for FVC estimation were R 2 = 0.96 and RMSE = 0.08, and for maturity monitoring, R 2 was 0.984.Dos Santos et al. [29] utilised a Red-Green-Near-Infrared (R-G-NIR) camera mounted on a UAV to estimate evapotranspiration and AGB of maize crops.They found that the soil-adjusted vegetation index (SAVI) exhibited a stronger correlation with AGB (R 2 = 0.74, RMSE = 0.092 kg m −2 ) compared to NDVI (R 2 of 0.69, RMSE = 0.104 kg m −2 ).Similarly, Ma et al. [30] employed deep learning (DL) techniques with VIs and colour indices (CIs) derived from UAV RGB and colour infrared (CIR) images to estimate rice LAI.The coefficient of determination (R 2 ) for CIs ranged from 0.802 to 0.947, with RMSE ranging from 0.401 to 1.13, while for VIs, R 2 ranged from 0.917 to 0.976, with RMSE ranging from 0.332 to 0.644.
To date, most studies have focused on using RS technology to monitor broad acre crops such as cereals, due to their economic importance and the ease of monitoring their aboveground canopy state variables and yield components [31][32][33].However, leguminous subterranean oilseed crops like groundnuts have been underrepresented in the literature.Groundnuts present unique challenges for yield prediction because their yield components are underground.Current methods rely on destructive sampling and manual inspections, which are labour-intensive and impractical for large-scale applications [34].Assessing factors like pod health, yield components (such as number of pods per plant, seeds per pod, seed size), and quality for groundnuts through surface observations requires complex modelling of canopy state variables, including LAI, canopy chlorophyll content (CCC), AGB, and FVC.While UAV-based RS and ML algorithms have shown effectiveness in monitoring these variables in cereals, their application to subterranean crops remains limited.Thus, this study investigated the use of a UAV-mounted digital camera and ML algorithms to monitor Bambara groundnut canopy state variables at various growth stages.The UAV allows for non-destructive, efficient, and frequent data collection over large areas.Moreover, ML algorithms are essential for processing and analysing large datasets, identifying patterns, and making accurate predictions, crucial for optimising agricultural practices [35].The continuous monitoring of these canopy state variables can help farmers make informed decisions on irrigation, fertilisation, pest management, and other agricultural practices, ultimately enhancing groundnut yield and quality.Finally, this research addresses a critical need in the agricultural sector by providing practical solutions for farmers.

Study Site
The research was conducted at the Field Research Centre of Crops for the Future, located in Semenyih, Selangor, Malaysia.(2 • 55 ′ 56.96 ′′ N, 101 • 52 ′ 33.59 ′′ E), at 560 m above mean sea level from April to September 2018 (Figure 1).The region exhibits a tropical climate, with an average annual temperature of 32.0 °C, an average annual precipitation of 1493 mm, and a typical photoperiod of 12 h day −1 .The environmental conditions recorded by the local weather station at the experimental site and the irrigation rates are shown in Figure 2. The region exhibits a tropical climate, with an average annual temperature of 32.0 • C, an average annual precipitation of 1493 mm, and a typical photoperiod of 12 h day −1 .The environmental conditions recorded by the local weather station at the experimental site and the irrigation rates are shown in Figure 2. The region exhibits a tropical climate, with an average annual temperature of 32.0 °C, an average annual precipitation of 1493 mm, and a typical photoperiod of 12 h day −1 .The environmental conditions recorded by the local weather station at the experimental site and the irrigation rates are shown in Figure 2.

Experimental Design
The field area covered approximately 0.2 ha with sandy clay loam soil with a pH value of 5.23.Three Bambara groundnut genotypes-NAV4 (genotype 1), IITA-686 (genotype 2), and CIVB (genotype 3)-were grown in a randomised complete block design across four blocks.Each block was subdivided into nine plots, with each genotype replicated three times within each block (Figure 1), resulting in 12 replications for each genotype (four blocks × three replicates per block).The gross plot size was 8 m by 7 m (56 m 2 ), while the net plot area was 30 m 2 (6 m by 5 m).Row-to-row (inter-row) spacing between rows was set at 40 cm, and plant-to-plant (intra-row) spacing was 30 cm.The seeding rate for each genotype was 300,000 seeds ha −1 .Prior to sowing, a starter fertiliser was applied at a rate of 20:60:40 kg of nitrogen, phosphorus, and potassium, respectively.Sowing was conducted using a precision planter on 25 April 2018.Fungicides and insecticides were applied at regular intervals to manage pathogens, and weeding was performed manually using hand hoes.Throughout the trial, soil moisture content was monitored weekly using a soil moisture PR2 probe (Delta-T Devices Ltd., Cambridge, UK), and irrigation was initiated when soil water content decreased to 50% of the plant-available water capacity in the root zone.

Agronomic Measurements
Measurements of LAI, FVC, and CCC were conducted at various growth stages-days after sowing (DAS)-namely: vegetative (41 DAS), flowering (58 DAS), podding (84 DAS), pod-filling (97 DAS), maturity (105 DAS), and senescence (114 DAS) from May to September 2018.AGB assessments were made during the vegetative, flowering, podding, and senescence stages.Leaf area was determined using the LI-3100C Area Meter from LICOR (Lincoln, NE, USA), and LAI calculated by dividing the green leaf area by the sampled area.AGB was determined by harvesting all crops within a 1 m 2 area in the central rows at ground level and then drying the clippings for 120 h at 70 • C until a constant weight was achieved.Chlorophyll concentration (Chl) was measured using a SPAD-502 Plus device (Konica Minolta Sensing Inc., Tokyo, Japan) with SPAD units calibrated using the same methodology as described in [36].CCC was estimated by multiplying the Chl per area values by the LAI.For FVC determination, digital images were captured and cropped to the area of interest (AOI) in ImageJ, with cropped images processed using the maximum likelihood supervised classification tool within ArcMap (ArcGIS ® by ESRI Inc., Redlands, AB, Canada).The zonal statistics tool was subsequently utilised to evaluate the number of "vegetated" pixels within each plot, and FVC was calculated by dividing the number of "vegetated" pixels by the total number of pixels in the AOI (refer to Table 1 for summary statistics).

UAV, Sensor and Remote Sensing Data Acquisition Missions
To enhance the reproducibility of our method, we utilised an affordable and readily available commercial UAV, the DJI Phantom 4 Pro (DJI Company, Shenzhen, China; website: https://www.dji.com/,accessed on 10 June 2020).This vertical take-off and landing (VTOL) quadcopter has a maximum payload capacity of 477 g and can sustain flight for 20-30 min, covering distances of up to 7 km.A Canon S100 camera (Canon, Tokyo, Japan), modified by MaxMax (LDP LLC, Carlstadt, NJ 07072, USA; website: www.maxmax.com,accessed on 10 June 2020), was mounted on the UAV using a two-axis aluminium-carbon-fibre gimbal set at nadir view (90 • downward angle).This modification enabled the capture of colour infrared (CIR) digital imagery spanning the Green-Red-NIR spectrum (520-880 nm).
Six flight missions were conducted to obtain high-resolution images during crucial crop growth stages, including vegetative, flowering, podding, pod filling, maturity, and senescence.These flights were carried out under stable weather conditions from 10:00 to 14:00 local time to minimise variations in illumination.The flight paths and settings were predetermined using the DJI Ground Station Pro software 2.0 v.2.0.16, employing a moving-box flight path planning approach.To facilitate future image georectification, four ground control points (GCPs) were permanently positioned at each corner of the field.The accurate GPS coordinates of the GCPs were surveyed using a real-time kinematic (RTK)enabled dual-frequency Leica 1200 Global Navigation Satellite System (GNSS) system with RTK precision (Leica Geosystems, Heerbrugg, Switzerland).Four 2 × 2 m calibration targets with nominal reflectance values of 10%, 20%, 50%, and 80% were utilised for radiometric calibration of the sensor.The spectral reflectance of these targets was measured using a handheld ASD spectroradiometer (FieldSpec, ASD, Boulder, CO, USA).Ultra highresolution images were captured at a pre-scheduled flight altitude of 10 m, in nadir view, with a low flight speed of 0.5 m s −1 , and with an intended overlap of 75% in both in-track and cross-track directions to ensure sufficient overlapping coverage.
A Dell ® Inspiron 7000 laptop (Dell Technologies, Round Rock, TX, USA) was utilised with the Ground Controlling Station software 2.0 v.2.0.16 to manage the autonomous UAV flight via wireless network.The Canon S100 camera, set to TV mode for consistent shutter speed and aperture settings, was autonomously controlled to maintain optimal exposure levels.To automate camera functionality, the Canon CHDK free development software kit version 1.6.1,available at www.chdk.wiki.com,accessed 12 June 2020 was employed.The CHDK script enabled the UAV autopilot system to transmit electronic control signals, automating camera shutter triggering for precise data recording.During the data collection flight, auto-triggering occurred every 3 s (at a frequency of 0.33 Hz), facilitating the capture of approximately 400 images covering all plots within a 20 min flight duration, with a ground resolution of approximately 4 mm per pixel.The captured images were stored in 16-bit local digital memory cards as raw Geographic Tagged Image File Format files for subsequent image processing.

Image Processing
Following the flights, the raw images underwent pre-processing to eliminate electromagnetic interference in the visible bands from the NIR band, using Remote Sensing Explorer Software version 1.0 (MaxMax LDP LLC, Carlstadt, NJ, USA; www.maxmax.com,accessed 10 June 2020).Further pre-processing involved correcting lens distortion, chromatic aberration, and gamma correction using the Digital Photo Professional image processing software version 4.17.20 (Canon Inc., Tokyo, Japan; http://www.canon.co.uk/ support/camera_software/, accessed 10 June 2020).The images were then imported into Agisoft Photoscan Pro Version 1.4.3 (Agisoft LLC, St. Petersburg, Russia) and mosaicked to create a single orthophotomosaic image for the entire study area for each flight date.The pixel size of the orthophotomosaics was approximately 4 mm per pixel to ensure high spatial resolution.Additionally, to minimise band-to-band misalignment, geometric correction was conducted on the orthophotomosaic using GCPs with surveyed GPS coordinates (Section 2.4).The georeferenced orthophotomosaic was subsequently brought into ArcMap Version 10.2.2 (ArcGIS ® , ESRI Inc., Redlands, AB, Canada; https://www.esri.com/en-us/arcgis/about-arcgis/overview,accessed 10 June 2020) for the co-registration of multi-temporal orthophotoimages.Radiometric correction was then applied to each orthophotomosaic for each flight date, band-by-band, to convert raw digital number (DN) values to reflectance values by employing a radiometric calibration equation developed between the known reflectance of the calibration targets and calibration target DN values in the empirical line correction method.Post-processing tasks included subsetting the area of interest (AOI) and digitising the quadrats using the ArcMap Editor tool.Geoprocessing steps also included resampling the images to a consistent pixel grid and correcting for any spatial misalignments between different flight dates.Finally, masks were applied to isolate the central portions of the canopy, excluding borders, and pure canopy pixels were utilised to compute average reflectance.

Vegetation Indices and Cumulative Vegetation Indices Calculation
After preprocessing, a set of five VIs was computed from the remotely sensed data (Table 2).These VIs were calculated for each orthophotomosaic using the ArcMap Raster Calculator tool, ArcMap Version 10.2.2 (ArcGIS ® , ESRI Inc., Redlands, AB, Canada; https:// www.esri.com/en-us/arcgis/about-arcgis/overview,accessed 10 June 2020).Subsequently, the VIs were extracted on a quadrat-by-quadrat basis, and the mean, standard deviation and other statistics were generated using the ArcMap Zonal Statistics tool.Green chlorophyll index (NIR/G) − 1 [27] Note: NIR, R, and G are the reflectance values of the near infrared, red, and green bands, respectively.
We evaluated three integration periods: from flowering to maturity, from flowering to senescence, and from vegetative to senescence.We calculated the ∑VIs over these integration periods using the same formula as in [8].

Model Selection and Modelling Strategy
We carefully selected the models to be used for estimating canopy state variables from VIs/ΣVIs.First, we assessed the dimensionality and multicollinearity of the dataset using principal component analysis (PCA) and variance inflation factor (VIF) analysis.Preliminary PCA results showed that a few principal components (PCs) captured most of the variance, suggesting the need for dimensionality reduction.This supported the decision to select models capable of handling high-dimensional data.High VIF values indicated the need for implementing models, which can deal effectively with multicollinearity.Furthermore, we examined residual plots and conducted tests to identify outliers and heteroscedasticity to inform our decision to select models, which are robust to outliers.Based on these assessments, we selected a diverse set of regression models that can effectively handle high dimensionality, multicollinearity, outliers, and overfitting.
We conducted Pearson's correlation analysis between VIs/∑VIs and canopy state variables.The VIs from the optimal integration interval, based on the correlation analysis, were then utilised as input for both linear parametric and non-linear non-parametric regression models, including Huber Regressor (HR), Theil-Sen Regressor (TSR), AdaBoost Regressor (ABR), CatBoost Regressor (CBR), Multiple Linear Regression (MLR), Random Forest Regression (RFR), Partial Least Squares Regression (PLSR), and Ridge Regression (RR).HR is known for its robustness to outliers, which can indirectly help in datasets where outliers might exacerbate multicollinearity issues [37].TSR, being a non-parametric method based on median calculations, is less sensitive to multicollinearity and offers robustness against certain types of data anomalies [38].ABR, an ensemble method, can enhance model accuracy by combining multiple weak learners, providing a safeguard against overfitting, especially when the base estimator is appropriately chosen [39].CBR, on the other hand, is designed to handle categorical features and high-dimensional data efficiently, with built-in mechanisms to prevent overfitting and address multicollinearity [40].MLR uses regularisation (e.g., Ridge or Lasso) for overfitting, VIF for multicollinearity, and feature selection for high dimensionality [41].RFR minimises overfitting with tree control, handles multicollinearity by subsampling, and addresses high dimensionality by feature selection [42].PLSR reduces overfitting and multicollinearity through latent variables and manages high dimensionality by capturing essential information [43].RR minimises overfitting using L2 regularisation, handles multicollinearity by redistributing variable influence, and reduces dimensionality effectively [44].
Prior to model implementation, the data were preprocessed for missing values, outliers capped, encoded, and normalised using Yeo-Johnson transformation.The dataset was split randomly into 80% training and 20% testing and scaled using StandardScaler from sklearn.preprocessing in a Python 3.9 environment (Python software version 3.9.0,Python Software Foundation, Wilmington, DE, USA).This was determined by experimenting with various splits: 50-50%, 60-40%, 70-30%, 80-20%, and 90-10%.The scikit-learn library was used for models such as HR, TSR, ABR, MLR, RFR, PLSR, and RR, utilising classes like HuberRegressor, TheilSenRegressor, AdaBoostRegressor, LinearRegression, RandomFore-stRegressor, PLSRegression, and RidgeRegression.For CBR, the CatBoostRegressor class from the catboost library was used.Hyperparameters were optimised using GridSearchCV from scikit-learn, which performs an exhaustive search over specified parameter values using cross-validation.Specifically, for each model, a range of hyperparameters was defined based on literature and preliminary results.The grid search involved training and validating models for each combination of hyperparameters in a five-fold cross-validation setup, ensuring robustness and preventing overfitting.The best hyperparameters were selected based on the highest average cross-validated score, ensuring optimal model performance (Table 3).The models were trained on the training data, and their performances were evaluated on the testing data using metrics such as the coefficient of determination (R 2 ), Root Mean Squared Error (RMSE), Root Mean Squared Error Percentage (RMSE%), Mean Absolute Error (MAE), Mean Squared Error (MSE), and Mean Absolute Percentage Error (MAPE).The predicted versus observed values and predictor importance plots were generated using the best model on the testing dataset.To identify the predictors that most significantly contribute to prediction models, we conducted a variable importance analysis.The importance of each predictor was determined by calculating its R 2 and ranking the indices from highest to lowest R 2 values.A workflow showing the main steps in modelling Bambara groundnut canopy state variables is shown in Figure 3.
The predicted versus observed values and predictor importance plots were generated using the best model on the testing dataset.To identify the predictors that most significantly contribute to prediction models, we conducted a variable importance analysis.The importance of each predictor was determined by calculating its R 2 and ranking the indices from highest to lowest R 2 values.A workflow showing the main steps in modelling Bambara groundnut canopy state variables is shown in Figure 3.

Correlation of Canopy State Variables with Vegetation Indices and Cumulative Vegetation Indices
Figure 4 shows strong correlations between LAI, CCC, AGB, and VIs at each growth stage, except during the vegetative and senescence stages.FVC exhibited strong correlations with VIs at all stages except senescence stage.
The low correlation observed during the vegetative stage can be attributed to low LAI and high background reflectance.LAI is crucial in linking VIs to agronomic measurements, particularly in the red and NIR spectral bands, which are sensitive to changes in aboveground biomass [45].The overall trend shows correlation increasing from the vegetative to the flowering stage, then declining from flowering to senescence, with peak values in LAI, CCC, FVC, and AGB observed at flowering.This trend is consistent with the findings of Tan et al. [46] who reported that VIs effectively estimate maize LAI, with the best results observed from the bell stage to the silking stage, when LAI experiences significant changes.Similarly, another study also highlighted the effectiveness of hyperspectral VIs in monitoring LAI across different growth stages of cotton [47].Similar variation in correlation between VIs and maize leaf chlorophyll content (LCC) throughout the growing season was reported by Yang et al. [48].Moreover, they reported that correlation between VIs and LCC varied significantly vertically in the upper and lower leaf layers during the early vegetative and maturity stages.

Correlation of Canopy State Variables with Vegetation Indices and Cumulative Vegetation Indices
Figure 4 shows strong correlations between LAI, CCC, AGB, and VIs at each growth stage, except during the vegetative and senescence stages.FVC exhibited strong correlations with VIs at all stages except senescence stage.The low correlation observed during the vegetative stage can be attributed to low LAI and high background reflectance.LAI is crucial in linking VIs to agronomic measurements, particularly in the red and NIR spectral bands, which are sensitive to changes in aboveground biomass [45].The overall trend shows correlation increasing from the vegetative to the flowering stage, then declining from flowering to senescence, with peak values in LAI, CCC, FVC, and AGB observed at flowering.This trend is consistent with the findings of Tan et al. [46] who reported that VIs effectively estimate maize LAI, with the best results observed from the bell stage to the silking stage, when LAI experiences significant changes.Similarly, another study also highlighted the effectiveness of hyperspectral VIs in monitoring LAI across different growth stages of cotton [47].Similar variation in correlation between VIs and maize leaf chlorophyll content (LCC) throughout the growing season was reported by Yang et al. [48].Moreover, they reported that correlation between VIs and LCC varied significantly vertically in the upper and lower leaf layers during the early vegetative and maturity stages.Our results show that FVC exhibited strong correlations with all VIs across most growth stages, confirming the effectiveness of VIs as indicators of vegetation cover.This finding aligns with a study on soybean growth dynamics, which also highlighted the potential of VIs in analysing vegetation cover and vigour [49].However, several complexities and confounding factors explain the weaker correlations between AGB and VIs compared to LAI, FVC, and CCC.Factors such as canopy structure, water content, and aboveground biomass distribution heterogeneity contribute to this weaker correlation.A recent study emphasised the importance in combining various VIs and canopy texture parameters for more accurate estimation of rice AGB.When both VIs and canopy structure parameters were combined, there was an improvement in the estimation accuracy [50].
The correlation declines as the canopy undergoes leaf withering and shedding from flowering to senescence.This results in reduced LAI, FVC, AGB, and CCC, alongside an increase in carotenoid content, thus decreasing correlation.Similar patterns were observed in wheat [51] and maize [52].Our study reveals that ∑VIs, from vegetative to senescence stages, exhibit stronger correlations with canopy state variables compared to single-stage VIs.Integrating VIs over time better captures canopy photosynthetic capacity, indicative of potential dry matter production.Models predicting maize yield based on ∑VIs (VI-SUMs) and VIs by area under the curve (VI-AUCs) demonstrated better stability and accuracy (R 2 = 60-65%) [53].This is comparable to Su et al. [54], who estimated rice yield using ∑VIs and leaf panicle abundance (R 2 = 0.73, RRMSE = 0.22 and R 2 = 0.75, RRMSE = 0.15, respectively).Additionally, our findings indicate that VIs composed of red and NIR bands exhibit stronger correlations with canopy state variables than VIs composed of NIR and green bands.This aligns with the well-known phenomenon of vegetation absorbing red bands more effectively and reflects NIR bands strongly.This observation is also supported by studies on other crops, such as wheat and soybean, which reported higher accuracy in biomass estimation using red and NIR bands [55][56][57].

Modeling the Relationship between Canopy State Variables and Cumulative Vegetation Indices
Using Machine-Learning Algorithms Table 4 and Figure 5 show performance results of predictive models in both training and testing.7 MAE 0.00 0.02 0.02 0.02 0.02 0.02 0.02 0.02 0.02 0.02 0.01 0.02 0.02 0.02 0.02 0.02 MSE 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 MAPE  For AGB, the gradient boosting CBR displayed an impressive training accuracy, but its performance in testing was moderate, suggesting potential overfitting.While our study employed CBR, other studies have also highlighted the efficacy of other gradient boosting algorithms such as the grasshopper optimisation algorithm-driven XGBoost (GOA-XGB) model which accurately predicted wheat AGB using multispectral bands and VIs.GOA-XGB outperformed other models, with RMSE of 0.226 kg m −2 and R 2 of 0.855 [58].MLR and TSR, on the other hand, showed a more consistent performance between training and testing, making them better performing models in terms of errors.Linear regression models, in general, have been foundational in many agronomy studies due to their interpretability.RFR and ABR both exhibited high training accuracies, but their testing performances were lower.Similarly, Wai et al. [59] employed Sentinel-2 Multispectral Instrument (MSI) derivatives and Shuttle Radar Topographic Mission (SRTM) Digital Elevation Model (DEM) together with field data and ML techniques to estimate AGB of evergreen and deciduous forest in Myanmar.They found that RFR and the gradient boosting algorithm provided moderate results (validation R 2 = 0.47, RMSE = 24.91 t/ha and R 2 = 0.52, RMSE = 34.72 t/ha).Their findings suggested that these models hold significant potential in estimating AGB.In a comparative study on AGB in plantation forests, Stochastic Gradient Boosting (SGB) outperformed RFR, especially when applied to a combined species dataset.This highlights the importance of appropriate model selection based on the specific dataset and ecological context [60].In our study, simpler models like PLSR, RR, and HR had varied performances, indicating they might not be the best fit for this particular dataset.RR performance can be influenced by its regularisation parameter, potentially leading to overfitting or underfitting.PLSR performance can be influenced by the number of latent variables used and the nature of relationships between predictors and the response variable.A study by Ohsowski et al. [61] highlighted these issues with RR and PLSR in estimating forest AGB.HR for instance is influenced by its epsilon parameter which determines its sensitivity to outliers in the dataset [37].These results indicate the importance of careful model tuning prior to modelling.For AGB, the gradient boosting CBR displayed an impressive training accuracy, but its performance in testing was moderate, suggesting potential overfitting.While our study employed CBR, other studies have also highlighted the efficacy of other gradient boosting algorithms such as the grasshopper optimisation algorithm-driven XGBoost (GOA-XGB) model which accurately predicted wheat AGB using multispectral bands and VIs.GOA-XGB outperformed other models, with RMSE of 0.226 kg m −2 and R 2 of 0.855 [58].MLR and TSR, on the other hand, showed a more consistent performance between training and testing, making them better performing models in terms of errors.Linear regression models, in general, have been foundational in many agronomy studies due to their interpretability.RFR and ABR both exhibited high training accuracies, but their testing performances were lower.Similarly, Wai et al. [59] employed Sentinel-2 Multispectral Instrument (MSI) derivatives and Shuttle Radar Topographic Mission (SRTM) Digital Elevation Model (DEM) together with field data and ML techniques to estimate AGB of evergreen and deciduous forest in Myanmar.They found that RFR and the gradient boosting algorithm provided moderate results (validation R 2 = 0.47, RMSE = 24.91 t/ha and R 2 = 0.52, RMSE = 34.72 t/ha).Their findings suggested that these models hold significant potential in estimating AGB.In a comparative study on AGB in plantation forests, Stochastic Gradient Boosting (SGB) outperformed RFR, especially when applied to a combined species dataset.This highlights the importance of appropriate model selection based on the specific dataset and ecological context [60].In our study, simpler models like PLSR, RR, and HR had varied performances, indicating they might not be the best fit for this particular dataset.RR performance can be influenced by its regularisation parameter, potentially For CCC, in training the ensemble models CBR and RFR, exhibited similar strong performance, however, in testing RFR achieved a higher R 2 than CBR.This was followed by ABR with a slightly lower performance in testing.Zhang et al. [62] employed several ensemble learning algorithms including RFR, gradient boosting decision tree (GBDT), extreme gradient boosting (XGBoost), light gradient boosting machine (LightGBM), and CBR to estimate chlorophyll-a content in water bodies using hyperspectral data.They reported that XGBoost outperformed the other models in estimating chlorophyll-a content (R 2 = 0.8351, RMSE = 6.6477 µg/L).The enhanced performances of ensemble models are due to their capacity to capture complex non-linear relationships, they are robust to noise, can handle multicollinearity, and maintain balance between bias and variance.Moreover, ensemble models combine predictions from multiple base models, hence providing better generalisation of unseen data [62].In our study, MLR, RR, TSR PLSR, and HR, had similar performance in predicting CCC.It is worth noting that these models displayed more stable performance when comparing training and testing.Ensemble models are often considered as "black-box" models as their decision-making processes are less interpretable than simpler models which provide a coefficient-based equation [63].Thus, it might be simpler to implement linear models to estimate groundnut canopy state variables.
In prediction of FVC, all models had similar performance (similar R 2 and error evaluation metrics) with RFR performing slightly better in testing.This is consistent with findings from a study that employed a combination of VIs derived from Sentinel-2 and UAV images together with ML algorithms for wheat and sugarcane FVC estimation in India.They reported that RFR and k-nearest neighbor (KNN) outperform support vector regression (SVR) and linear regressor.The RFR model achieved highest R 2 values ranging from 0.862 to 0.873 when different VIs are used as input features.Another study employed UAV imagery and three models including RFR, artificial neural network (ANN), and MLR to estimate maize FVC under varying irrigation levels [64].Results indicate that RFR was the most accurate, especially for different growth stages and water stress conditions, while MLR performed poorly for high FVC levels.However, the authors pointed to the need to test the developed models on the same crop and other crops in different locations across several growing seasons.
For LAI, the ensemble models CBR, RFR, and ABR had impressive training performance while testing performance was lower.This is indicative of potential overfitting.This is consistent with the findings of Martinez et al. [65] where they explored the use of satellite-derived VIs together with linear, non-linear, decision trees and RFR algorithms to model mangrove forests.Although their results indicated the ensemble RFR model to be the best performing with lowest RMSE, they highlighted the potential of overfitting.As observed previously, TSR, MLR, PLSR, HR, and RR had more stable performance comparing training and testing evaluation metrics.Although MLR performed best in testing there was not much difference.This suggest that these models generalise better on unseen data in testing.This corroborates with a study on soybean phenotyping where multimodal UAV hyperspectral, multispectral, and LiDAR, data collected at three growth stages were used as input to six algorithms including MLR, RFR, XGBoost, SVM, and back propagation (BP) for the estimation of LAI [66].Although their results indicate that XGBoost and RFR had better performances in validation (R 2 of 0.762, RMSE of 0.236 and R 2 of 0.737, RMSE of 0.277, respectively) it was not much different from MLR (R 2 of 0.649 and RMSE of 0.253).
To conclude, models like CBR, RFR, and ABR showed highest performance in training for estimation of groundnut canopy state variables, indicating their enhanced capacity to capture complex relationships and patterns in the data.However, as highlighted before, the performance was significantly lower in testing which might indicate the need for careful hyperparameter tuning and feature engineering to minimise potential overfitting.Moreover, high model complexity, insufficient training data, noisy data, and lack of crossvalidation might also contribute to overfitting.Models like MLR, TSR, PLSR, HR, and RR showed more consistent performance between training and testing, suggesting their reliability in estimating groundnut canopy state variables.These models are reliable choices which might be worthy of consideration in future modelling of agronomic variables.

Predictor Importance for Estimating Bambara Groundnut Canopy State Variables
As shown in Figure 6, for AGB, the ranking of feature importance is GNDVI (most important), CI green , NDVI, SR, and EVI2 (least important).In CCC estimation, SR leads, followed by NDVI and GNDVI, with EVI2 and CI green being least important.For FVC, CI green is most crucial, then NDVI and GNDVI, with SR and EVI2 as least important.Lastly, in LAI estimation, GNDVI is most important, followed by EVI2 and NDVI (similar importance), then CI green , and SR (least important).
In accordance with our findings, Gitelson and Merzlyak [67], highlighted that GNDVI differs from NDVI by substituting the red band with the green band, resulting in increased sensitivity to changes in chlorophyll concentration.This substitution of bands has led to enhancements in estimation.Our findings are consistent with the work of Sankaran et al. [68], who demonstrated strong correlations between GNDVI and seed yield, as well as biomass.Additionally, Macedo et al. [69] found a robust relationship between GNDVI and corn crop productivity, further emphasizing the importance of GNDVI in AGB estimation.In addition to GNDVI, other VIs have also been extensively studied for AGB estimation.For instance, Liu et al. [70] demonstrated the estimation of potato AGB based on UAV red-green-blue images using different texture features and crop height.Tao et al. [71] found significant correlations between VIs, red-edge parameters, and AGB, highlighting the potential for improved AGB estimation using a combination of these parameters.Furthermore, the integration of RS data, such as LiDAR and optical sensors, has been shown to enhance AGB estimation [72].
ful hyperparameter tuning and feature engineering to minimise potential overfitting.Moreover, high model complexity, insufficient training data, noisy data, and lack of crossvalidation might also contribute to overfitting.Models like MLR, TSR, PLSR, HR, and RR showed more consistent performance between training and testing, suggesting their reliability in estimating groundnut canopy state variables.These models are reliable choices which might be worthy of consideration in future modelling of agronomic variables.

Predictor Importance for Estimating Bambara Groundnut Canopy State Variables
As shown in Figure 6, for AGB, the ranking of feature importance is GNDVI (most important), CIgreen, NDVI, SR, and EVI2 (least important).In CCC estimation, SR leads, followed by NDVI and GNDVI, with EVI2 and CIgreen being least important.For FVC, CIgreen is most crucial, then NDVI and GNDVI, with SR and EVI2 as least important.Lastly, in LAI estimation, GNDVI is most important, followed by EVI2 and NDVI (similar importance), then CIgreen, and SR (least important).
In accordance with our findings, Gitelson and Merzlyak [67], highlighted that GNDVI differs from NDVI by substituting the red band with the green band, resulting in increased sensitivity to changes in chlorophyll concentration.This substitution of bands has led to enhancements in estimation.Our findings are consistent with the work of Sankaran et al. [68], who demonstrated strong correlations between GNDVI and seed yield, as well as biomass.Additionally, Macedo et al. [69] found a robust relationship between GNDVI and corn crop productivity, further emphasizing the importance of GNDVI in AGB estimation.In addition to GNDVI, other VIs have also been extensively studied for AGB estimation.For instance, Liu et al. [70] demonstrated the estimation of potato AGB based on UAV red-green-blue images using different texture features and crop height.Tao et al. [71] found significant correlations between VIs, red-edge parameters, and AGB, highlighting the potential for improved AGB estimation using a combination of these parameters.Furthermore, the integration of RS data, such as LiDAR and optical sensors, has been shown to enhance AGB estimation [72].
Similar to our results, SR was identified as having the best predictive performance for estimating canopy level chlorophyll content [73].Additionally, the red-edge region of the electromagnetic spectrum showed strong potential for estimating CCC [74].Similar to our results, SR was identified as having the best predictive performance for estimating canopy level chlorophyll content [73].Additionally, the red-edge region of the electromagnetic spectrum showed strong potential for estimating CCC [74].Furthermore, the red-edge chlorophyll index (CI red-edge ) was found to be a good and linear estimator of CCC in different cropping systems [75].The relevance of these indices in estimating CCC was further supported by studies that utilised in situ measurements of spectral reflectance, biophysical characteristics, and ecosystem CO 2 fluxes to estimate CCC [76].It is important to note that the effectiveness of these indices in estimating CCC was studied not only at the canopy level but also at the leaf scale.For instance, the relationship between leaf chlorophyll content and canopy reflectance was explored, indicating the potential for accurately extrapolating leaf-scale indices to the canopy level [77].
In the estimation of FVC, CI green emerges as the most important index, outperforming others in robustness and suitability, similar to the results of Thanyapraneedkul et al. [78].NDVI and GNDVI also played significant roles, showing strong linear relationships with FVC across various environmental conditions [76,77].Conversely, indices like SR and EVI2 were found to be less effective in predicting Bambara groundnut FVC.Similar studies by Viña et al. [79] and Luscier et al. [80] indicated that these indices had lower predictive performance and are recommended less frequently for monitoring vegetation biophysical characteristics.This suggests that the selection of appropriate VIs is critical for accurate and reliable FVC estimation.
Our results revealed that GNDVI was the best index for estimating Bambara groundnut LAI.The importance of GNDVI for LAI estimation is supported by Dai et al. [81], who found that GNDVI, along with other indices such as EVI and Soil Adjusted Vegetation Index (SAVI), exhibited higher correlations with LAI compared to NDVI.EVI2 and NDVI have also been highlighted as important indices for LAI estimation.Fang et al. [82] discussed the widespread use of NDVI for LAI estimation, indicating its significance in this context.Kang et al. [83] emphasised the comparable performance of EVI2 and EVI in LAI estimation at global scales, further underlining the importance of these indices.On the other hand, our results indicate CI green and SR were identified as the least important indices for Bambara groundnut LAI estimation.This is supported by the study of Stenberg et al. [84] who found that CI green and SR exhibited poor sensitivity to changes in LAI due to saturation.In conclusion, the importance of VIs for LAI estimation varies, with GNDVI being the most important, followed by EVI2 and NDVI, with CI green and SR being the least important.

Limitations and Future Perspectives
A limitation of our study is the restricted temporal scope, focusing on data from 2018-2019.Key crop growth information is influenced by multiple factors such as environment and temperature, necessitating a longer period of data collection to capture inter-annual variability.Moreover, expanding the study to cover more regions and a wider range of genotypes is crucial for comparative analysis and generalisation of the findings.Such comprehensive data would allow for a more robust understanding of crop growth dynamics and the development of models that are resilient across different environmental conditions and genetic variations.
Another significant challenge in our study is the underlying issue of overfitting.To mitigate overfitting in ML modelling, several techniques have been proposed.One of the main methods is k-fold cross-validation, which assesses model performance on different subsets of the data to ensure generalisation [85].Feature selection is another technique whereby irrelevant features are eliminated to reduce model complexity [86].Regularisation techniques, such as L1 (Lasso) or L2 (Ridge), which penalise large coefficients aiming to simplify models are also effective.Novel methods like L1/4 regularisation can address both overfitting and underfitting [87].Another technique is monitoring model performance on a validation set and employing early stopping when performance deteriorates.Data augmentation techniques, for example synthetic data augmentation for tabular data (SMOTE), generative adversarial networks (GAN), and combined use of SMOTE and GAN, can effectively increase training dataset size thereby increasing model variations and handling concerns of overfitting and loss of generalisation [88].Moreover, ensemble models, such as RFR and gradient boosting regressor, have been shown to offer high resilience against overfitting.However, opting for simpler models or reducing the complexity of the models are beneficial, as observed in our results.Finally, fine-tuning hyperparameters can effectively minimise both underfitting and overfitting.
Future studies should explore the use of Explainable artificial intelligence to understand the decision-making mechanisms models and interactions among various indices and algorithms.The integration of data from various in situ and other sources, such as meteorological data, field ancillary data, and real-time data, has the potential to enhance understanding and improve prediction accuracy.In addition, the utilisation of advanced deep learning methodologies, such as convolutional neural networks, has promise in the extraction of features from RS data, thereby improving prediction accuracy.Moreover, multispectral and hyperspectral data, in conjunction with the abovementioned data, have the potential to provide a comprehensive understanding of crop growth and development.Future studies should explore the integration of radiative transfer models with deep learning for gauging the estimation of crop state variables.Such hybrid methods combine both the strengths of process-based and data-driven models for scalable, accurate, and robust prediction models.

Conclusions
The study investigated correlations between LAI, AGB, FVC, and CCC with VIs/∑VIs obtained at single stages and over combinations of stages.The highest single-stage correlation was observed at flowering, while ∑VIs spanning from vegetative to senescence showed the strongest correlations with all canopy state variables.Although ensemble models performed well in training across all canopy state variables, they exhibited signs of potential overfitting during testing.In contrast, simpler models consistently provided reliable results in both training and testing.Our results show that MLR is particularly effective for AGB and LAI estimation, while RFR shows superior performance for CCC and FVC estimation.Our findings highlight GNDVI as the most crucial index for estimating Bambara groundnut AGB and LAI, while SR proves optimal for CCC and CI green for FVC, underlining the importance of careful VI selection in canopy state variable estimation.

Horticulturae 2024 , 20 Figure 1 .
Figure 1.Location of the study area at Field Research Centre of Crops for the Future.The experimental layout of plots was digitised on an image acquired with the integrated DJI Phantom 4 Pro camera at a height of 10 m on flowering stage.B1G1R1 means plot is in block 1; genotype is genotype 1, and replicate is the first replicate.

Figure 1 .
Figure 1.Location of the study area at Field Research Centre of Crops for the Future.The experimental layout of plots was digitised on an image acquired with the integrated DJI Phantom 4 Pro camera at a height of 10 m on flowering stage.B1G1R1 means plot is in block 1; genotype is genotype 1, and replicate is the first replicate.

Figure 1 .
Figure 1.Location of the study area at Field Research Centre of Crops for the Future.The experimental layout of plots was digitised on an image acquired with the integrated DJI Phantom 4 Pro camera at a height of 10 m on flowering stage.B1G1R1 means plot is in block 1; genotype is genotype 1, and replicate is the first replicate.

Figure 2 .
Figure 2.This depicts the environmental parameters and irrigation levels throughout the 2018 cultivation period spanning from May to September.The arrows indicate distinct growth phases in the life cycle of Bambara groundnut.These stages include SOW (sowing), VEG (vegetative), FLO

Figure 2 .
Figure 2.This depicts the environmental parameters and irrigation levels throughout the 2018 cultivation period spanning from May to September.The arrows indicate distinct growth phases in the life cycle of Bambara groundnut.These stages include SOW (sowing), VEG (vegetative), FLO (flowering), POD (podding), PF (pod filling), MAT (maturity), SEN (senescence), and HAR (harvest).The asterisk (*) denotes data collection time points.

Figure 3 .
Figure 3. Workflow for modelling Bambara groundnut canopy state variables.Figure 3. Workflow for modelling Bambara groundnut canopy state variables.

Figure 3 .
Figure 3. Workflow for modelling Bambara groundnut canopy state variables.Figure 3. Workflow for modelling Bambara groundnut canopy state variables.

Figure 5 .
Figure 5. Plot comparing predicted versus observed values using the top-performing models for each canopy state variable.The solid black line represents the best-fit line, while the dashed grey line corresponds to the line y = x.

Figure 5 .
Figure 5. Plot comparing predicted versus observed values using the top-performing models for each canopy state variable.The solid black line represents the best-fit line, while the dashed grey line corresponds to the line y = x.

Figure 6 .
Figure 6.Predictor importance plots ranking ΣVIs for estimating Bambara groundnut canopy state variables, where higher feature importance values indicate greater importance in the model.

Figure 6 .
Figure 6.Predictor importance plots ranking ΣVIs for estimating Bambara groundnut canopy state variables, where higher feature importance values indicate greater importance in the model.

Table 1 .
Summary statistics for Bambara groundnut canopy state variables.Min denotes the minimum value; Mean is the average value; Max indicates the maximum value; SD represents the standard deviation in the data; and CV (%) signifies the percent coefficient of variation.The sample size (N) for LAI, CCC, and FVC was 216 while N for AGB was 144.

Table 2 .
Vegetation indices (VIs) used in this study.

Table 3 .
Hyperparameters evaluated for optimising the ML models.