Evaluation of Leaf Area Index (LAI) of Broadacre Crops Using UAS-Based LiDAR Point Clouds and Multispectral Imagery

—Leaf area index (LAI) is an established structural variable that reﬂects the three-dimensional (3-D) leaf layering of vegetation in response to environmental inputs. In this context, unmanned aerial system (UAS) based methods present a new approach to such plant-to ﬁeld-scale LAI assessment for precision agriculture applications. This article used UAS-based light detection and ranging (LiDAR) data and multispectral imagery (MSI) as twomodalitiestoevaluatetheLAIofasnapbeanﬁeld,towardeven-tualyieldmodelinganddiseaseriskassessment.LiDAR-derivedandMSI-derivedmetricswerefedtomultiplebiophysical-based andregressionmodels.ThecorrelationbetweenthederivedLAIandﬁeld-measuredLAIwassigniﬁcant.SixLiDAR-derivedmet-ricswereﬁtineightmodelstopredictLAI,amongwhichthesquarerootofthelaserpenetrationindexachievedthemostaccurate predictionresult( R 2 = 0.61, n RMSE = 19%). The MSI-derived models, which contained both structural features and spectral sig- natures, provided similar predicting effectiveness, with predicted R 2 ≈ 0.5 and n RMSE ≈ 22%. We furthermore observed variation in model effectiveness for different cultivars, different cultivar groups, and different UAS ﬂight altitudes, for both the LiDAR and MSI approaches. For data collected at a consistent ﬂight altitude, MSI-derived models could even exceed LiDAR-derived models, in terms of accuracy. This ﬁnding could support the possibility of replacing LiDAR with more cost-effective MSI-based approaches. However,LiDARremainsaviablemodality,sinceaLiDAR-derived3-Dmodelonlyrequiredasinglepredictorvariable,whilean MSI-derivedmodelreliedonmultipleindependentvariablesinourcase.


I. INTRODUCTION
R EMOTE sensing systems have exceedingly advanced precision agriculture applications such as growth stage classification, yield estimation, and harvest scheduling in the past few decades [1]- [3]. Among variables that affect the mentioned practices to a large extent is leaf area index (LAI) [4]- [7]. LAI is commonly defined as one-half of the total green leaf area per unit horizontal ground surface area [8], [9]. It is a critical variable that governs multiple canopy-light processes and essentially quantifies the amount of photosynthetic area in an ecosystem [10]. Several LAI measurement devices have been developed for LAI field measurements based on Beer-Lambert's law; these include digital cover photography, LAI-2200 (LI-COR Inc., Lincoln, NE, USA), and AccuPAR LP-80 ceptometer (Decagon Devices Inc., Pullman, WA, USA). However, the use of such devices results in a significant labor and time cost, especially for wide crop fields.
Over the past few decades, satellite-and airborne-based remote sensing methods thrived and have been used to estimate LAI across coarser spatial scales. Existing techniques can be divided into two categories, based on how the data are collected. The first is termed passive remote sensing and is based on spectral imagery, such as color (RGB), near infrared (NIR), multispectral, and hyperspectral sensors. The second encompasses active remote sensing and typically is based on light detection and ranging (LiDAR) or synthetic aperture radar systems [11]. While promising results have been demonstrated in many studies, the temporal and spatial resolutions are still lacking when the plot-scale evaluation is needed. It is in this context that the unmanned aerial systems (UAS) based remote sensing methods have proliferated as a value-adding option to satellite or airborne-based methods. Both passive sensors and active sensors can be mounted on a UAS and simultaneously capture different types of data for crops. There are two main types of UAS-based methods for LAI assessment: spectral features derived from reflectance imagery, and structural features derived from the three-dimensional (3-D) point clouds, i.e., LiDAR or structure-from-motion (SfM) based point clouds.
As for the models derived from spectral features, Yao et al. [12] used UAS-based narrowband multispectral imagery to estimate wheat LAI at the middle-to-high LAI levels. The authors concluded that the LAI model, based on the modified triangular vegetation index (MTVI2), provided the highest coefficient of determination (R 2 ≈ 0.8) and lowest normalized root-mean-squared error (nRMSE = 24%). Another study [13], in turn, applied UAS-based multispectral camera imagery for sorghum phenotypic assessment and determined that the normalized difference vegetation index (NDVI) and the enhanced vegetation index (EVI) correlated well with LAI, both during the vegetative growth phase (preanthesis) and at maximum canopy cover shortly after anthesis (R 2 = 0.66 − 0.85). Researchers from a related effort [14] also developed a regression model for evaluating LAI from NDVI with a significant correlation of R 2 = 0.77 via a UAS-based multispectral camera. Such spectral approaches to LAI assessment often have been augmented by 3-D sensing inputs, i.e., 3-D point clouds generated by LiDAR or SfM, given the ability of the latter to directly measure plant structure.
As for structural feature-based methods, Comba et al. [15] used UAS-LiDAR data to explore the effect of leaf occlusion on LAI inversion for maize crops. They determined the optimal voxel size for inverting LAI and demonstrated better results by designing a flight direction that is perpendicular to the maize rows. A related study [16] used a UAS-LiDAR system to evaluate structural parameters of a mixed-species restoration plantation experiment. Three structural variables were analyzed from LiDAR data, including canopy height, gap fraction, and LAI, and a significant correlation (R 2 = 0.84; nRMSE = 15.5%) between LiDAR-derived canopy height/LAI and field-observed aboveground biomass was found. A related study [14] used UAS-based SfM point clouds to visualize and quantify vineyard canopy LAI, with a moderate R 2 value of 0.567 being reported. Finally, Comba et al. [15] evaluated the LAI of 704 vines in a vineyard using SfM point clouds from UAS imagery and obtained an solid R 2 value of 0.82. The performance metrics for the mentioned studies might sound promising, but they were constrained by various factors.
The most well-studied environments for LAI and UAS-based methods are forests or relatively tall crops, such as maize, wheat, and grapevines [12]- [18] (see Table I). It is worth noting that for smaller/shorter crop types, many methods for extracting structural parameters failed [11], [19]. To our knowledge, only a few studies have focused on point cloud-based methods that are applied to short (<0.5 m) crops. It has become evident that unique problems emerge in such scenarios. While the vegetation canopy becomes smaller/shorter, the associated biophysical features become smaller, and noise from the background material could cause challenges [20]. Vegetation points near the ground are difficult to separate from ground points. Mild disturbance such as weeds, rocks, depressions, and even windy weather could add significant noise to analyses. In other words, difficulty in analysis increases as the "signal-to-noise ratio" decreases [7], [21].
These complications require a robust approach that could take into account various characteristics of the environment being studied. One such method is data fusion processes. Recent studies have shown that data fusion, which utilizes data from different sensor platforms, often can achieve more accurate results. A related study [32] found that the fusion of RGB image and thermal image could improve the evaluation of LAI in a soybean field by reducing the nRMSE by 0.4% when compared to only using a RGB camera. In a later study [33], the authors showed that fusion of high resolution RGB, multispectral, and thermal data could provide an accurate estimation of biochemical and biophysical parameters and improve the yield prediction accuracy from a R 2 value of between 0.26 and 0.52 for a single imaging sensor, to 0.67 for multiple sensors. Sankey et al. (2017) [34], in turn, demonstrated that LiDARhyperspectral image fusion performed more accurately (88% overall accuracy) than either data source alone for representing a gradient of vegetation and topography in northern Arizona, USA. In a follow-up study in arid and semi-arid land vegetation monitoring [35], the fusion of LiDAR-derived plant height estimates and hyperspectral images derived spectral signatures resulted in an overall accuracy of 84-89% for vegetation species classification, which outperformed either data type alone.
Our primary goal in this article, therefore, was to estimate the LAI of a short-crop, namely snap bean (Phaseolus vulgaris L.), using UAS-based SfM and LiDAR point clouds, as well as the spectral information from multispectral imagery. We combined output products from both the structural (3-D) and spectral datasets to improve their preprocessing pipelines. LAI predictive models then were derived for both modalities. We pursued three main objectives: 1) assess whether LiDAR-derived and multispectral imagery-derived SfM models can provide reliable estimates of LAI, 2) identify the most suitable models, either as a single or fused approach, for structural characterization of snap bean crops, and 3) compare the LiDAR-derived and multispectral imagery-derived SfM models and explore the potential of using the more cost-effective SfM point clouds to replace LiDAR in precision agriculture applications, using snap bean as a proxy crop.

