Predicting Table Beet Root Yield with Multispectral UAS Imagery

: Timely and accurate monitoring has the potential to streamline crop management, harvest planning, and processing in the growing table beet industry of New York state. We used unmanned aerial system (UAS) combined with a multispectral imager to monitor table beet ( Beta vulgaris ssp. vulgaris ) canopies in New York during the 2018 and 2019 growing seasons. We assessed the optimal pairing of a reﬂectance band or vegetation index with canopy area to predict table beet yield components of small sample plots using leave-one-out cross-validation. The most promising models were for table beet root count and mass using imagery taken during emergence and canopy closure, respectively. We created augmented plots, composed of random combinations of the study plots, to further exploit the importance of early canopy growth area. We achieved a R 2 = 0.70 and root mean squared error (RMSE) of 84 roots (~24%) for root count, using 2018 emergence imagery. The same model resulted in a RMSE of 127 roots (~35%) when tested on the unseen 2019 data. Harvested root mass was best modeled with canopy closing imagery, with a R 2 = 0.89 and RMSE = 6700 kg/ha using 2018 data. We applied the model to the 2019 full-ﬁeld imagery and found an average yield of 41,000 kg/ha (~40,000 kg/ha average for upstate New York). This study demonstrates the potential for table beet yield models using a combination of radiometric and canopy structure data obtained at early growth stages. Additional imagery of these early growth stages is vital to develop a robust and generalized model of table beet root yield that can handle imagery captured at slightly different growth stages between seasons.


Introduction
New York is a sentinel center of production for table beet (Beta vulgaris spp. vulgaris: Family Chenopodiaceae) in the USA, ranking second behind Wisconsin [1], and is undergoing exponential industry growth. The table beet industry in New York is diverse in enterprise size and scale, ranging from small, diversified growers that supply farm markets and roadside stands to broad-acre fields up to approximately 45 ha in size. Broad-acre production is for processing into cans and jars, direct fresh market sales, and feedstock as value-added beet products, such as snack packs and juices. The growth of the industry is fueled by an enhanced awareness of the health benefits of consuming table beets and beet-based products. These well-documented effects range from nutritional intervention to disease mitigation [2]. The most commonly reported physiological changes associated with table beet consumption include improvements in cardiovascular health and sugar metabolism attributed to the rich source of dietary nitrate [3,4], the antioxidant and antiinflammatory effects of betalains [5], and the potential to improve stamina in exercise [6]. 2017, with an R 2 = 0.816, during the sugarbeet V10 growth stage (approximately 28 June in western Minnesota) [19].
Canopy area is especially important for imagery obtained at early growth stages before the neighboring canopies have merged. Canopy height is also important but can be difficult to assess before the plants are several centimeters in height, depending on the ground sample distance (GSD) of the imagery. Overall, the early stages are well suited to assess the vigor of early crop growth via both the average radiometric properties of the canopy and the extent of the canopy growth in terms of canopy area and/or height [22][23][24][25].
In this work, we exclusively investigate the predictive power of UAS multispectral imagery for table beet yield components (root mass, count, and shoulder diameter size distributions) at harvest. This is a challenging research task, since, for many crops with aboveground fruits, the health or yield can be directly estimated with imagery, but for root vegetables, we can only infer their yield with imagery exclusively via their canopy structure and reflectance. We build on previous work by creating area-augmented data from our experimental plots in order to exploit the importance of structural metrics derived from imagery and to evaluate model performance on independent data.

Study Area
The Rochester Institute of Technology UAS and Cornell University agronomy teams collected multispectral imagery and ground truth agronomic data from table beet fields located in Batavia, New York, USA during July and August of 2018 and 2019 (see Figure 1). The commercial fields contained table beet cv. Merlin, planted in accordance with typical industry management practices, i.e., twin rows with approximately 10 cm between-plant spacing and~60 cm between-row spacing.
Remote Sens. 2021, 13, x FOR PEER REVIEW 3 of 20 and sites [20]. The best sugar beet root yield model was for combined NDVI and crop height at one site in 2017, with an = 0.816, during the sugarbeet V10 growth stage (approximately 28 June in western Minnesota) [19].
Canopy area is especially important for imagery obtained at early growth stages before the neighboring canopies have merged. Canopy height is also important but can be difficult to assess before the plants are several centimeters in height, depending on the ground sample distance (GSD) of the imagery. Overall, the early stages are well suited to assess the vigor of early crop growth via both the average radiometric properties of the canopy and the extent of the canopy growth in terms of canopy area and/or height [22][23][24][25].
In this work, we exclusively investigate the predictive power of UAS multispectral imagery for table beet yield components (root mass, count, and shoulder diameter size distributions) at harvest. This is a challenging research task, since, for many crops with aboveground fruits, the health or yield can be directly estimated with imagery, but for root vegetables, we can only infer their yield with imagery exclusively via their canopy structure and reflectance. We build on previous work by creating area-augmented data from our experimental plots in order to exploit the importance of structural metrics derived from imagery and to evaluate model performance on independent data.

Study Area
The Rochester Institute of Technology UAS and Cornell University agronomy teams collected multispectral imagery and ground truth agronomic data from table beet fields located in Batavia, New York, USA during July and August of 2018 and 2019 (see Figure  1). The commercial fields contained table beet cv. Merlin, planted in accordance with typical industry management practices, i.e., twin rows with approximately 10 cm betweenplant spacing and ~60 cm between-row spacing.

