Forecasting Table Beet Root Yield Using Spectral and Textural Features from Hyperspectral UAS Imagery

: New York state is among the largest producers of table beets in the United States, which, by extension, has placed a new focus on precision crop management. For example, an operational unmanned aerial system (UAS)-based yield forecasting tool could prove helpful for the efﬁcient management and harvest scheduling of crops for factory feedstock. The objective of this study was to evaluate the feasibility of predicting the weight of table beet roots from spectral and textural features, obtained from hyperspectral images collected via UAS. We identiﬁed speciﬁc wavelengths with signiﬁcant predictive ability, e.g., we down-select >200 wavelengths to those spectral indices sensitive to root yield (weight per unit length). Multivariate linear regression was used, and the accuracy and precision were evaluated at different growth stages throughout the season to evaluate temporal plasticity. Models at each growth stage exhibited similar results (albeit with different wavelength indices), with the LOOCV (leave-one-out cross-validation) R 2 ranging from 0.85 to 0.90 and RMSE of 10.81–12.93% for the best-performing models in each growth stage. Among visible and NIR spectral regions, the 760–920 nm-wavelength region contained the most wavelength indices highly correlated with table beet root yield. We recommend future studies to further test our proposed wavelength indices on data collected from different geographic locations and seasons to validate our results.


Introduction
Table beet (Beta vulgaris spp. vulgaris: Family Chenopodiaceae) consumption has increased recently, primarily due to an enhanced awareness of the potential health benefits [1], e.g., beets are a good source of dietary fiber and potassium [2]. Furthermore, Betalains from table beet roots represent a new class of dietary cationized antioxidants with excellent antiradical impact and antioxidant activity, and they have been linked to cancer prevention [3,4]. With the ever-increasing popularity and increasing production of table beets, our ability to predict root yield before harvest has become essential for both logistical planning and within-season crop inputs to manipulate yield.
Researchers have identified the suitability of unmanned aerial systems (UAS) to spectrally map large crop areas in a short duration, specifically for applications requiring high spatial resolution. The increased availability of UAS and miniaturized sensors [5] has also catalyzed various applications for UAS in agriculture [6]. Uses include, but are not limited to, the assessment of nutrient content [7,8], disease [9,10], above-ground biomass [11], leaf area index [12], evapotranspiration [13], water stress [14], weed presence [15], and yield [16][17][18] of various crops.
Hyperspectral sensors sample a wide range of narrow, contiguously spaced spectral information. As such, hyperspectral imagery (HSI) contains more detailed spectral In this study, our objective was to assess the predictability of beet root yield (weight per unit length) entirely from HSI to determine the benefits of high-dimensional spectral data. A secondary objective was to quantify the effect of using both narrow-band spectral and textural features individually and in combination. Finally, we evaluate table beet root yield at five different growth stages during the season to assess the ideal flight timing and its impact on model performance.

Data Collection
Table beets (cv. Ruby Queen) were planted at Cornell AgriTech, Geneva, New York, USA. The beets were planted in a 60 m × 100 m field on 20 May 2021, representing 18 different plots where the root weight of the table beets was collected. Each of the plots was 1.52 m (5 ft) long. The roots were manually harvested on August 5, 2021, after which root weight was measured.
We collected hyperspectral imagery of the entire field using a Headwall Photonics Nano-Hyperspec imaging spectrometer (272 spectral bands; 400-1000 nm), which was fitted onto DJI Matrice 600 UAS. The ground sampling distance of our imagery was 3 cm. A hyperspectral image cube, captured from the UAS for all 18 plots, is shown in Figure 1. A MicaSense RedEdge-M camera [54] also was fitted to the drone, which captured five bands (475, 560, 668, 717 and 840 nm) of multispectral images. We flew our UAS at five different stages during the growing season, from the sowing date to harvest. The synopsis of the data collection is shown in Table 1. The first flight date was selected to align with canopy emergence. The second flight was chosen to correspond with canopy closure (last phase of the leaf development). The rest of the dates were chosen based on weather conditions. s. 2022, 14, x FOR PEER REVIEW