A. Study Site and Data Collection
The experimental field is located in Geneva, NY, USA (42°49 53 N, 77°00 50 W; Fig. 1) with an area of 2060 m 2 (20 m × 103 m). The field has a ∼ 2 • slope from the westto-east. Four replications of six different snap bean cultivars, namely Venture, Huntington, Colter, Cabot, Flavor Sweet, and Blevet were planted, resulting in 24 plots, with four rows for each plot. We used 11 flights to collect data over six days at different growth stages, from July 8-August 20, 2019. There were two flights at different altitudes (∼28 m and ∼51 m) each day, except for the single flight on July 8, which was for generating a digital elevation model (DEM) when the field only contained snap bean seedlings, i.e., the bare soil was mostly visible. The selection of subsequent dates largely followed snap bean phenological events: August 5 was selected because the snap beans were blooming and changing rapidly; on August 12, the snap beans were close to maturity, and the blooming period had ended; on August 14, the snap beans were almost ready for harvest and pods were relatively mature; and finally, some beans were harvested on August 16 and August 20, i.e., these dates represented fully mature, harvest stages.
The UAS-based mapping system consists of a DJI Matrice 600 Pro hexacopter, which carries a Global Navigation Satellite System (GNSS)/inertial measurement unit (IMU) unit, a Velodyne VLP-16 PuckTM (Velodyne, San Jose, CA, USA) LiDAR, and a MicaSense RedEdgeTM (Micasense, Seattle, WA, USA) multispectral camera. The GNSS/IMU unit records geolocation and GPS time during flights. The VLP-16 Puck generates up to ∼600 000 points/second in dual return mode. It holds a ±15°v ertical field-of-view and a 360°horizontal field-of-view. The LiDAR's laser wavelength is 903 nm, and the range accuracy is ±3 cm [36]. The MicaSense RedEdge captures imagery from five discrete spectral bands: blue, green, red, red edge, and near infrared, centered at 475 nm, 560 nm, 668 nm, 717 nm, and   [37]. Sampling measurements within each plot were then averaged to be a representative of the plot on the date. Therefore, for each flight on each date, we have 24 samples of ground truth canopy height and LAI. The statistics of the field measured LAI is shown in Table III.
We used every two adjacent crop rows and the between-row gap as an elementary sampling unit (ESU) [10] when processing the LiDAR point clouds, which corresponded to the field samples. We calculated the predicted LAI of each ESU and then averaged the three ESUs in each plot to obtain the predicted LAI of the plot.
It is worth noting that the VLP-16 is a widely used industriallevel LiDAR. The system has a relatively larger beam divergence (a horizontal beam divergence of 3.0 mrad and a vertical beam divergence of 1.5 mrad), when compared to survey-level Li-DAR systems, but it comes at a lower acquisition and operational cost [38]. According to the descriptions in the VLP-16 manual [26], the footprint size would be 150 × 85 mm in flights at 50 m altitude and 81 × 47 mm for flights at a 25 m altitude, while the snap bean leaves can grow 6-15 cm long and 3-11 cm wide. We intentionally selected the Velodyne VLP-16 as our LiDAR modality in the hope of eventual operational implementation of our methods. The AccuPAR LP-80, on the other hand, boasts a proven track record as an effective in situ LAI measurement device, as evidenced by a variety of previous studies [10], [37], [39]- [41]. These considerations led to the selection of our respective airborne and ground-based 3-D assessment tools.