Data Collection
Fifty 3 m long plots were marked soon after crop emergence (5 July 2018), as well as 20, 1.5 m long plots during the following season (15 July 2019). Four and two data collections of agronomic ground truth and UAS canopy reflectance data were collected in 2018 and 2019, respectively. An additional ground truth harvest collection, without imagery,

Data Collection
Fifty 3 m long plots were marked soon after crop emergence (5 July 2018), as well as 20, 1.5 m long plots during the following season (15 July 2019). Four and two data collections of agronomic ground truth and UAS canopy reflectance data were collected in 2018 and 2019, respectively. An additional ground truth harvest collection, without imagery, was also conducted in 2019. Crop stands were assessed for both years by counting the plants in the entire length of each of the double-row plots on each occasion. Canopy reflectance imagery were obtained as close as possible to the ground truth collection dates via a DJI Matrice 100 UAS vehicle, mounted with a MicaSense RedEdge-M five-band multispectral visible-near-infrared (VNIR) camera system. Due to instrument availability, canopy reflectance imagery was obtained using a DJI Matrice 600 UAS equipped with a Headwall Photonics Nano-Hyperspec imaging spectrometer (272 spectral bands; 400-1000 nm) for the second observation in 2018, in lieu of the multispectral imagery. At harvest, in addition to plant/root counts, a number of yield parameters were collected in the 2018 season by hand removal of all plants within the plots (Table 1). These included root number, root mass (kg), root shoulder diameters (mm), and the dry weight of the foliage separated from the roots. Foliage was weighed in the field and a subsample (approximately 10% of the total weight) was dried at 60 • C for 48 h to calculate dry weight. Roots over 20 mm in diameter were counted and weighed, and returned to Geneva, NY for storage at 10 • C for up to 10 days, until shoulder diameters were measured using digital calipers (Mitutoyo, USA). Only root count was conducted for the 2019 harvest. The MicaSense RedEdge-M camera captured images in five narrowband spectral channels (blue, green, red, red edge, and near-infrared [26]; Table 2). The flight ground sample distances ranged from 0.5 to 10 cm/pixel. We down-sampled the hyperspectral imaging canopy reflectance spectra for each pixel to simulate the standard multispectral imagery for the second observation in 2018. We integrated the pixel spectra, multiplied with a Gaussian of peak central wavelength and full width at half maximum (FWHM) to down-sample the spectra (Table 2).

Data Preprocessing
We registered, orthorectified, and calibrated all UAS multispectral imagery obtained with the MicaSense camera using the Pix4D software (V.4.4.12) [27]. The software processing ingests the raw frames from each spectral band and outputs co-registered and radiometrically-corrected reflectance orthomosaics. GPS ground control points were deployed in the fields for manual identification of tie-points to verify and improve band and image registration.
The calibration procedure for the MicaSense camera uses sensor and image specific metadata parameters for conversion from raw digital counts to radiance (W/m 2 /sr/nm). These include parameters for a vignette correction function, radiometric calibration coefficients, sensor gain, exposure time, and black level values [28,29]. Additionally, the MicaSense camera was paired with a 5-band downwelling light sensor (DLS 2) to capture the downwelling radiance in sync with the radiance imagery of the crop canopy [30]. The Pix4D software uses the primary camera parameters, sun irradiance measured with DLS, and sun angle from the DLS inertial measurement unit for radiometric correction [29]. We then stacked the resulting reflectance mosaics for each band into one multispectral image, using ENVI software (V.5.5) [31].
The 2018 observations with hyperspectral imagery were pre-processed for orthorectification and radiance calibration/conversion using Headwall's Hyperspec III SpectralView software [32]. The spectral imagery captured black, dark gray, light gray, and white calibration panels deployed on the north side of the 2018 site. We extracted the calibration panel spectra using ENVI (V.5.5 [31]). The panels have well-characterized ground truth reflectance spectra that were used to convert the pre-processed radiance flight imagery into reflectance. We used the Empirical Line Method (ELM) tool in ENVI [33] to determine the relationship between the calibration panel UAS radiance spectra and their true reflectance spectra and determined linear coefficients to convert between radiance and reflectance at each wavelength. This process removes both solar irradiance and atmospheric path radiance.
Plots in the field were marked with red plates to rapidly locate and define regions of interest (ROIs) for each plot using ENVI. The plots differ in pixel dimensions for each set of imagery, depending on the ground sample distance of the UAS-based MicaSense imagery. We extracted 5-band plot images from each of the ROIs and used the Spectral Python library SPy (V.0.20) [34] to import the ENVI-formatted plot images to NumPy Arrays in Python for further analysis. Figure 2 displays a sample of the extracted plot imagery from the 9 July 2018 flight (shortly following crop emergence). We noted that plots 41, 42, 43, and 46 exhibited substantial weed growth at this stage and were removed from the analysis. Additionally, plot 27 had an anomalously high recorded foliage mass (>2× of other plots) and was removed as an outlier/erroneous recording. The next step required accurate identification of vegetation pixels for further analysis.