Data Preprocessing
The raw digital count of the HSI was orthorectified and converted to r

Data Preprocessing
The raw digital count of the HSI was orthorectified and converted to radiance using Headwall's Hyperspec III SpectralView V3.1 software [55]. Simultaneously, an orthorectified, spatially registered mosaic of the multispectral image was formed using Pix4D V4.6.3 [56]. Ten georegistered ground control points (AeroPoints) were deployed in the field to ensure accurate spatial registration of the multispectral images. These multispectral mosaics were then used to register the HSI spatially. All the study plots were marked in the field with orange plates, staked above the table beet canopy, which helped visually identify the plots from the HSI mosaics. The study plots were clipped from the mosaic via a rectangular per-plot shapefile. The shapefiles' dimensions were 0.67 m (distance between two adjacent crop rows) across the row and 1.52 m (the length of the study plot) along the row. The box and HSI were spatially registered to ensure image analysis of the same crops across all the flights.
Since radiance data are susceptible to changes in illumination conditions, the HSI data were converted from radiance to reflectance [19]. Four panels (black, dark gray, medium gray, and light gray) were placed in the field to serve as calibration targets. The reflectance of each panel was measured using an SVC spectrometer. We also extracted the radiance of each panel from the HSI mosaic of the UAS data. The empirical line method tool in ENVI V5.3 [57] was used to convert the radiance HSI of the study plots to reflectance. We ignored the last 32 wavelength bands in the HSI imagery, given their noisy nature toward the detector fall-off range.

Denoising
Most UAS-based HSI data contain thermal noise, quantization noise and shot noise [58]. We suppressed as much noise from the image as possible to draw conclusive results before performing any analysis.
Our high-level denoising procedure was similar to Chen et al. [59]. This method hinges on the fact that most noise is contained in principal components with lower eigenvalues. As a result, a principal component analysis (PCA) was performed on the scene containing n bands, and the top-k components were retained and kept unaltered, with the remaining components being passed for denoising; k was determined by the minimum average partial (MAP) test [60]. Consecutively, n-k components were passed to dual-tree complex wavelet transformations [61] for denoising, after which an inverse PCA was performed to obtain the denoised image. We refer the reader to [18] for more information regarding this approach.

Feature Extraction
The high dimensionality of HSI is a well-documented problem across scientific literature [20]. Therefore, we carefully selected the bands, or a combination of bands, to derive a predictive model to reduce the complexity of analysis when dealing with a high number of bands.  [62] was used to mask the vegetation from the soil in the imagery, after which the mean vegetation spectra for each plot were evaluated. Figure 2a shows the mean vegetation spectra across all flights. We observed the increase in reflectance in the NIR region with each successive flight, which was expected since as vegetation matures, canopy layering typically increases. [69]. The shortcomings of using NDVI have also been well documented [70][71][72][73][74]. NDVI is restricted to just two specific band pairs, so instead, we evaluated normalized difference ratio indices (NDRI) of every possible combination of narrow-band reflectance, as follows: This method of identifying wavelength pairs has been used in [75] to understand the relationship between evapotranspiration, to evaluate biophysical characteristics of cotton, potato, soybeans, corn and sunflower crops [76], and to estimate biomass and nitrogen concentrations in winter wheat [77].

Texture Features
The texture of a plot is another possible source of features for our model. Haralick et al. [25] documented the most widely used texture feature extraction method, which was also used in this study. It involves the calculation of the gray-level co-occurrence matrix (GLCM), followed by extracting descriptive statistics from the matrix. GLCM is a Huete et al. [63] examined cotton canopy spectra and concluded that the spectra of soil and vegetation mix interactively, and thus, normalization is crucial to remove the background effect. Normalized difference indices of various wavelength bands have been found to normalize soil background spectral variations [64] while simultaneously reducing the effect of sun and sensor view angle [65]. Additionally, the well-known normalized difference vegetation index (NDVI) has historically been used to map vegetation [66] while also being used to predict leaf area index (LAI) [67], crop biomass [68], and yield [69]. The shortcomings of using NDVI have also been well documented [70][71][72][73][74]. NDVI is restricted to just two specific band pairs, so instead, we evaluated normalized difference ratio indices (NDRI) of every possible combination of narrow-band reflectance, as follows: This method of identifying wavelength pairs has been used in [75] to understand the relationship between evapotranspiration, to evaluate biophysical characteristics of cotton, potato, soybeans, corn and sunflower crops [76], and to estimate biomass and nitrogen concentrations in winter wheat [77].