B. LiDAR and MSI Data Preprocessing
The raw LiDAR data were provided directly by our data collection team, with each LiDAR point containing data/header fields for x, y, z, intensity, return number, and the GPS timestamp. The multispectral imagery was processed using the Pix4DMapper (V.4.4.12) for generating SfM point clouds, image mosaics of the crop field, and vegetation indices. This section will describe details of both datasets.
1) Registering MSI With LiDAR: SfM point clouds and multispectral mosaic images were significantly shifted and inclined from LiDAR point clouds due to differences in platform design and flight collection parameters, making it challenging to utilize a fusion of the LiDAR and multispectral data directly. The spatial disparities were attributed to the following issues: 1) the ground control points (GCPs) were not well utilized for mosaicking images or generating the SfM point cloud; 2) different horizontal and vertical georeference coordinate systems were used for different sensors; and 3) differences between the ellipsoidal height and the orthometric height. Therefore, we registered the multispectral data with LiDAR data via an independent step.
Six valid GCPs (AeroPoints) were used to optimize the process of mosaicking images and generating SfM point clouds. The GCPs used NAD83 as the geographic coordinate reference system (horizontal) datum and NAVD88 as the vertical coordinate reference system datum. The LiDAR point cloud accordingly used WGS84 and EGM1996. Different reference coordinates resulted in several meters difference of the x/y coordinates. Second, the default z value in SfM point clouds corresponded to orthometric height, i.e., the distance of a point on the earth's surface to the reference geoid, which is determined by the earth's gravity and approximated by the global mean sea level. In contrast, the z value in LiDAR point clouds was ellipsoidal height, meaning the distance of a point on the earth's surface to the ellipsoid that approximates the earth's surface [42]. The difference between the ellipsoidal height and the orthometric height is called "geoid-ellipsoid separation" [43], calculated by where H o is the orthometric height, H e is the ellipsoidal height, and N is a signed number [44]. At our experimental field, N = −34.9 m.
We, therefore, first converted the coordinates of the GCPs from NAD83/NAVD88 to WGS84/EGM1996 by using the free software tool VDatum 4.0.1 [45]. We then used the transformed GCPs to reoptimize the SfM point clouds. This was followed by calculating geoid height at our experimental field and shifting the SfM point clouds to match the LiDAR point clouds vertically. Fig. 2 shows examples of the cross-sections of the 3-D point cloud from the ninth flight; note that each ridge represents part of a row crop. The valleys in the middle of each two adjacent ridges represent the between-row space. Judging from the peaks and valleys, we found that the LiDAR and SfM point clouds matched well in all 3-D. The LiDAR point clouds were generally higher than the SfM point clouds and exhibited more variation, especially at the top of the canopy rows. We attributed this to the fact that LiDAR operates by actively recording the reflected signal from objects in space, while the SfM algorithm is based on feature extraction. As the leaves at the top of canopy were small and erectophile in nature, it was challenging for the SfM approach to capture top leaves.
2) LiDAR Data Preprocessing: LiDAR sensing enables probing of both the horizontal and vertical structures of crops by actively emitting high-frequency laser pulses toward the object and recording the reflected responses as a function of time [46], [47]. While this design enables such active systems to capture more in-depth structural information than classical passive imagery sensors, it also creates specific challenges, such as complex data preprocessing, noise points, nonuniformity across different scanning angles, and a rapid increase in footprint size as the detection range increases. Moreover, different vegetation characteristics, in terms of height, layering, and leaf characteristics, limit the effectiveness and generality of evaluation models based on LiDAR-derived metrics. Therefore, a robust assessment requires proper preprocessing of the raw point cloud and careful selection and validation of predictive models.
We implemented an amended version of the UAS-based Li-DAR data preprocessing pipeline described in [48] (see Fig. 3). First, we retrieved the flight trajectory from IMU recordings. We then cropped the point clouds according to the valid spatial boundaries of x/y/z and temporal boundaries of GPS time, using the las2las function embedded in the LAStools software [49]. Subsequently, the outliers in the point cloud were removed using the statistical outlier removal filter, which computes the average distance of each point to its k nearest neighbors and then eliminates the points that have larger distance to their neighbors than the average distance plus n σ times the standard deviation [50]. Then, the second return points were removed. The LiDAR operated in dual-return mode, and thus, could distinguish two returns on a per-pulse basis only if the distance between the two returns was greater than 1 m [36]. However, since snap bean plants only grow 0.3-0.6 m tall, the first return points were adequate for representing most of the crop-specific structural information for the LiDAR data. The next step involved the removal of duplicate points, which were ≤0.001 m away and may result during preprocessing of raw data. These duplicates could inflate the pulse density and consume additional disk space [51].
Next, we retrieved the scan angle of each 3-D point in the point cloud by associating it with a flight trajectory point with the nearest timestamp. Thus, from (1), we calculated the horizontal and vertical scan angles α 1 and α 2 of each LiDAR point where x l , y l , z l , are the coordinates of the LiDAR point, and x t , y t , z t are the coordinates of the flight trajectory points. Points with large scan angles were removed because 1) the point density decreases rapidly at large angles, and 2) measurement error could increase substantially as the laser beams propagate further (the footprint size of a single laser beam increases by a factor of two as the propagation distance increases). We, therefore, set a threshold of ±20 • to only retain points with small enough vertical and horizontal scan angles. While this threshold was ±22.5 • in [48], we considered that this threshold was adjustable for different data characteristics and objectives. The flight line overlap percentage can also impact the selection of this threshold. We then generated a DEM to create the normalized height to ground dataset by 1) applying the cloth simulation filter (CSF) [52]- [55] via CloudCompare (version 2.11) to identify pure ground points from the point cloud collected on July 8, when the bare soil was mostly visible, 2) subsampling the ground points to a 0.02 m grid, given that the ground points are even denser than the target grid data, and 3) ordinary kriging interpolation [56] to fill the empty cells. Then, we normalized all z coordinates from the other point clouds to heights above ground using the DEM [26]. Most LiDAR-derived predictive models rely heavily on tallying the number of points or feature values of all points per class (vegetation versus ground). Therefore, the final models' efficacy and accuracy largely depend on the quality of the segmentation of ground points versus vegetation points. Many studies for vineyard or forestry applications [14], [57], [58] used a manually selected threshold of normalized height, which could be 0 or a small positive number, to achieve the segmentation. Points lower than the threshold are classified as ground points, and points above the threshold are classified as nonground points. However, crops like snap bean have a much smaller height than vineyards (meters) or trees (tens of meters), ending up with severe blending of near-ground vegetation points and ground points. A random threshold could easily bring bias from either ground side or vegetation side. We, thus, resorted to the segmentation results from 2-D image mosaics of the whole field. By combining the vertical traits from LiDAR and the horizontal textures, we hope to minimize ground versus nonground return segmentation error. We used the ENVI 5.5.2 to create a 2-D segmentation mask for each multispectral image. First, we generated mosaic images of the whole field, corresponding to the five bands, using the Pix4DMapper (V. 4.4.12). Second, we stacked the five single-band mosaic images and then applied the spectral angle mapper (SAM) [59] to classify vegetation pixels versus ground pixels. SAM operates by comparing the spectra of each pixel to known reference spectra or endmembers. It is insensitive to illumination change, since only the vector direction is used and not the vector length. Then, we used the "mask" to assess results from a list of different threshold valuesẐ th1 ,Ẑ th2 , . . . ,Ẑ thN . From the 3-D segmentation results, corresponding to eachẐ thi , we used stratified sampling on both ground and vegetation points according to their x/y values and matched the sampled points to the nearest pixels in the 2-D mask. A contingency table, thus, could be constructed.
Evaluation parameters such as recall, precision, total accuracy, and Cohen's Kappa coefficient (κ) [60], [61] were calculated (see Table IV). Since we want the highest accuracy for both ground points and vegetation points, we determined the bestẐ th by selecting the threshold where the highest Cohen's Kappa occurred, i.e., the classification results of the 2-D SAM and of the z-threshold-based method reached the maximum overall agreement. We then applied the bestẐ th as a threshold for differentiating ground and nonground points in the LiDAR point cloud. Fig. 4 displays an example of determining the idealẐ th for the LiDAR point cloud and the values of differentẐ th for LiDAR point clouds from different flights. The idealẐ th ranged from 0.06-0.14 m for different flights. Fig. 5 shows an example of a top-down view of the segmented nonground points. Apart from a few points that were 1 m above the ground (in-field plate markers) and a small portion of misclassified ground points, most of the points were within 0.1-0.4 m above ground, i.e., resulting from the snap bean plants. The segmented point clouds then were separated into ESUs by 1) manually creating polygons of the ESUs' boundaries, 2) intersecting the ground points and the nonground points of each ESU with its associated coverage polygon, and 3) allocating the points to different groups by their polygon ID. Points for each ESU were used to derive one observation (x i ) for evaluating LAI from LiDAR. Our next step involved the addition of spectral information to augment the 3-D datasets.
3) MSI Data Preprocessing: Multispectral imagery can provide rich information about the crop from two general perspectives. The first involves the generation of a 3-D point cloud of the field through a photogrammetric range imaging technique, i.e., SfM [62]- [64]. The other focuses on spectral information, i.e., deriving vegetation indices from the multiple bands per image [65]- [68]. We implemented and combined both methods in this section to generate extended SfM point clouds that contain both structural (3-D) features from the x/y/z coordinates and spectral features from vegetation indices. Fig. 6 displays the processing workflow for multispectral imagery. First, all multispectral images from each flight were fed into Pix4Dmapper to generate an SfM 3-D point cloud and five mosaic reflectance images (one for each band) of the whole field. We next separated the vegetation points from nonvegetation points via the following steps: 1) classifying the  2-D mosaic image, stacked from five single band mosaic images, into vegetation versus nonvegetation pixels using the SAM classification in ENVI and 2) "masking" the 3-D points according to the classification labels and the geographical information (x/y coordinates) in the 2-D mosaic image. We then normalizedthe the z coordinates in the SfM point cloud to height-above-ground and derived normal vectors for vegetation points; additionally, we derived vegetation indices from the five mosaic images. Finally, by using QGIS (version 3.12.2), we intersected the 3-D vegetation points with the mosaic images and polygons of the ESU boundaries. Thus, we integrated the structural and spectral information, as well as the ESU ID, into the final point cloud.
One possible caveat of this 2D-classification-to-3Dsegmentation method might be that some ground points directly beneath vegetation points could be mistaken as vegetation points. However, we considered the impact of this drawback as negligible in our project. Given the fact that the images were collected at a nadir angle, the dense snap bean canopy hindered the camera from capturing the ground underneath leaves, which resulted in the SfM point cloud only having one "visible" layer for the imagery in the study area.
We mainly utilized structural features that consisted of normalized heights and projected normals. The z coordinates in the SfM point clouds were normalized to height above the LiDAR-derived DEM. We used CloudCompare to calculate the point normals. We hypothesized that the normal vectors could explain a plant canopy's inherent structural variability. In other words, normals that were uniformly distributed across a plant canopy were indicative of a closed, even canopy, not dissimilar to an umbrella's form. In contrast, highly variable normal angles represent a gap-filled, uneven plant canopy. CloudCompare outputs normals as dip angle (zenith) and dip direction angle (azimuth). These values were calculated using Dip angle = arccos n z n Dip direction angle = arctan n x n y where n refers to the unit normal at a certain point, and n x , n y , n z are components of the normal. A typical vegetation spectrum shows absorption features in the red and blue wavelengths, slight reflectance in the green wavelength range, and a strong NIR reflectance plateau, with water absorption features at the longer shortwave-infrared wavelengths. Vegetation indices (VIs) are broadly used to reflect and enhance the relationships between these spectral absorptionreflectance spectral features and, thus, provide meaningful information about plant growth stage, plant health, and among other important characteristics [69]. Considering the five bands in our data, we selected five VIs by referencing previous relevant studies (see Table V). They are the normalized difference vegetation index (NDVI), simple ratio index (SR) Anthocyanin Reflectance Index 2 (ARI2), green NDVI (GNDVI), and normalized difference red edge index (NDRE).