Canopy Pixel Segmentation
Several methods may be used to identify the pixels containing vegetation in plot imagery. We clearly distinguished the crop canopy from background soil, prior to the growth reaching a state of overlapped row canopies, due to the high spatial resolution of the UAS imagery. Common classification practices include applying an appropriate threshold to the spectral angle map (SAM) of an endmember spectrum [35] or using k-means clustering to identify and classify spectral endmembers [36]. However, band reflectance can change substantially with canopy growth, causing these methods to provide inconsistent results amongst the various datasets. We thus opted for a relatively simple VI threshold to identify vegetation pixels in each plot. Similar to Barzin et al. [16], we determined that while a Normalized Difference Vegetation Index (NDVI) threshold is commonly used for soil masking [37,38], it is not capable of producing a consistent result amongst the imagery from varying growth stages using a single threshold value. We used

Canopy Pixel Segmentation
Several methods may be used to identify the pixels containing vegetation in plot imagery. We clearly distinguished the crop canopy from background soil, prior to the growth reaching a state of overlapped row canopies, due to the high spatial resolution of the UAS imagery. Common classification practices include applying an appropriate threshold to the spectral angle map (SAM) of an endmember spectrum [35] or using kmeans clustering to identify and classify spectral endmembers [36]. However, band reflectance can change substantially with canopy growth, causing these methods to provide inconsistent results amongst the various datasets. We thus opted for a relatively simple VI threshold to identify vegetation pixels in each plot. Similar to Barzin et al. [16], we determined that while a Normalized Difference Vegetation Index (NDVI) threshold is commonly used for soil masking [37,38], it is not capable of producing a consistent result amongst the imagery from varying growth stages using a single threshold value. We used the triangular vegetation index (TVI), because the TVI of the canopy exhibited a limited variation over the cropping season compared to NDVI.
The traditional formula for TVI, introduced by Broge and Leblanc [39], incorporates the green, red, and NIR bands as: with reflectances = 550 nm, = 670 nm, and = 750 nm. This index essentially measures the area of the triangle, as defined by the peak green and red absorption trough, due to chlorophyll absorption, and the NIR plateau (inter-cellular leaf structure) in vegetation spectra [39]. We used the green, red, and red edge bands for vegetation classification ( Table 2). The traditional TVI values are, therefore, expected to represent an increase in the triangular area as the canopy grows and NIR reflectance increases with more layers of vegetation. We chose to use the red edge band, instead of NIR, in order to limit the difference in the index between growth stages. The median reflectance for each band The traditional formula for TVI, introduced by Broge and Leblanc [39], incorporates the green, red, and NIR bands as: with reflectances ρ green = 550 nm, ρ red = 670 nm, and ρ N IR = 750 nm. This index essentially measures the area of the triangle, as defined by the peak green and red absorption trough, due to chlorophyll absorption, and the NIR plateau (inter-cellular leaf structure) in vegetation spectra [39]. We used the green, red, and red edge bands for vegetation classification ( Table 2). The traditional TVI values are, therefore, expected to represent an increase in the triangular area as the canopy grows and NIR reflectance increases with more layers of vegetation. We chose to use the red edge band, instead of NIR, in order to limit the difference in the index between growth stages. The median reflectance for each band in all plot canopy imagery from 2018 are shown in Figure 3 (left) along with the complete average spectra (right). The canopy near-infrared reflectance is clearly the most variable throughout the growing season. The red edge band reflectance observed here only decreases slightly over the cropping season, as the red edge inflection point shifts to longer wavelengths with an increase in vegetation density [40]. This choice is convenient for the purpose of classification, since it allows use of a consistent threshold for every dataset. Other common indices, like NDVI, increase for canopy-level sensing throughout the growing season and require an adjustment to the threshold for each dataset to consistently classify all vegetation appropriately.
throughout the growing season. The red edge band reflectance observed here only decreases slightly over the cropping season, as the red edge inflection point shifts to longer wavelengths with an increase in vegetation density [40]. This choice is convenient for the purpose of classification, since it allows use of a consistent threshold for every dataset. Other common indices, like NDVI, increase for canopy-level sensing throughout the growing season and require an adjustment to the threshold for each dataset to consistently classify all vegetation appropriately. We included a comparison of the TVI threshold, SAM, k-means, and NDVI threshold methods on a sample plot over the course of the 2018 season to further reinforce our choice of the modified TVI for classification ( Figure 4). The plot's RGB composite and vegetation classification masks of each method are shown for each growth stage. The k-means and SAM masks were determined using the kmeans and msam [35] functions in the SPy library. The NDVI and SAM methods both classified the entire plot image as vegetation in the two final observations, but a stricter (higher) threshold for either method failed to classify the first observation's vegetation pixels, due to the less-developed canopy. The kmeans method is strongly dependent on the number of clusters chosen. In the later observations, it was optimal to use fewer clusters, because vegetation covered nearly the entire plot area, resulting in less spectral diversity. However, in early observations there was more variation in the plot image due to increased soil visibility. The image from 9 July 2018 required at least three clusters to properly classify only vegetation cover as actual vegetation. Thus, these three methods failed to classify the plot imagery of all four growth stages using a standard set of method parameters. Only the TVI threshold effectively classified all four growth stage images with a single threshold value for all four sets of imagery. We included a comparison of the TVI threshold, SAM, k-means, and NDVI threshold methods on a sample plot over the course of the 2018 season to further reinforce our choice of the modified TVI for classification ( Figure 4). The plot's RGB composite and vegetation classification masks of each method are shown for each growth stage. The k-means and SAM masks were determined using the kmeans and msam [35] functions in the SPy library. The NDVI and SAM methods both classified the entire plot image as vegetation in the two final observations, but a stricter (higher) threshold for either method failed to classify the first observation's vegetation pixels, due to the less-developed canopy. The k-means method is strongly dependent on the number of clusters chosen. In the later observations, it was optimal to use fewer clusters, because vegetation covered nearly the entire plot area, resulting in less spectral diversity. However, in early observations there was more variation in the plot image due to increased soil visibility. The image from 9 July 2018 required at least three clusters to properly classify only vegetation cover as actual vegetation. Thus, these three methods failed to classify the plot imagery of all four growth stages using a standard set of method parameters. Only the TVI threshold effectively classified all four growth stage images with a single threshold value for all four sets of imagery.