Texture Features
The texture of a plot is another possible source of features for our model. Haralick et al. [25] documented the most widely used texture feature extraction method, which was also used in this study. It involves the calculation of the gray-level co-occurrence matrix (GLCM), followed by extracting descriptive statistics from the matrix. GLCM is a frequency table representing the number of times a specific tone occurs next to one another. The tone is the quantized level of pixel values for a particular spectral band. One of the many descriptive statistics is the mean of GLCM, and it is given by Equation (2), where N g is the quantization level and P(i, j) is the probability of occurrence of pixel j next to i.
In several previous studies [11,35,37,38,78], the top-performing models predominantly consisted of the mean of GLCM. Furthermore, the GLCM mean represents the value of the more frequently occurring mean tonal level within a particular window size. Thus, the GLCM means contains both spatial and spectral variation information.
The texture features were calculated over all 240 narrow wavelength bands. An 11 × 11 kernel size was used to calculate the GLCM mean for each kernel over each study plot. Both small [11] and large kernel sizes [35] have yielded promising results for their specific scenarios. For our case, using a smaller kernel size meant that the feature became more susceptible to noise, while too large a value could lead to over-smoothing of the data and loss of texture information. The GLCM was evaluated in eight directions (N, S, E, W, N-E, N-W, S-E, S-W) for each kernel. A quantization level of 16 was used to evaluate the GLCM mean for each kernel within a plot, after which the mean of GLCM means was taken as the texture feature for each plot. The texture features of an image depend on the canopy coverage of a plot [11,79], so we did not crop out the vegetation during the calculation. Figure 2b is a plot of the mean of GLCM means over a particular plot at different flight times. We note the dip in spectra at around 680 nm due to the contrast in reflectance of vegetation and non-vegetation pixels at that wavelength band. Note that the contrast at 480 nm was not high enough; thus, there was no significant dip, even though there were absorption spectra in reflectance due to chlorophyll absorption. Vegetation shadow pixels result in variation of reflectance within the green and the NIR region, leading to a lower GLCM mean. We note that the relative uniformity of GLCM mean values for all the wavelength regions for flight 1, compared to other flights, was due to a greater number of soil pixels in the image.
In a study by Zheng et al. [11], textural measurements from individual wavelengths did not exhibit a high predictive capability of rice AGB. The authors therefore proposed a normalized difference texture index (NDTI). The texture feature extracted in this paper is also the NDTI, as defined by Equation (3). Here, T 1 and T 2 are the means of GLCM means at two different wavelengths. We also evaluated NDTI with each possible combination of narrow-band texture wavelengths, as follows:

Model Evaluation and Feature Search
Our goal was to identify the pair (or pairs) of wavelength bands that produced the most predictive NDRI and NDTI features, for an accurate and precise model for table beet root yield. Due to low sample size, we performed leave-one-out cross-validation (LOOCV) to evaluate our model [80]. Here, the dataset is divided into a number of subsets, which equals the number of instances in that dataset. In each subset, all but one instance was used to train models, while the omitted sample was used to test models. In our case, a model with the same feature was trained on 17 occasions and tested on one occasion with 18 combinations. Our model evaluation parameters were the coefficient of determination R 2 and root mean square error (RMSE). The LOOCV R 2 and RMSE for a model were calculated from these 18 unique occasion predictions. We investigated the model performance for each flight individually, and the final models for each flight were evaluated in multiple ways. Firstly, we noted the top 10 LOOCV R 2 features obtained from NDRI and NDTI. We also noted the models obtained using a combination of one NDRI and one NDTI feature. Finally, we fed the top-10 NDRI and Remote Sens. 2023, 15, 794 7 of 21 NDTI features into a random forward-selected stepwise linear regression to obtain models that contained multiple features.