C. Predicting LAI From LiDAR and MSI
1) LiDAR LAI Models: Multiple models have been proposed to retrieve LAI from LiDAR point clouds, with [10] providing an extensive overview. The models can be divided into three main categories: 1) based on Beer-Lambert's law, LAI is correlated with the gap fraction, which usually is derived through various laser-based metrics, such as the laser penetration index (LPI) [72], [73] and the all echo cover index (ACI) [74]; 2) LAI is linearly correlated with laser-based metrics, in which the correlation coefficients are determined by regression [57], [58], [75], and 3) LAI is evaluated through its allometric relationship with other LiDAR-derived biophysical parameters, such as vegetation height [72] and canopy cover [76].
We selected eight models to test our data (see Table VI). In model 1-4, we used two different methods to determine the β value: 1) simply assuming a spherical leaf angle distribution, which gave us β = 2 for all plots [77]; and 2) calculate β via its relationship with the extinction coefficient k [78], i.e., β = 1 k , where k can be decided by the solar zenith angle and the leaf angle distribution parameter χ. The zenith angle in our study was retrieved from time and geo-coordinates. The leaf distribution parameter χ refers to the ratio of the length of the horizontal semi-axis to the vertical semi-axis of an ellipsoid, described by the leaf angle distribution function of the canopy [78], [79]. The default value for χ is 1.0 for the AccuPAR LP80, which assumes a spherical canopy angle distribution. While the value of LAI is not strongly dependent on the value of χ, we adjusted this value as the ratio of the measured row width to the measured canopy height, assuming that the wider a snap bean plant becomes, the more horizontal its leaves will be. In models 5-8, a and b were empirical parameters derived from linear regression between the medium terms (p or CH) and the measured LAI.
2) MSI Evaluation Models: We used both the structural and spectral descriptors per plot to derive multiple statistical metrics, including the mean (μ), standard deviation (std), coefficient of variation (CoV), first quartile (q1), median, third quartile (q3), interquartile range (IQR), and trimmed mean (trMean). These metrics (7 descriptors × 8 metrics = 56 predictors) then were used to predict LAI through multivariate regression. We input these independent variables to three feature selection and regression algorithms using a widely used classical algorithmstepwise regression [57], [58], two common machine learning algorithms-Lasso cross-validation (LassoCV) [80], and recursive feature elimination (RFE), which used linear support vector regression (SVR) [81] as the estimator. For stepwise regression, we used backward selection when the sample size was larger than the number of features (n>p) and forward selection in the opposite case (n<p) [82]. We first used LassoCV and RFE-SVR to select the most significant features, after which collinear variables were removed from the selected subset, if their variance inflation factor > 5-10 [83], [84]. Multivariate linear regression subsequently was implemented for the remaining features to develop prediction models for LAI. All three methods were implemented in Python 3, and LassoCV and RFE-SVR were implemented by using the scikit-learn package [85]. The software and packages used in this project along with their objectives can be found in Table IX.