Feature Choice
Using the vegetation classification masks for each plot, we calculated a range of statistical and morphological features. To avoid overfitting issues caused by small sample

Feature Choice
Using the vegetation classification masks for each plot, we calculated a range of statistical and morphological features. To avoid overfitting issues caused by small sample sizes with high dimensionality [41], we limited our study to individual reflectance bands or biophysically-sensible VI features and single radiometric features paired with canopy area. We tested only VIs known to exhibit a priori phenological correlations with crop vigor, and from these constructed our final models only with those features exhibiting optimal posteriori correlations.
In addition to the five-band multispectral imagery, we calculated several VI images for each plot (Table 3; Figure 5). The indices used here were specifically chosen to capture different aspects of the plots' spectra via incorporation of different combinations of the five camera bands. The basic Difference Vegetation Index (DVI) effectively separates the vegetation from the soil, but is sensitive to background soil reflectance [42]. We also included the similar Green Difference Vegetation Index (GDVI), designed for predicting nitrogen requirements in corn [43]. The most commonly-used vegetation index, NDVI, is related to canopy structure, Leaf Area Index (LAI), and photosynthesis [44], but is adversely sensitive to the background soil, atmospheric conditions, and shadowing. The Enhanced Vegetation Index (EVI) [45] attempts to correct for the soil and atmospheric effects in NDVI and offers more contrast in developed canopies, where NDVI may saturate due to a higher Leaf Area Index (LAI). The Modified Soil Adjusted Vegetation Index 2 (MSAVI2) [46] reduces noise due to soil and increases the dynamic range of the canopy signal, while the Green Chlorophyll Index (GCI) provides a better relation to chlorophyll content than NDVI [47]. The Modified Triangular Vegetation Index (MTVI), used for LAI estimation [48], is slightly different from the TVI we used earlier for vegetation classification, effectively swapping the NIR band for the red edge. We included the Visible-Band Difference Vegetation Index (VDVI) [49], which incorporates the blue, green, and red channels, to assess the effectiveness of using only visible channels. Finally, we included the Red Edge Normalized Difference Vegetation Index (RENDVI) [50] to capture variation in the vegetation spectra along the red edge.
Green Chlorophyll Index [47] N IR We next determined the median reflectance of each band and median value of each VI for the vegetation-classified pixels in each plot. Although the pixel reflectance of each band did not deviate significantly from a Gaussian distribution, we used the median in lieu of the mean to eliminate the potential of an anomalous pixel skewing a plot's reflectance or index values. We calculated the physical area of vegetation pixels for each plot, by multiplying the number of TVI > 7.5 threshold vegetation-classified pixels by the squared ground sampling distance. This yielded 15 total features (five bands, nine indices, and one area) for each plot to be used in data analysis for each growth stage imagery set.  We next determined the median reflectance of each band and median value of each VI for the vegetation-classified pixels in each plot. Although the pixel reflectance of each band did not deviate significantly from a Gaussian distribution, we used the median in lieu of the mean to eliminate the potential of an anomalous pixel skewing a plot's reflectance or index values. We calculated the physical area of vegetation pixels for each plot, by multiplying the number of TVI > 7.5 threshold vegetation-classified pixels by the squared ground sampling distance. This yielded 15 total features (five bands, nine indices, and one area) for each plot to be used in data analysis for each growth stage imagery set.

Data Analysis
We investigated multiple linear regression (MLR) models using individual radiometric features combined with canopy area to model the various yield parameters with datasets from each growth stage. We used leave-one-out cross-validation (LOOCV) and for each MLR model (scikit-learn version 0.21.3) to determine the best-per-

Data Analysis
We investigated multiple linear regression (MLR) models using individual radiometric features combined with canopy area to model the various yield parameters with datasets from each growth stage. We used leave-one-out cross-validation (LOOCV) R 2 and RMSE for each MLR model (scikit-learn version 0.21.3) to determine the best-performing growth stage and radiometric feature for each yield component. We also summarize LOOCV results for individual feature linear models of canopy area and the best-performing radiometric feature for comparison with the optimal MLR model. Finally, we built upon the most promising growth stage models by including additional augmented plots with larger areas.
We created augmented area plots for the emergence and canopy closing growth stages using the original N plots. An additional N plots (the double set) consisted of two randomly-chosen and combined plots from the original set, and a third set of N plots (the triple set) consisted of a randomly-chosen combination of three original plots. This yielded three times as many plots at each growth stage, with enhanced canopy area variability. The radiometric data were extracted from these double and triple plots in the same manner as the original plot data. While the original plot models naturally favored canopy area in some cases, the use of augmented area plots forced the area feature to be most important. We repeated the LOOCV assessment of individual radiometric features paired with canopy area using the area-augmented models for the previously identified optimal growth stage. We applied the final root count model to the 2019 plots nearest in growth stage for direct testing and the root mass model to the whole field for a final coarse assessment of the model's performance on unseen data, comparing the outcome to a standard yield estimate in terms of kg/ha for this region.
Prior to using MLR, we centered and scaled the dataset by removing the mean and scaling to unit variance, so that the variety of scalar value ranges of the individual input features did not skew the model coefficients higher for certain features [52]. In addition to reporting LOOCV we provided visual assessments of model performance (measurements vs. predictions) by training models with the optimal feature combination using 90% of the plot samples in the selected dataset. We also report RMSE for the 10% testing data prediction of these models.