Single Feature Approach
To identify the wavelength bands that were best predictors of table beet root yield, the NDRI and NDTI for each possible wavelength pair in the spectrum were fitted with root yield in a simple linear regression (SLR). The predictive power from single features was identified for two purposes: to narrow down the feature search space for a general model and to evaluate the performance of a model that predicts root yield from just two wavelength bands. Additionally, we also wanted to monitor the change in predictive capabilities at each growth stage.

Double Feature Approach
We selected one feature from NDRI and one from NDTI and evaluated the model performances. Each model formed from all possible dual feature combinations of the top 10 NDRI and NDTI features was evaluated. Both feature extraction methods were used, as each signifies different properties of an image. NDRI is a placeholder for spectral characteristics, while NDTI provides information related to image texture. We selected one index from each, as the goal of the final model was to select just four wavelength bands, to ensure compatibility with inexpensive sensors.
Zheng et al. [11] showed that using spectral and textural measurements improved their rice AGB prediction model. Furthermore, Zheng et al. [78] showed significant improvements in the estimation of leaf nitrogen concentration when using VIs combined with NDTIs, compared to the indices alone. Yue et al. [79] also reported improved results when using vegetation indices in combination with image textures.

Modified Stepwise Regression
Finally, to find the best model derived from NDTI and NDRI features during each flight, stepwise regression offers a simple way to explore regression models over a large feature space without exhaustively interrogating every possible combination. In order to mitigate some of the concerns raised in the literature [81][82][83][84], we introduced some additional procedures in the stepwise regression algorithm.
Smith [81] argues the probability of choosing an explanatory variable decreases if the number of candidate variables is too large. Therefore, we restricted our procedure to pool 20 candidate variables. Grossman et al. [82] highlighted that the stepwise procedure features depend on the samples in the dataset, so LOOCV p-values were chosen to mitigate this issue. The LOOCV p-values are the mean of the p-values of a feature for each instance within LOOCV. Simple forward selection of features also does not account for interdependencies within the features [85], which can lead to the entry of redundant features into the final model, thus inflating its performance [86]. Therefore, we performed backward elimination after adding each new feature to discard all features with a p-value greater than 0.05.
In this process, the order of parameter entry affects the selection of features [83], while the stepwise regression algorithm only accounts for the interaction of a particular set of features at a fixed permutation [84]. The first feature for the stepwise regression algorithm was thus randomly selected, and each new feature was randomly added to the algorithm. The algorithm was repeated multiple times, ensuring that a broad range of sub-models was evaluated. Finally, multiple sets of features were obtained that could lead to predicting table beet root yield. Among the multiple models, we selected those exhibiting an R 2 adj value greater than 0.75 and the VIF for each feature as less than five. VIF is the variance inflation factor, a measure of the degree of inter-correlation of one feature with the rest of the features in a model [87]. A VIF greater than five shows high inter-correlation [88]. Models with features exhibiting the above criteria have the potential to be operational.

Results
In this section we report all the LOOCV R 2 , RMSE and relevant models' performance parameters for each approach.

Single Features
The R 2 values for each pair of narrow band wavelengths of NDRI and NDTI were recorded and shown in    Table A1 in Appendix A. The top 10 R 2 values for the first two flights were between 0.49-0.61. Flights 3 and 4 had more predictive features, ranging from 0.63-0.73, while the top-10 R 2 values of Flight 5 were all around 0.6. The most predictive features for all flights appear primarily in the NIR region (750-900 nm).
The predictive regions for NDTI are sparsely distributed, as seen in Figure 3b. For flight 1, the predictive features appear in the NIR pairs, with R 2 values of around 0.4 and 0.5. With a maximum R 2 of 0.36, the NDTI features of flight 2 perform poorly. There are some relatively high R 2 regions around the 500-700 nm wavelength pairs for flights 3 and 4. Moreover, the features from flight 5 have a low predictive ability with table beet root yield, with the top R 2 ranging in the range 0.25-0.30. Table A2 in Appendix A lists the top 10 NDTI features.