3) Modeling From Different Combinations of Datasets:
In both datasets, we first derived models from all the data from flights 2-10. Field measurements on each day were reused as corresponding true responses for both flights. Next, we assumed that different snap bean cultivars might require different evaluation models, since their physical characteristics varied as they matured. The six cultivars were separated into three groups according to their morphological characteristics, i.e., Venture and Huntington belong to the "Large cultivar variety (L)" group; Colter and Cabot belong to the "Four sieve cultivar (F)," and Flavor sweet and Blevet belong to the "Whole sieve cultivar (W)" [86].
Flight altitude could also affect the point density in LiDAR point clouds and the laser pulse average footprint size. Denser points and smaller footprint size arguably could lead to more structural crop detail. Lower flight altitudes, in general, lead to higher LiDAR point density and smaller footprint sizes. We, therefore, also assumed that flight altitude could impact the evaluation results. We evaluated two sets of prediction models to test this assumption for flights at two different altitudes, i.e., the 3rd, 7th, 9th, 10th flights as a group at 28 m altitude and the 2nd, 4th, 6th, 8th, 11th flights as the other group at 51 m altitude.

D. Evaluation Metrics
We used the coefficient of determination (R 2 ) the root-meansquared error (RMSE), and the normalized (nRMSE) to evaluate the accuracy and precision of the predicting models from LiDAR and MSI. Moreover, since MSI models were derived from multivariate regression, we calculated the adjusted R 2 [14], [90], the predicted R 2 [91]- [93]. The adjusted R 2 essentially is the R 2 that has been adjusted to compare models for different numbers of predictors. It increases only if a predictor improves a model by more than what is expected by chance [90]. The predicted R 2 reflects the ability of the regression model in predicting the response of a new observation and reflects how much of the I canopy refers to the sum of points' intensity values. R ≥T refers to the number of returns above or equal to the height threshold, i.e., the idealẐ th in Fig. 4(b) and R <T refers to the number of returns below the height threshold. R total is the number of total returns. Res 0. 95 and Res 0.05 refer to the 95% quantile and the 5% quantile, respectively, of the height above the DEM. variance in the dependent variable is explained by the model. If a model is overfitting, the predicted R 2 will drop significantly. The equations for calculating five evaluation metrics are as follows.
where y i is the true response and y i is the predicted response for the i th observation.
, where n is the number of observations.
where N is the number of observations, p is the number of independent variables in the model. R 2 is calculated from 1). 4) P red.
, where SS total is the total sum of squares and PRESS is the predicted residual error sum of squares (PRESS) statistic. PRESS is essentially a form of leave-one-out crossvalidation: after removing one observation each time, the model is refitted from the remaining observations and then predict response to the removed observation (ŷ i, −i ), then the PRESS statistic is calculated as the sum of square of the resulting error [94], [95].