Table Beet Root Count
The LOOCV scores for root count modeled with the standard plot imagery from each growth stage are compiled in Table 4. The highest R 2 and lowest RMSE were found for the emergence imagery. The canopy area (R 2 = 0.23, RMSE = 25 roots) alone was only slightly informative at this growth stage and did not improve the predictive ability of the NIR band (R 2 = 0.40, RMSE = 22 roots) when combined in the MLR model (R 2 = 0.38, RMSE = 22 roots). The optimal RMSE of 22 roots is 13% of the mean root count across all plots. The LOOCV scores for each band and VI paired with canopy area are shown in Figure 6 for the area-augmented emergence imagery. The VIs generally outperformed the individual reflectance bands, and nearly all combinations provided an improvement to the canopy area alone (R 2 = 0.65, RMSE = 91 roots). The best pair was DVI and canopy area with R 2 = 0.70 and RMSE = 84 roots (24% of the mean root count for all area-augmented plots). A sample of DVI + Area prediction vs. measured root counts is shown in the left panel of Figure 7 using 90% training and 10% testing data. While the radiometric features individually provided little predictive ability for the area-augmented dataset, their combination with canopy area clearly provides an improvement on area alone. The predictions for larger augmented plots are promising and distributed normally, but the predictions of the original plots are skewed to overpredicting in the majority of cases.     The model is trained with 90% (red) and tested with the remaining 10% (green) area augmented study plot imagery. RMSE are displayed for the training data using leave-one-out cross-validation and for the unseen test data. The DVI + Area model is trained with 100% of the 2018 emergence imagery and tested to predict the root count (blue) for 100% of the 2019 area-augmented study plot imagery occurring closest to emergence (right).
We further trained the DVI + Area model using the 2018 area-augmented emergence imagery and tested on the area-augmented imagery from 2019 obtained on 16 July, closest to emergence for that season. The predictions for this model are shown in the right panel of Figure 7. The model underpredicts several of the larger augmented plots with an = 127 roots or 35% of the mean root count for the 2019 area-augmented plots. The model is trained with 90% (red) and tested with the remaining 10% (green) area augmented study plot imagery. RMSE are displayed for the training data using leave-one-out cross-validation and for the unseen test data. The DVI + Area model is trained with 100% of the 2018 emergence imagery and tested to predict the root count (blue) for 100% of the 2019 area-augmented study plot imagery occurring closest to emergence (right).
We further trained the DVI + Area model using the 2018 area-augmented emergence imagery and tested on the area-augmented imagery from 2019 obtained on 16 July, closest to emergence for that season. The predictions for this model are shown in the right panel of Figure 7. The model underpredicts several of the larger augmented plots with an RMSE = 127 roots or 35% of the mean root count for the 2019 area-augmented plots.

Beet Root Mass 2018
The LOOCV scores for root mass modeled with the standard plot imagery from each growth stage are compiled in Table 5. The highest R 2 and lowest RMSE were found for the canopy closing imagery. However, the improvements to RMSE between growth stages were minimal. The canopy area (R 2 = 0.36, RMSE = 0.9 kg) alone was most informative at the canopy closing growth stage while VDVI (R 2 = −0.03, RMSE = 1.2 kg) when combined in the MLR model (R 2 = 0.37, RMSE = 0.9 kg) provided meager improvement. The optimal RMSE of 0.9 kg is 10% of the mean root mass for all plots (or approximately 4800 kg/hectare). The LOOCV scores for each band and VI paired with canopy area are shown in Figure 8 for the area-augmented canopy closing imagery. The best pairing was VDVI and canopy area, but all reflectance bands and VIs offered no improvement to canopy area alone (R 2 = 0.89, RMSE = 2.5 kg per area-augmented plot or approximately 6700 kg/hectare). The total root mass amongst the plots exhibited less variation than the root count and appear as three distinct clusters after augmentation-original plots, double combo-plots, and triple-combo plots (Figure 9). Capturing the relationship between radiometric data and root mass may require plots with a broader range of root mass. The LOOCV scores for each band and VI paired with canopy area are shown in Figure 8 for the area-augmented canopy closing imagery. The best pairing was VDVI and canopy area, but all reflectance bands and VIs offered no improvement to canopy area alone ( = 0.89, = 2.5 kg per area-augmented plot or approximately 6700 kg/hectare). The total root mass amongst the plots exhibited less variation than the root count and appear as three distinct clusters after augmentation-original plots, double comboplots, and triple-combo plots (Figure 9). Capturing the relationship between radiometric data and root mass may require plots with a broader range of root mass.