Double Features
The top five models (in terms of R 2 values) obtained from combinations of each feature are reported in Table 2. Incorporating both NDRI and NDTI features increased the R 2 values of the model for all five flights. We also noted that each model's VIF values are low, signifying the absence of correlation between the features. The dual features had better-performing models for flights 3 and 4 than the other flights. While the wavelength bands were different, there were some frequently occurring features; for example, for flight 4, the NDRI features were either 895 nm and 785.9 nm or 895 nm and 783.7 nm, while the NDTI features were around 850 nm and 760 nm. For flight 1, both NDRI and NDTI features were from wavelength pairs in the NIR, except for one pair of blue wavelengths. For flight 2, the NDRI features are from NIR pairs, but the NDTI features were at different wavelength pairs. For flight 3, we again observed the high R 2 for blue wavelength pairs. Finally, for flight 5, we observed some texture features from blue and red pairs.

Multiple Features
Random stepwise regression with an alpha value of 0.05 was implemented. The summary of obtained models are shown in Table 3. Most of the models identified had an R 2 value of greater than 0.80 and an RMSE of lower than 15.27% of the estimate, which is equivalent to around 0.28 kg/m. The top-performing models for each flight had similar R 2 values, ranging from 0.85 to 0.90, with the best performing model occurring at flight 4, with an R 2 of 0.90 and 10.81% (0.2 kg/m) RMSE. Flights 3 and 5, and 1 and 2 had similar performances compared to each other, albeit each with their own unique set of predictive features. Table 3. Summary of results using multiple features. Random stepwise regression also generated models already found in the double feature case we did not report those models in this We identified select operational models as defined by our previously mentioned performance criteria (R 2 adj > 0.75 and VIF < 5) for all flights. However, we note the scarcity of models for flight 2 that meet this performance criteria. We emphasize that for the conditions we set to filter our models, flights 1 and 2 contained four features, while the rest of the top-performing models were obtained from just three features. For flight 1, most of the models are focused in the NIR region. For the fifth flight, however, there was a blue texture feature, while for the fourth flight, a red feature was present. We observed red and blue NDRI pairs, with NIR texture feature pairs, for the third flight. The fit of the best-performing models for each flight is shown in Figure 4. formance criteria (R adj >0.75 and VIF < 5) for all flights. However, we note the scarcity of models for flight 2 that meet this performance criteria. We emphasize that for the conditions we set to filter our models, flights 1 and 2 contained four features, while the rest of the top-performing models were obtained from just three features. For flight 1, most of the models are focused in the NIR region. For the fifth flight, however, there was a blue texture feature, while for the fourth flight, a red feature was present. We observed red and blue NDRI pairs, with NIR texture feature pairs, for the third flight. The fit of the bestperforming models for each flight is shown in Figure 4.

Extrapolating Yield Modeling Results to the Field Scale
Yield maps are helpful tools for farmers to plan and make logistic decisions, and here we apply our model to provide insight into one such yield map. Figure 5 shows the predicted table beet root yield for plots in the northern region of the field. Each prediction was made by training the best-performing model on all 18 plots. In the plot, some boxes were empty, since the predicted root yields from the model were outside the labeled range. This occurred as some of the images within the plot appeared blurred in the mosaic, leading to an anomalous prediction. We observed similar regions or locally specific results, based on similar colors around the same plot. However, there was also some variability in growth factors across the field and across the growing season, which is understandable, since soil variability came into play, and the time-specific model performances, respectively, were different. The projection of the mosaic was not always uniform, and there were also missing crop locations in the plot. This resulted in some of the plots having anomalous results, which are shown as empty blocks in the image. We notice that, specifically for flights 3 and 4, the predicted table beet root yield is (spatially) similarly distributed.