III. RESULTS
We derived multiple models for predicting LAI from both LiDAR and extended SfM point clouds. By providing the evaluation metrics mentioned previously, we compared not only models for the same method or under the same condition, but also models across different methods and conditions. Fig. 7 shows the scatter plots of LiDAR-predicted LAI versus field-measured LAI and the R 2 , RMSE, and nRMSEof the model 1-8 in Table VI. Except for model 8, which exhibited some correlation between the predicted LAI and measured LAI with R 2 = 0.05, all the other models provided moderate accuracy with R 2 values ranging between 0.55-0.61, RMSE ranging between 0.38-0.59 m 2 /m 2 , and nRMSE ranging between 19% and 31%. The model based on LPI achieved the highest R 2 of 0.61 and the lowest RMSE of 0.38 m 2 /m 2 .

A. LiDAR-Based LAI Estimation
Figs. 8 and 9 list the R 2 and RMSE of the models generated from data on a percultivar and percultivar group basis, respectively. The models were sensitive to the differences in cultivars and within cultivar groups.
Models for the Flavor Sweet cultivar tended to yield the highest R 2 and lowest RMSE, while models for the Venture cultivar resulted in the lowest R 2 and the highest RMSE. The highest R 2 = 0.89 was achieved by ACI, β = 1 k model on the Flavor Sweet. The lowest RMSE = 0.26 and nRMSE = 0.14 were achieved by the f c model and √ LPI model for the Flavor Sweet cultivar. Models based on ABRI, f c , and LPI (5, 6, 7) provided consistently better predictions for nearly all cultivars. The CH models did not exhibit apparent correlations between the predicted LAI and the measured LAI for most cultivars, and were inferior to other models. However, for the Flavor Sweet cultivar, it gave satisfactory predictions with R 2 = 0.82 and an RMSE = 0.31. We attributed this difference to Flavor Sweet's strong resistance to lodging, even during harvest time.
In Fig. 9, similar to the percultivar-polylines, the models based on ABRI, f c , and LPI continued to perform better than other models. The CH models again showed no predictive value (close to zero R 2 for cultivar group L and F), but moderate accuracy and precision for cultivar group W (R 2 = 0.51 and RMSE = 0.45). Fig. 10 shows a comparison of the evaluation metrics for all the models derived for two flight altitude groups. Generally, all the models from the 28 m flights yielded better results than  Table VI. their corresponding models from the 51 m altitude flights, except for the CH models. Table VIII shows the results of the models for evaluating LAI using stepwise, LassoCV, and RFE-SVR methods. For the stepwise regression, we set α = 0.05 as the include (forward selection) or exclude (backward selection) threshold [14], [96]. For the LassoCV and the RFE-SVR methods, we used the default setting of key parameters in the scikit-learn package. Three regression/feature selection methods achieved similar results based on the adjusted R 2 , predicted R 2 , and nRMSE. For the models built from data from all flights, the adjusted R 2 and the predicted R 2 were 0.48-0.53, and nRMSE values between 0.21 and 0.22. The derived model using the LassoCV method only Fig. 8. LiDAR regression results per cultivar. Note the apparent differences between results for the various cultivars. All the models achieved similar results, except for the CH-based model. Fig. 9. LiDAR regression results per cultivar group. Across all models, the whole bean cultivar group resulted in higher accuracy than the other two groups.   11. Adjusted R 2 , predicted R 2 , and nRMSE of the L, F, and W cultivar group. According to Green's rule-of-thumb [98] for samples per variable, since the sample size decreased to 80 for each model, only a maximum of four predictor variables were included in each model.  [98] required three variables to achieve nearly the same results as the stepwise and RFE-SVR models. The third quantile of NDVI and the standard deviation of NDRE were shared variables across all three models. From Fig. 10, we found that the 28 m flight group achieved obviously better results than the 51 m flights, with predicted R 2 = 0.72 to 0.73 and nRMSE = 15% to 16%. The variables shared between models also changed. For the 28 m flight group, these were the standard deviation of GNDVI and median of NDRE and for the 51 m flight group, the common variable was the standard deviation of the dip angles. Fig. 11 shows the results of MSI models derived from cultivar groups. We found that the models for the W group achieved relatively higher adjusted and predicted R 2 values than models in the other two groups, although their nRMSE were similar. We did not derive regression models for flights per cultivar as in the LiDAR data analysis, since Green's rule-of-thumb states that when selecting samples per variable, the sample size should not be smaller than 50 in order to avoid overfitting in multivariate regression [97], [98].