Figure 8. The area-augmented table beet root mass 2 (left) and
(right) leave-one-out cross-validation performance (y-axis) for individual multiple linear regression models. Each model consists of one reflectance band or vegetation index shown on the x-axis combined with canopy area. The VDVI paired with canopy area was the best-performing model with 2 = 0.89 and = 2.5 kg. However, the vertical axis scaling shown for both plots indicates that the radiometric data provide no significant improvement to the table beet root mass model derived from the area-augmented imagery. Figure 8. The area-augmented table beet root mass R 2 (left) and RMSE (right) leave-one-out cross-validation performance (y-axis) for individual multiple linear regression models. Each model consists of one reflectance band or vegetation index shown on the x-axis combined with canopy area. The VDVI paired with canopy area was the best-performing model with R 2 = 0.89 and RMSE = 2.5 kg. However, the vertical axis scaling shown for both plots indicates that the radiometric data provide no significant improvement to the table beet root mass model derived from the area-augmented imagery.
Remote Sens. 2021, 13, x FOR PEER REVIEW 13 of 20 Figure 9. Measured vs. predicted table beet root mass for the 2018 canopy closing VDVI + Area multiple linear regression model based on area-augmented plot imagery. The model is trained with 90% (red) and tested with the remaining 10% (green) area-augmented study plot imagery. RMSE are displayed for the training data using leave-one-out cross-validation and for the unseen test data.
The average attainable table beet root yield in New York is estimated at 40,000 kg/hectare (J . Henderson, personal communication). We, therefore, applied our area-augmented table beet root mass model to new data for a large section (~0.3 ha) of the 2019 field imagery to assess the yield for the entire image coverage. Again, the 2019 data used here did not coincide with the exact growth stage of the 2018 canopy closing dataset used to build the area-augmented beet root mass model. Regardless, the results of our model applied to the 2019 field are shown in Figure 10. We predicted a total yield of 41,000 kg/ha, where the average grid cell (~4 m 2 ) had table beet roots totaling 17 ± 2.5 kg, by tallying the root mass from each grid cell and dividing by the total area. Considering the differ- Figure 9. Measured vs. predicted table beet root mass for the 2018 canopy closing VDVI + Area multiple linear regression model based on area-augmented plot imagery. The model is trained with 90% (red) and tested with the remaining 10% (green) area-augmented study plot imagery. RMSE are displayed for the training data using leave-one-out cross-validation and for the unseen test data.
The average attainable table beet root yield in New York is estimated at 40,000 kg/hectare (J . Henderson, personal communication). We, therefore, applied our area-augmented table beet root mass model to new data for a large section (~0.3 ha) of the 2019 field imagery to assess the yield for the entire image coverage. Again, the 2019 data used here did not coincide with the exact growth stage of the 2018 canopy closing dataset used to build the area-augmented beet root mass model. Regardless, the results of our model applied to the 2019 field are shown in Figure 10. We predicted a total yield of 41,000 kg/ha, where the average grid cell (~4 m 2 ) had table beet roots totaling 17 ± 2.5 kg, by tallying the root mass from each grid cell and dividing by the total area. Considering the differences in the canopy growth used to construct and test the model, this performance should merely be taken as a moderate prediction, rather than proof of concept.
90% (red) and tested with the remaining 10% (green) area-augmented study plot imagery. RMSE are displayed for the training data using leave-one-out cross-validation and for the unseen test data.
The average attainable table beet root yield in New York is estimated at 40,000 kg/hectare (J . Henderson, personal communication). We, therefore, applied our area-augmented table beet root mass model to new data for a large section (~0.3 ha) of the 2019 field imagery to assess the yield for the entire image coverage. Again, the 2019 data used here did not coincide with the exact growth stage of the 2018 canopy closing dataset used to build the area-augmented beet root mass model. Regardless, the results of our model applied to the 2019 field are shown in Figure 10. We predicted a total yield of 41,000 kg/ha, where the average grid cell (~4 m 2 ) had table beet roots totaling 17 ± 2.5 kg, by tallying the root mass from each grid cell and dividing by the total area. Considering the differences in the canopy growth used to construct and test the model, this performance should merely be taken as a moderate prediction, rather than proof of concept.

Beet Root Diameter 2018
The LOOCV scores for root mass modeled with the standard plot imagery from each growth stage are compiled in Table 6. While the emergence data had the highest correlations to root diameter the correlations were lower than for root count and mass. The NIR band was most correlated with root diameter (R 2 = 0.22, RMSE = 3.2 mm), but the canopy area (R 2 = 0.02, RMSE = 3.5 mm) offered no improvement in the combined MLR model (R 2 = 0.22, RMSE = 3.2 mm). The average root diameter was not explored further in this study with area-augmented imagery because the average root diameter does not scale in the same manner as table beet root count and mass. Table 6. Table beet root diameter growth stage and band/VI assessment report with leave-one-out cross-validation R 2 and root mean squared error (RMSE) for the best pairing of a radiometric feature with canopy area MLR model and the linear regression scores of the individual features in New York in 2018. The average root diameter exhibited poor correlation with radiometric features and the canopy area. However, when we separate the table root counts into standard diameter ranges (John Henderson, private communication), the three diameter ranges nearest the center of the diameter distribution (10-25 mm, 25-44 mm, and 44-63 mm), had the highest NIR + Area MLR model correlations (LOOCV R 2 = 0.24, 0.37, and 0.26, respectively) with emergence plot imagery. This is expected, as they make up the largest percentages of the total root count. We show the percentage of roots falling in each of the processing diameter ranges in Figure 11, along with target percentages. The actual numbers from the test plots indicated over-production of smaller diameter and less production of larger diameter beets. Additional emergence data may be used to develop a root size distribution model in the future, but currently the low numbers for the 63-75 mm range prevent meaningful size distribution modeling in the three ranges that are key to a typical grower target. with black open circles and fall outside the error bars located at Quantile 3 + 1.5 IQR (above) and Quantile 1-1.5 IQR (below). The optimal percentages of table beet roots falling in each range, as provided by Love Beets USA (a commercial beet grower), is plotted with red filled circles.