Significance of Obtained Features
An imaging system's limitations must be considered before analyzing the prominent wavelength bands and the spectra for any given application. A number of critical factors include the band centers and the FWHM of the sensors in the imager, as well as flight shifts in airborne hyperspectral imaging systems due to vibrations [89]. These effects are further exacerbated by SMILE and keystone effects of the sensor [90], which lead to spectral absorption features appearing at slightly shifted locations. Furthermore, our sensor has a FWHM of around 6 nm, which is relatively broad compared to our spectral resolution of 2.2 nm [55]. These factors led to our observed absorption peaks appearing as relatively broad spikes, as seen in Figure 2(a). We mention this as background to the discussion of the specific wavelength features and pairs, especially insofar as their specificity of location and combination is concerned.
The wavelength bands of each feature for the best-performing model (in terms of R 2 ) of each flight are shown in Figure 6. We marked oxygen absorption at 760 nm and water vapor absorption at 820 nm and 910 nm [19] in the plot. Most of the spectral/texture features for all flights were located around the 900 and 800 nm pairs, attributable to the broad spectral spikes caused by oxygen and water absorption features. This has potential physiological implications. Table beet roots primarily consist of water [91]. The amount of water in the foliage could be an indication of the relative hydration of the roots and thus indicative of root weight. We contend that this is a type of normalized water absorption feature. Finally, by analyzing observed vs. predicted plots in Figure 4, we observe a higher error variability for our models, especially for flights 1, 2, and 5. Thus, we do not see a similar trend of predicted values that we observe for flights 3 and 4. As we do not possess the actual table beet root yield of the plot; i.e., it is difficult to verify the entire map's accuracy, but we assume at least some validity in the table beet root yield map, since we hedged against the overall performance of flight 3 (growth stage).

Significance of Obtained Features
An imaging system's limitations must be considered before analyzing the prominent wavelength bands and the spectra for any given application. A number of critical factors include the band centers and the FWHM of the sensors in the imager, as well as flight shifts in airborne hyperspectral imaging systems due to vibrations [89]. These effects are further exacerbated by SMILE and keystone effects of the sensor [90], which lead to spectral absorption features appearing at slightly shifted locations. Furthermore, our sensor has a FWHM of around 6 nm, which is relatively broad compared to our spectral resolution of 2.2 nm [55]. These factors led to our observed absorption peaks appearing as relatively broad spikes, as seen in Figure 2a. We mention this as background to the discussion of the specific wavelength features and pairs, especially insofar as their specificity of location and combination is concerned.
The wavelength bands of each feature for the best-performing model (in terms of R 2 ) of each flight are shown in Figure 6. We marked oxygen absorption at 760 nm and water vapor absorption at 820 nm and 910 nm [19] in the plot. Most of the spectral/texture features for all flights were located around the 900 and 800 nm pairs, attributable to the broad spectral spikes caused by oxygen and water absorption features. This has potential physiological implications. Table beet roots primarily consist of water [91]. The amount of water in the foliage could be an indication of the relative hydration of the roots and thus indicative of root weight. We contend that this is a type of normalized water absorption feature. There also are feature pairs around 650 nm and 500 nm, even if only one feature appears around this region for the best-performing models for NDRI. However, from the heat map in Figure 3, we noted the significance of this region for predicting table beet root yield. Strong chlorophyll a and b absorptions can be found at around 440 nm and 680 nm, and 460 and 635 nm, respectively, while carotenoids have been shown to express as an absorption feature at around 470 nm [19,92].

ER
Additionally, similar wavelength regions have been shown to have utility for the development of a broad range of crop physiological models. Thenkabail et al. [76] demonstrated the importance of 800 nm and 900 nm wavelength pairs, which significantly correlated with the wet biomass for soybean and cotton crops. The 672 nm and 448 nm combination also has been shown to be an indicator of evapotranspiration by Marshall et al. [75]. Finally, 900-700 nm and 750-400 nm wavelength pairs are predictive of AGB and LAI for paddy rice [93].