B. MSI-Based LAI Estimation
Moreover, minimum sample size requirements increase with an increased number of predictor variables (see Table VII). The "all flights" group contains all the observations for the 24 plots during flight numbers 2-11, i.e., 240 observations. When they were split into a 28 m flight group and a 51 m flight group, the sample sizes decreased to 96 (flight number 3/7/9/11) and 144 (flight number 2/4/5/6/8/10), respectively. It followed that these models would be restricted to up to six predictors. However, when the observations were split into six cultivars, each sample only had 40 observations, which violates the rule-of-thumb.

A. LiDAR-Derived Metrics Modeling
Among the LiDAR-derived models, the models based on ABRI f c , and LPI consistently provided the strongest correlation and the smallest error. This result corroborated the findings from studies in [57] and [58]. We observed that the LiDAR-derived models (models 1-7) generally account for around 60% of the variance in measured LAI (see Table VI and Fig. 7). The models based on Beer-Lambert or the linear/logarithmic regression achieved similar R 2 , and the latter ones yielded relatively lower RMSE values. The model based on an allometric relationship could not explain the variance adequately. We attributed this to the model only considering one structural metric, i.e., the canopy height, while no other structural information was included. If additional structural metrics, such as plant width and aboveground biomass were measured, this model's accuracy and precision likely could be improved, as mentioned in the work of [99].
Two phenomena seemed counterintuitive in Fig. 7-the comparison between models 1 and 3 with 2 and 4, showed that a more precisely retrieved β from row width and canopy heights did not perform better than the naive assumption of β = 2. This could be attributed to the uncertainty from field measurements [100] and systematic inconsistency between different flights. By comparing models 1 and 2 with 3 and 4, we found the sum of intensity ratio (IRI) models yielded slightly lower R 2 and higher RMSE values than the point number ratio (ACI) models, which implies that integrating point intensity in LAI evaluation did not improve the evaluation accuracy. This result could be due to the radiometric conditions changing across different flights, which in turn hints at the need for additional research into the impact of illumination conditions on LiDAR-based models.
Results listed in Figs. 7 and 8 showed that 1) models for the Venture cultivar exhibited worse predictions than the general TABLE VIII LAI PREDICTIVE MODELS GENERATED FROM THE EXTENDED SFM POINT CLOUD DATA * Features were normalized by subtracting the mean and divided by standard deviation before modeling. * * Since the sample size for the 28 m flight group models was only 96, according to Table VII, to avoid overfitting, we cannot legitimately use more than six independent variables. Therefore, the number of features was limited to six by limiting the number of best predictor features. Note: The naming convention for model variables follow a format of "[abbreviation of statistical metric]_[structural or spectral feature]." For example, μ_h implies the mean of the height above DEM, CoV_dip refers to the coefficient of variation of the dip angle of normals, and q3_NDRE means the third quantile of NDRE. models for the whole field (highest R 2 < 0.40, lowest RMSE > 0.40); 2) models for the Cabot cultivar yielded nearly the same prediction performance; and 3) models for the other four cultivars, i.e., Huntington, Colter, Flavor Sweet, and Blevet achieved better results (highest R 2 > 0.65, lowest RMSE < 0.30). This was attributed to the morphological (biophysical) characteristics that varied across different cultivars. For example, in our observation, leaves of cv. Venture expanded parallel to the ground with maturity. The canopy of other snap bean cultivars began to lodge with crop maturity. In contrast, snap bean cv. Huntington (also a member of the L group), had leaves, which grew closer to the main stem and the entire canopy was upright and remained so until maturity. Cultivars with a tendency to lodge as crops matured tended to affect the accuracy and precision of the LiDAR models. Cultivars subject to lodging may be strongly affected by weather such as heavy rainfall. Since data were collected over the full growing stages of the snap bean crop, the tendency to lodge along with the influence of external variables likely impacted modeling results across days. These external factors could help explain the variation in model results for the three cultivar groups in Fig. 9. We concluded that leaf distribution angles for different cultivars may contribute to the observed differences, and therefore, we suggest that future studies include an assessment of leaf angles. The models for the whole bean cultivar (W) group achieved relatively better predictions than the other two groups, with the largest R 2 = 0.79 and the lowest RMSE = 0.30.
Additionally, when comparing Figs. 7 and 10, we observed a general relationship among prediction accuracy and precision for models 1 to 7 at different flight altitude groups, where the 28 m flight group > all-flight-group > 51 m flight group. This observation supports our assumption that flight altitude is a significant factor in predicting LAI from UAS-based LIDAR data. The lower flight altitude enhances the capture of structural detail of the crop by generating LiDAR point clouds of higher density and smaller footprint size.