Foliage Mass 2018
The LOOCV scores for root mass modeled with the standard plot imagery from each growth stage are compiled in Table 7. The harvest growth stage provided the best correlation with the dry weight of foliage. The combination of VIs and canopy area produced the best MLR model ( = 0.32, = 63 g), but neither the individual RENDVI ( = 0.06, = 74 g) or canopy area ( = 0.13, = 71 g) model had comparable performance.
The canopy area for all plots varied by the smallest margin at harvest (0.35 ). Ultimately, all plots from 2018 reached maximum canopy area by the canopy-closed growth stage. Therefore, the canopy area loses value for predicting yield components. Moreover, any variation explained by canopy area at later growth stages is likely due to differences in pixel shading of the plots, rather than actual area. Later growth stage models should rely most heavily on the radiometric features, but these also reach a point of saturation or lack of variation between plots in most cases, as a function of closed-canopy table beets. It has been shown that typical red-near-infrared vegetation indices, such as NDVI, saturate at high Leaf Area Index (LAI) values [53]. Although we did not measure LAI explicitly, it is likely that closed canopies result in a reduction in predictive power for such red-nearinfrared vegetation indices. Table 7. Table beet growth stage and band/VI assessment report with leave-one-out cross-validation R 2 and root mean squared error (RMSE) for the best pairing of a radiometric feature with canopy area MLR model and the linear regression scores of the individual features in New York in 2018. Figure 11. The median percentage of all 2018 plots' table beet roots, falling in each diameter range, are plotted in green inside of a box spanning the inter quantile range (IQR). Outliers are marked with black open circles and fall outside the error bars located at Quantile 3 + 1.5 IQR (above) and Quantile 1-1.5 IQR (below). The optimal percentages of table beet roots falling in each range, as provided by Love Beets USA (a commercial beet grower), is plotted with red filled circles.

Foliage Mass 2018
The LOOCV scores for root mass modeled with the standard plot imagery from each growth stage are compiled in Table 7. The harvest growth stage provided the best correlation with the dry weight of foliage. The combination of VIs and canopy area produced the best MLR model (R 2 = 0.32, RMSE = 63 g), but neither the individual RENDVI (R 2 = 0.06, RMSE = 74 g) or canopy area (R 2 = 0.13, RMSE = 71 g) model had comparable performance. The canopy area for all plots varied by the smallest margin at harvest (0.35 m 2 ). Ultimately, all plots from 2018 reached maximum canopy area by the canopy-closed growth stage. Therefore, the canopy area loses value for predicting yield components. Moreover, any variation explained by canopy area at later growth stages is likely due to differences in pixel shading of the plots, rather than actual area. Later growth stage models should rely most heavily on the radiometric features, but these also reach a point of saturation or lack of variation between plots in most cases, as a function of closed-canopy table beets. It has been shown that typical red-near-infrared vegetation indices, such as NDVI, saturate at high Leaf Area Index (LAI) values [53]. Although we did not measure LAI explicitly, it is likely that closed canopies result in a reduction in predictive power for such red-near-infrared vegetation indices.
We did not test foliage models on the 2019 data, because we did not obtain imagery of the canopy-closed and harvest growth stages or gather the foliage mass ground truth for that year. Additionally, measured foliage mass was only modestly correlated with root mass and root count (R 2 = 0.38 and 0.23, respectively). We see limited value in the development of a precise foliar biomass model without first expanding our dataset to more clearly exploit any potential relationship between foliar biomass and root yield.