Model Performance
Several different models were developed to identify possible wavelength combinations that could lead to effective scalable solutions for the prediction of table beet root yield at each growth stage. However, we are aware of the relatively small sample size on which we performed our analysis. To mitigate the small number of samples, we evaluated our models by LOOCV. All our reported models arguably rely on a limited number of predictive features, thereby avoiding the challenge of over-fitting in the context of a limited sample set. Hyperspectral indices are highly susceptible to noise [58], so the top-reported wavelength indices may not be extrapolated to other locations and times. However, a broad range of models with LOOCV R 2 greater than 0.80 and RMSE of less than There also are feature pairs around 650 nm and 500 nm, even if only one feature appears around this region for the best-performing models for NDRI. However, from the heat map in Figure 3, we noted the significance of this region for predicting table beet root yield. Strong chlorophyll a and b absorptions can be found at around 440 nm and 680 nm, and 460 and 635 nm, respectively, while carotenoids have been shown to express as an absorption feature at around 470 nm [19,92].
Additionally, similar wavelength regions have been shown to have utility for the development of a broad range of crop physiological models. Thenkabail et al. [76] demonstrated the importance of 800 nm and 900 nm wavelength pairs, which significantly correlated with the wet biomass for soybean and cotton crops. The 672 nm and 448 nm combination also has been shown to be an indicator of evapotranspiration by Marshall et al. [75]. Finally, 900-700 nm and 750-400 nm wavelength pairs are predictive of AGB and LAI for paddy rice [93].

Model Performance
Several different models were developed to identify possible wavelength combinations that could lead to effective scalable solutions for the prediction of table beet root yield at each growth stage. However, we are aware of the relatively small sample size on which we performed our analysis. To mitigate the small number of samples, we evaluated our models by LOOCV. All our reported models arguably rely on a limited number of predictive features, thereby avoiding the challenge of over-fitting in the context of a limited sample set. Hyperspectral indices are highly susceptible to noise [58], so the top-reported wavelength indices may not be extrapolated to other locations and times. However, a broad range of models with LOOCV R 2 greater than 0.80 and RMSE of less than 15% for table beet root crops have been documented in this study.
In our previous study [53], the best performing model had an R 2 of 0.89 and RMSE of 2.5 kg at the canopy closing stage using VDVI (Visible Band Difference Vegetation Index) and canopy area obtained from multispectral imagery. Here, we obtained similar and sometimes better results for each growth stage. The R 2 for predicting carrot yield by Wei et al. [94] using satellite imagery was around 0.80. The best performing model for predicting potato yield for Li et al. [49] had an R 2 of 0.81, using spectral indices collected from UAS-based HSI. Therefore, our results outperform those analogous studies. However, we are aware of the sample size limitation of our study, as we only used 18 plots to establish our model. We therefore recommend that our models be verified with table beets grown at different geographic locations and across years, with a larger dataset encompassing multiple fields.
Although the top-performing models for each flight exhibit similar performances, the best performing model was from flight 4. In contrast, both flights 3 and 4 have the largest number of high performing models (12 and eight models with R 2 value greater than 0.80, respectively). Flights 3 and 4 were performed within five days of each other, and as such, they exhibit similar feature sets, as well as model performances.
All the feature sets from each flight (as reported in Tables 2 and 3) were assessed on the remaining four flights' data, and their LOOCV R 2 values are shown in Figure 7.
Only the top five model performances across all flights are shown. From these plots, the top-performing feature for one flight does not necessarily perform well on other flight data, meaning the features' foresight depended on the day the data was captured. Although the R 2 values for flight 3 is relatively high with feature sets from all flights. This signifies the presence of abundant and broader wavelength features that potentially exhibit table beet root yield predictive ability during this period. number of high performing models (12 and eight models with R 2 value greater than 0.80, respectively). Flights 3 and 4 were performed within five days of each other, and as such, they exhibit similar feature sets, as well as model performances.
All the feature sets from each flight (as reported in Table 2 and Table 3) were assessed on the remaining four flights' data, and their LOOCV R 2 values are shown in Figure 7.
Only the top five model performances across all flights are shown. From these plots, the top-performing feature for one flight does not necessarily perform well on other flight data, meaning the features' foresight depended on the day the data was captured. Although the R 2 values for flight 3 is relatively high with feature sets from all flights. This signifies the presence of abundant and broader wavelength features that potentially exhibit table beet root yield predictive ability during this period.
Moreover, the features obtained from flight 5 show relatively high coefficient of determinations for all flights, with R 2 values ranging from 0.4-0.5 for flights 3 and 4, and 0.1-0.2 for flights 1 and 2. The features are close to the best-performing features of flights 3 and 4, leading to relatively higher R 2 scores. In addition, some of the features were similar to the features of flights 1 and 2. However, their R 2 values are lower, signifying that during the earlier stages of growth, the predictive ability was more in tune with the specific narrow-band indices, while as the crop matures and approaches harvest, a broader range of indices shows prognostic capability.