B. MSI-Derived Metrics Modeling
The models derived from the three feature selection methods (see Table VIII) exhibited similar accuracy and precision (Adj. R 2 , P red. R 2 , and nRMSE values), demonstrating the robustness of our data processing methods. While we observed one or two common features among the models from different methods, it is noticeable that the models apparently varied based on selected features. We attributed this to the fact that feature selection methods use different metrics to eliminate less important features. More specifically, stepwise regression uses the p-value as its criterion, while LassoCV retains only nonzero estimated coefficients in its sparse solution in L1 norm, and RFE-SVR prunes features by ranking the coefficients of features and recursively considering smaller and smaller sets of features.
Moreover, as shown in Fig. 12, features from the different models were correlated. For example, RFE-SVR selected the covariance of ARI2 (square A) and median of ARI2 (square B) as the two most significant features. Although these two features did not show up in the other two models within the same group, they are significantly correlated with the covariance of SR (stepwise model) and the third quartile of NDVI (LassoCV model), respectively. We, therefore, concluded that although different methods selected different features, many of these  Table VI. Yellow lines represent the boundaries of features from the same models. Note that pixels at A and B exhibited a significant correlation. features in fact were correlated. Such correlation is to be expected for a dataset of 96 samples × 56 features. If our methods were to be extended to larger datasets, we speculate that more common features will be found among the models, and the between-feature correlation will be reduced. Judging from the performance and the number of features of the models, we recommend the LassoCV-based models based on performance and the required number of predictor variables.

C. LiDAR Versus MSI
For the general models across all flights and all cultivars, the best LiDAR-derived model achieved slightly better prediction results than the best multispectral imagery-derived model, with R 2 lidar = 0.61 versus R 2 msi = 0.53 and nRMSE lidar = 19% versus nRMSE msi = 22%. However, when we separated the data according to their flight altitudes, the best MSI-derived models surpassed the best LiDAR-derived models. For the 28 m flight group, we observed R 2 msi = 0.75 > R 2 lidar = 0.65, and nRMSE msi = 15% < nRMSE lidar = 19%; for the 51 m flight group, we observed R 2 msi = 0.54 > R 2 lidar = 0.50, and nRMSE msi = 21% < nRMSE lidar = 22%. This provided credence to our claim that an SfM approach, based on multispectral imagery, could be a competent alternative to an costly LiDAR sensor for measurement of LAI of snap beans under the same flight settings. However, it is noticeable that the MSI-derived models required more independent variables than the LiDAR-derived models. Such complexity in models could undermine their robustness and generality when research expands to alternative crops. While some studies have shown that the performance of different UAS-LiDAR systems varies under different settings, such as flying height and speed [101], we acknowledge that the conclusion was based on a set of limited conditions and therefore should not be generalized too broadly.
The methods based on LiDAR data and multispectral imagery shared some common findings: 1) the 28 m flight group generally performed better than 51 m flight group, which attributed to a denser spatial sampling of crop structure at lower altitudes; and 2) the models for the W cultivar group exhibited higher prediction accuracy than models for the L and the F cultivar group, which was attributed to more consistent growth patterns within a crop and during the season. These findings could offer guidance for flight planning when predicting LAI of other crops. First, without causing undue impacts on the operational flight execution, it is preferable to operate the UAS at a lower altitude if the environment and UAS battery power allow. Second, crops' biophysical variance characteristics could result in additional complexity (and poorer performing models) when assessing LAI remotely. The more uniform and regular the crop growth patterns are, i.e., leaf angles and leaf direction, the more likely an accurate and precise LAI modeling outcome.
It is informative that despite the challenges associated with the relatively low canopy height and high structural/temporal growth variance in snap bean crops, our model performance was comparable to other recently published studies. A related study [12] evaluated LAI on winter wheat. Their models, based on the modified triangular vegetation index, achieved R 2 = 0.79 and nRMSE = 24%. Another related study [17], in turn, used voxel-based methods from UAS-LiDAR data to invert LAI of maize and obtained nRMSE values of 10.8%, 12.4%, 42.8% for the upper, middle, and lower canopy layers, respectively. Maimaitijiang et al. [102] compared LiDAR and photogrammetry for LAI retrieval for sorghum fields and reported relatively better results from LiDAR metrics, with R 2 = 0.41 and RMSE = 0.32 for a tall, dense field and R 2 = 0.75 and RMSE = 0.27 for a lower, sparse field. These findings match our results concerning the apparent improvement of models' performance on the W cultivar group and the Flavor Sweet cultivar, and in general bodes well for UAS-based sensing of crop LAI. Major contributions of our study are the comparison of models across cultivars and especially sensing modalities. The use of image-based SfM versus more costly (financial and processing costs) LIDAR-based structural sensing is also a key finding.

V. CONCLUSION
The application of the UAS-based LiDAR and multispectral imagery for predicting LAI has distinct advantages in terms of temporal and spatial resolution and accuracy for field-tolandscape-scale crop management. Our study demonstrated the utility of LiDAR-derived models and MSI-derived 3-D models, via SfM, in predicting LAI of short broadacre crops like snap beans.
The relatively low canopy height of snap bean plants and foliage architecture made the internal canopy details hard to capture via either LiDAR or photogrammetry (SfM). Despite the challenges, this article strongly supported the potential of UAS-based LiDAR and multispectral imagery to estimate LAI of short broadacre crops. Furthermore, since snap beans only grow up to 0.3-0.6 m in height, the methods in our study should be extensible to other short broadacre crops, such as sugar beets, soybeans, and winter wheat (see Table I). Temporal transferability of the models was not fully explored due to limited, one-year sample data.
Results from this article are encouraging for the translation of an eventual operational solution to crop structural assessment, and potential extension to, for example, yield models, especially given that image-based SfM approaches performed nearly as well as LiDAR, an active sensing modality. Our objective here was to contrast LiDAR-and MSI-based LAI estimation, especially given that latter modality's cost benefits, i.e., SfM-based structural information, coupled to spectral content. However, we recommend that future studies evaluate the fusion of LiDAR and multispectral imaging, in order to assess the potential benefits of such a coupled modality approach. Larger sample size, multitemporal data, and more diversity in vegetation growth patterns all could contribute to the development of even more robust models. One more interesting topic is utilizing LiDAR indices to mitigate the saturation problem of VIs when LAI > 3 [103], which we would like to include in future studies.