Discussion
The performance of the 2018 area-augmented root count and mass models assessed on the 2019 data is promising, especially considering the slight differences in growth stages, in both cases. In fact, the median DVI values for the 2019 plots ranged from 0.28 to 0.49 compared to 0.20 to 0.32 in 2018, indicating that the earliest 2019 imagery captured a more advanced stage of canopy growth than the 2018 emergence data used to develop the root count model. Additionally, the smallest plots used to construct the model were twice as long as the 2019 plots. We thus recommend that additional imagery, collected on several days following emergence, could help pinpoint an optimal time to capture imagery for table beet root count predictions. Alternatively, more data may allow the use of a temporal feature representing the growth stage in our prediction model.
Similarly, there is significant potential to develop improved root mass models with additional data collected at different days during the emergence and canopy closing growth stages. Even though the Band/VI + Area root mass models did not improve upon the Area-only model, using the radiometric properties of the canopy to quantify area before closure is an effective method to predict table beet root yield.
The physical size of the canopy growth arguably should have the strongest correlation with yield, when measured close to harvest. However, by the canopy-closed growth stage observations in 2018, every plot canopy had merged with neighboring rows, thereby resulting in virtually identical canopy area. The radiometric (image) data also lacked the significant variability needed to predict a variety of root yield components for the plots by the canopy-closed growth stage. In other words, it appears that eventually the plots reach a state of maturity where they become less distinguishable in terms of vegetated area and their radiometric properties. This finding also has been demonstrated in previous studies, where especially red-NIR vegetation indices saturate and thus become less sensitive under high Leaf Area Index (LAI) or canopy-closed scenarios [53,54].
The need for early growth stage imagery is highlighted by Figure 12, where we show the radiometric variance of each input feature for the 2018 study, divided by the mean of the feature for each of the growth stages. The emergence and canopy closing growth stages exhibited the highest ratio of feature variance, normalized by the mean. Many of the radiometric features showed minimal variance (<2%) across plots. All of the canopy-closed and harvest stage radiometric features exhibited minimal variation amongst the plots. The canopy area showed the clearest trend of decreasing variation as the canopy grew into the between-row space. We, therefore, contend that monitoring the early growth stages more closely and determining when each plot reaches maturity may be more relevant to root development and ultimately, yield component prediction. In addition to acquiring imagery more frequently, ground truth data should be obtained from a larger sample size. A larger number of smaller scale samples covering an identical area could be accomplished with equivalent sampling labor to reach a higher base sample size. The area-augmented procedure used here could then be applied to model scalable yield components. A future study may also use hyperspectral data derived from an UAS or handheld spectrometer, with the goal of determining a set of spectral bands that exhibit optimal association with yield components [55], thereby improving upon the multispectral bands used in this study.
Plots with obvious visual evidence of weeds were removed as outliers before regression modeling. We, therefore, would need to incorporate a method for identifying and mapping weeds between rows in the early growth stages, as implemented by Pérez-Ortiz et al. [56] for sunflower crops, if we were to apply our models to a table beet field. An early, targeted elimination of weeds in problem areas may prevent the low root yields that we observed in four test plots that contained obvious examples of weed growth between rows in the emergence 2018 imagery (i.e., plots 41-43 and 46; Figure 2).

Conclusions
We evaluated the use of five-band visible-near-infrared UAS imagery to assess New York table beet root yield (root count and mass) across four airborne campaigns spanning the 2018 cropping season. Overall, our investigations of table beet root yield uncovered the need to obtain additional data, at different times throughout the early growth stages, and to prioritize flight plans for also obtaining canopy heights from our imagery using structure-from-motion (photogrammetry) or LiDAR data. We observed modest predictive ability for both the root count ( = 0.38, = 22 roots) and root mass models ( = 0.37, = 0.9 kg) using UAS multispectral imagery. However, extending 2018 singleyear models to independent data, e.g., full 2019 field or 2019 plot-level imagery, was successful after incorporating area-augmented plots (vegetation coverage/plot). These models took advantage of the obvious relationship between canopy area and root yield/count. The inclusion of 3D structural features, such as canopy heights from either light detection and ranging (LiDAR) or via image-based structure-from-motion (photogrammetry) may improve predictions for future studies. We might effectively exploit the imaging data to determine both canopy area and height, supplemented with the radiometric data when it is most variable early in the season to provide early predictions. Even if the canopy area and vegetation indices have reached saturation, canopy height may continue to exhibit variability between plots at the later growth stages and provide a more direct measure of the foliar biomass.
In addition to acquiring imagery more frequently, ground truth data should be obtained from a larger sample size. A larger number of smaller scale samples covering an identical area could be accomplished with equivalent sampling labor to reach a higher base sample size. The area-augmented procedure used here could then be applied to model scalable yield components. A future study may also use hyperspectral data derived from an UAS or handheld spectrometer, with the goal of determining a set of spectral bands that exhibit optimal association with yield components [55], thereby improving upon the multispectral bands used in this study.
Plots with obvious visual evidence of weeds were removed as outliers before regression modeling. We, therefore, would need to incorporate a method for identifying and mapping weeds between rows in the early growth stages, as implemented by Pérez-Ortiz et al. [56] for sunflower crops, if we were to apply our models to a table beet field. An early, targeted elimination of weeds in problem areas may prevent the low root yields that we observed in four test plots that contained obvious examples of weed growth between rows in the emergence 2018 imagery (i.e., plots 41-43 and 46; Figure 2).

Conclusions
We evaluated the use of five-band visible-near-infrared UAS imagery to assess New York table beet root yield (root count and mass) across four airborne campaigns spanning the 2018 cropping season. Overall, our investigations of table beet root yield uncovered the need to obtain additional data, at different times throughout the early growth stages, and to prioritize flight plans for also obtaining canopy heights from our imagery using structure-from-motion (photogrammetry) or LiDAR data. We observed modest predictive ability for both the root count (R 2 = 0.38, RMSE = 22 roots) and root mass models (R 2 = 0.37, RMSE = 0.9 kg) using UAS multispectral imagery. However, extending 2018 single-year models to independent data, e.g., full 2019 field or 2019 plot-level imagery, was successful after incorporating area-augmented plots (vegetation coverage/plot). These models took advantage of the obvious relationship between canopy area and root yield/count.
In reality, any application of a crop yield model based on radiometric information from satellites, UAS, or over-the-row tractor equipment needs to be robust to a multitude of scales and conditions. A model developed solely from radiometric data of a standard plot size likely will not be ideal if the goal is to apply these results in the field. Our area-augmented modeling demonstrates potential to create a model robust to this shortcoming by using a combination of radiometric and spatial (coverage or structural) data. We thus recommend that more data from an array of growth stages and seasons be included in future studies to develop more confidence in the operational UAS-based yield modeling applications.