Conclusions
The goals of this study were to assess the utility of HSI to accurately/precisely model beet root yield and to identify the wavelength indices that exhibit a predictive capability The features are close to the best-performing features of flights 3 and 4, leading to relatively higher R 2 scores. In addition, some of the features were similar to the features of flights 1 and 2. However, their R 2 values are lower, signifying that during the earlier stages of growth, the predictive ability was more in tune with the specific narrow-band indices, while as the crop matures and approaches harvest, a broader range of indices shows prognostic capability.

Conclusions
The goals of this study were to assess the utility of HSI to accurately/precisely model beet root yield and to identify the wavelength indices that exhibit a predictive capability for beet root weight. We developed a varied selection of wavelength indices at various growth stages, which farmers could utilize to design a system for predicting table beet root yield with parameters best suited to their means.
The potency of predicting root yield from NIR pairs, especially for spectral features, is highlighted in this study. These results are promising for developing cost-effective and operational sensing platforms, as most of the spectral/texture features were located in the spectral range that can be detected by relatively cheap silicon (Si) detector material. Furthermore, the addition of texture features showed significant improvement to our model, showing that HSI texture features could serve as an important tool for analysis. We explored its effectiveness with some specific hyperparameters for texture feature extraction. However, variation of hyperparameters leads to varying levels of performance, which are also dependent on the growth stage of a plant [35,78]. Therefore, this is a potential area for future research, i.e., to find the hyperparameters suitable for a particular study at a particular growth stage.
The flight that was performed 21 days prior to harvest (or equivalently, 56 days after planting) exhibited a more significant number of features that led to good-performing (LOOCV R 2 > 0.80) models. Additionally, we noted that the features identified from other flights showed significant fit (R 2 around 0.4) with data obtained from this flight. While flight 4 (which was performed five days after) had the best-performing model, with an R 2 of 0.90 and a 10.81% RMSE. We thus surmise, with caution, that the time period of roughly 16 to 21 days from harvest (Rosette growth stage) is ideally the best time to fly for predicting table beet root yield. However, we also acknowledge the results to depend on the weather and growing condition of the 2021 harvest at Geneva, New York, USA. As such, the results are likely to vary for a different year and at different geographic locations. We thus recommend that future studies (i) expand the number of field samples, (ii) assess the efficacy of modeling for different regions and even different beet cultivars, and (iii) consider the inclusion of structure-related features, such as those derived from light detection and ranging (LiDAR) and structure-from-motion techniques. Nevertheless, we were encouraged by our models'performance in the variability in yield explained (R 2 ) and the relatively low RMSE values. Acknowledgments: We credit Nina Raqueno and Tim Bauch of the RIT drone team, who were responsible for capturing the UAS imagery. We also thank Imergen Rosario and Kedar Patki for assistance with data collection.

Conflicts of Interest:
The authors declare no conflict of interest.