UAV‐based imaging platform for monitoring maize growth throughout development

Abstract Plant height (PH) data collected at high temporal resolutions can give insight into how genotype and environmental variation influence plant growth. However, in order to increase the temporal resolution of PH data collection, more robust, rapid, and low‐cost methods are needed to evaluate field plots than those currently available. Due to their low cost and high functionality, unmanned aerial vehicles (UAVs) provide an efficient means for collecting height at various stages throughout development. We have developed a procedure for utilizing structure from motion algorithms to collect PH from RGB drone imagery and have used this platform to characterize a yield trial consisting of 24 maize hybrids planted in replicate under two dates and three planting densities. PH data was collected using both weekly UAV flights and manual measurements. The comparisons of UAV‐based and manually acquired PH measurements revealed sources of error in measuring PH and were used to develop a robust pipeline for generating UAV‐based PH estimates. This pipeline was utilized to document differences in the rate of growth between genotypes and planting dates. Our results also demonstrate that growth rates generated by PH measurements collected at multiple timepoints early in development can be useful in improving predictions of PH at the end of the season. This method provides a low cost, high throughput method for evaluating plant growth in response to environmental stimuli on a plot basis that can be implemented at the scale of a breeding program.


| INTRODUC TI ON
Plant height (PH) serves as a major growth indicator and can be used for assessing crop productivity and making crop management decisions. PH has been shown to be linked to nitrogen (N) nutrition during vegetative development in maize (Yon, Jaja, Mcclure, & Hayes, 2011;Freeman et al., 2007), making it useful for assessing spatial variability in crop response to N. PH can, therefore, help guide precision agriculture practices such as variable-rate N applications within the field. PH has also been linked to plant biomass in maize and barley crops (Bendig et al., 2014;Freeman et al., 2007;Li et al., 2016). Studies have found that PH measurements at early to mid-developmental stages in maize are correlated to grain yield Katsvairo, Cox, & Van Es, 2003), and can be useful for forecasting crop yields  or improving yield predictions (Sharma & Franzen, 2014). These findings suggest that tracking PH at these earlier stages and throughout development can be useful for identifying superior cultivars in plant breeding programs and for developing management practices that account for spatial heterogeneity in production fields. However, there is still substantial variation across these studies in the level of utility of early-season PH in predicting end of season traits such as yield making this a topic that needs to be further evaluated. To maximize the utility of PH data, and to address some of the uncertainties regarding the utility of PH, data for large-scale experiments at various stages throughout development needs to be gathered.
With current methods for measuring height, this can be challenging. Similar to many phenotypic traits, current practices of gathering height in field settings involve physically measuring PH with a large ruler, which is time consuming and difficult to implement on a large scale. Ruler measurements are also subject to user bias and error, decreasing the accuracy and utility of these measurements.
Remote sensing technology has been proposed as an efficient means for collecting rapid and objective measurements for PH (Araus & Cairns, 2014;Madec et al., 2017;Yang et al., 2017).
Unmanned aerial vehicles (UAVs) have proven to be a promising platform for collecting remote sensing data as they are inexpensive yet capable of achieving very high spatial and temporal resolutions.
Many modern UAVs also have the functionality of flying automatically along a path specified by the user through mission planning software, increasing reproducibility through time. Commercial UAV platforms come equipped with RGB cameras that can be utilized to capture images along the field to create 3D field reconstructions using structure from motion (SfM) and multiview stereo (MVS) algorithms which rely on estimating the 3D structure of a scene utilizing a set of 2D images (James & Robson, 2012). The speed and ease with which UAVs can be flown for data acquisition provide advantages in scale and temporal resolution over other remote sensing methods utilized for creating 3D models, such as light detection and ranging (LiDAR), which involve high costs in data acquisition and processing (Díaz-Varela, Rosa, León, & Zarco-Tejada, 2015). However, there is substantial improvement that is needed in developing tools and approaches for simple implementation of UAV technology to generating phenotypic data in field settings.
There are examples of success in using RGB sensors to extract a number of traits that, if taken through time, can give us insight into plant growth and the relative effects of genetic and environmental variables on these traits. Canopy coverage estimates have been gathered for soybeans using digital cameras and these have been found to highly correlate to canopy light interception measurements (Purcell, 2000) and grain yield (Xavier, Hall, Hearst, Cherkauer, & Rainey, 2017). Previous studies have estimated PH from UAV imagery in sorghum (Chang, Jung, Maeda, & Landivar, 2017;Shi et al., 2016;Watanabe et al., 2017), wheat (Holman et al., 2016;Madec et al., 2017;Michalski et al., 2018), cotton (Feng, Sudduth, Vories, Zhang, & Zhou, 2018), barley (Bendig et al., 2014), and maize (Anderson II et al., 2019;Anthony, Elbaum, Lorenz, & Detweiler, 2014;Geipel, Link, and Claupein, 2014;Grenzdörffer, 2014;Su et al., 2019;Varela et al., 2017). Although these studies have reported high correlation in PH measurements for various dates pooled together to manual measurements, they lack estimates of how daily correlations of imagery-derived PH (PH UAV ) and ruler-derived PH (PH R ) vary throughout different growth stages and how this compares to the inherent error in PH R measurements.
Taking into account past efforts in height estimation from UAV platforms, the present work aims to implement SfM algorithms for the estimation of PH UAV for a maize field throughout development and compare these measurements to daily PH R measurements collected using the traditional manual ruler-based method. We tested the error in PH R measurements and assessed the correlation of PH UAV to PH R measurements. We also evaluated applications in which this method could be used such as in estimating terminal height from early-season PH UAV data.

| Experimental design
Two replicates of 12 maize hybrids were planted in 4-row plots at two dates (14 May 2018 and 21 May 2018) and three densities (60 k, 90 k, and 120 k plants per hectare) following a randomized complete block design blocked by planting date and replicate within planting date in St Paul, MN in the summer of 2018 ( Figure S1). Each row was 15 ft long center-to-center with 12 ft of plots and 3 ft of alleys and 30 in spacing between rows. There were a total of 144 four-row plots planted for this experiment. The 12 hybrid genotypes utilized were generated by crossing ex-PVP lines selected due to their past use in production settings and the availability of seed. This experiment was grown on two acres. Nine 1 × 1 m ground targets were distributed around the border and internal alleys of the field for use as ground control points (GCPs) based on previously developed optimization algorithms (Gomez-Candon, De Castro, & Granados, 2014).

| Drone image data collection
The field was imaged approximately weekly from planting until plants reached terminal height using a DJI Phantom 4 Advanced drone flown 30 m above ground to achieve a ground sampling distance (GSD), measured as the ground distance represented between the centers of two neighboring pixels, of approximately 0.82 cm.
SfM algorithms for 3D reconstruction of field points from 2D images rely on finding correspondences between images. Wind can cause issues by causing the object to move between two consecutive recordings. This can be reduced by using a high measuring repetition rate (Paulus, 2019). Images were collected in a grid pattern using 85% frontlap and 85% sidelap to maximize reconstruction efficiency ( Figure 1b). A total of 615 images were gathered per mission and a total of twelve missions were conducted throughout the season.

| Manual measurement data collection
Manual measurements for PH (PH R (2)) were collected on the same day as eight of the drone flights by measuring two plants per row from each of the middle two rows of each four-row plot using a ruler.
PH R was measured as the distance between the ground and the top-

| Data processing workflow
Agisoft software (Agisoft Photoscan Professional Edition v1.4.4, 2018) was used to process the collected imagery and generate crop surface models (CSM) and RGB orthomosaics for each date following the developed workflow ( Figure 1a). Processing steps implemented in Photoscan include feature matching, solving for camera intrinsic and extrinsic orientation parameters, building a dense point cloud, building a field orthomosaic, and building a digital elevation model as specified by Photoscan manual recommendations (Agisoft LLC, 2018; Figure 1c). Parameters set for image alignment included: high accuracy for obtaining camera position estimates, reference pair selection so that overlapping pairs of photos are selected based on measured camera locations, a key point limit of 40,000 points, and a tie point limit of 4,000 points. GCPs were then used to optimize camera alignment. For generating the dense point cloud, the quality value was set to High with Moderate depth filtering to optimize processing time and model quality.
QGIS software (QGIS v2.18.9, 2017)  as ground height estimates for a subset of plots and the resulting plot height estimates compared to each other to select the optimal ground height extraction procedure. All three methods utilizing the different percentiles represented plot ground height values in the alleys; however, the 3rd percentile was selected as this value reduced the amount of error in ground height extraction compared to the 1st percentile due to outlier points present in the 3D models and the subsequent DEMs that were generated. This percentile also produced plot height estimates that had a higher correlation to ground truth measurements compared to utilizing the 5th percentile. The 97th percentile was selected following the same logic and to keep the procedure consistent. To extract the PH UAV value for each bin, Only the middle 12 bins for the middle two rows in each plot were assessed when calculating an average PH value to prevent the average height for a given plot to be biased by the height of plants near alleys or different density treatments. A piecewise polynomial spline of order 20 was fit using the MATLAB splinefit function to better represent the continuous plant growth that occurs throughout the growing season (Lundgren, 2019). To assess how growth changes through time, we then extracted the first derivative, or slope, of the curve to get change in height through time.

| Statistical analyses
A total of 1,152 mean plot PH UAV values were extracted across eight timepoints. Outliers were identified for each timepoint as those values that were less than or larger than the mean PH UAV value for that date plus or minus the standard deviation multiplied by a factor of 21 2. A total of 17 outlier points across all dates were identified and removed from the data set for subsequent analysis. It should be noted that while there was a strong wind event that caused substantial lodging in the experiment (see below), these outlier points were random and not enriched for data points collected at or shortly after the lodging event. To measure the strength of the linear relationship between different PH measurements, Pearson's correlation was utilized.
Variation in mean plot PH values was partitioned into genotype, planting date, density, genotype-by-density interaction, genotypeby-planting date interaction, and replicate all treated as fixed effects with the linear model where y ijkl is the phenotype value of the ith genotype in the jth planting date in the kth planting density of the lth replication; u is the phenotypic mean across planting dates and planting densities; g i is the ith genotype effect; e j is the jth planting date effect, d k is the kth planting density effect; r (e) l(j) is the lth replication effect nested within the jth planting date; (ge) i * j is the interaction effect of the ith genotype by the jth planting date; (gd) i * k is the interaction effect of the ith genotype by the kth planting density; and e ijkl is the residual effect. To test the significance of the various effect variables of the linear model on PH, we used the ANOVA function of the R stats package (R Core Team, 2012).

| Model development for predicting terminal height
A piecewise polynomial function of 20 degrees was fit to PH UAV values for each plot through time to create a smooth curve with optimal fit. The first derivative of the fitted curve was then extracted to get F I G U R E 2 Method for extracting mean PH UAV values for a given plot. The plot is segmented and broken down into 20 bins. The 3rd and 97th percentiles for each bin are extracted and used to calculate the average plot height by subtracting the minimum 3rd percentile of all bins from the 97th percentile of each bin change in height through time (Figure 6b). The slope curve was then broken up into 100 equidistant timepoints for each plot. The contri- provided the largest contribution, in this case more than 40%, of the observed variation for terminal hand height. The final model contained the terminal height as the response variable and the slope at nine timepoints as the predictor variables ( Figure 6b). This model was derived from 1/4 of the data points and tested to predict the terminal height of the other 3/4 of the plots.

| RE SULTS AND D ISCUSS I ON
Our goal was to implement a robust low-cost platform for efficient estimation of PH for maize plots. High-throughput collection of PH measurements throughout the growing season can provide a better understanding of how different environments influence maize growth and can help to dissect the underlying cause of end of season genotype × environment (G × E) interactions that are commonly observed. In order to develop a robust PH collection pipeline and to begin to consider the factors that influence plant growth rates, we planted a set of hybrids in a design that incorporated multiple planting densities and planting dates such that plants would be exposed to similar environments but at different growth stages ( Figure S1). To decrease the impact of neighboring plots, each plot was grown as a four-row plot and only plants in the middle two rows were measured.
Manual measurements were taken using a ruler on the dates that aerial images were obtained. Analysis of variance showed significant variation in final PH within this experimental design due to genotype, planting date, and density (Table 1).
Aerial images of the field were collected throughout the growing season until plants reached terminal height (approximately weekly, weather permitting). The images were used to create high resolution 3D models, field orthomosaics, and DEMs using existing SfM algorithms ( Figure 1). There are multiple methods that have been utilized for extracting PH from UAV images. A common approach used is

| Error within ruler height measurements
A critical question for any application using PH UAV is the quality of the height estimates that are obtained through this high throughput measurement. However, there are multiple complicating factors in truly assessing the quality of PH UAV estimates. An important factor is that the manual measurements of PH (PH R (2)) that are used as the "ground truth" for determining the accuracy of PH UAV estimates can include some level of error. In order to begin to understand the quality of the plot mean PH R (2) revealed variation in the correlations and root mean square errors (RMSEs) (  Table S1).
Another potential difference between PH R and PH UAV estimates of plot height is that manual measurements are often collected for a subset of plants in a row while the UAV measurements are often generated as a per plot mean height that samples a larger number of plants.
To investigate the accuracy in plot mean PH R (2) Table S1). These correlations are comparable to those obtained when correlating plot mean PH UAV measurements to plot mean PH R (2) measurements ( Figure 4a).
This indicates that PH R measurements can have substantial variation that is impacted by the number of plants measured and the selection of plants to be measured.

| Correlation of ruler height to UAV-derived height measurements
There are limitations in accuracy of ruler-estimated heights for specific plants as well as for estimating the height of a plot. However, these measurements can still be used to assess the quality of PH UAV measurements, if we accept that the quality of the correlations between PH UAV and PH R likely cannot be improved beyond the inherent errors in PH R measurements. Compared to our manual height estimates, analysis of variance showed significant variation in terminal PH UAV within this experimental design due to genotype and planting date; however, density had a smaller effect (Table S2). Comparison of PH R (2) and PH UAV showed the overall correlation across all data collection dates was strong (adj. r-square value of .96 for our dataset; If we compare the error in PH UAV measurements to that in PH R (2) replicated measurements for the subset of plots where two individual plants were hand-measured twice, we see that the correlations of PH R (2) and PH UAV measurements for the same plots increased throughout development, but were lower than the correlation between PH R (2) replicated measurements (Table S1). By terminal height, the correlation between PH R (2) and PH UAV measurements was high (adj. r-square value of .72) but had a substantial amount of error ( Figure 3c). Unlike the high correlation and low RMSE obtained from plot mean PH R (2) measurement replicates calculated by measuring the same two plants per row, plot mean PH R (2) measurements obtained by sampling different sets of plants within a row showed a much lower correlation and larger error across most timepoints ( Figure 4c; Table S1). These correlations between plot mean PH R (2) calculated by measuring varying plants in a row, however, are comparable to those obtained when correlating plot mean PH UAV measurements to plot mean PH R (2) measurements (Figure 4a-c) correlation to PH R (2) and high error rate throughout most individual timepoints with an RMSE ranging from 15.32 to 48.73 cm (Figure 4d).
These lower correlations and higher error estimates compared to the PH UAV measurements obtained through our approach show the value in binning and extracting an average to represent height for plots, where height variation is present rather than extracting a single value.

| UAV-estimated heights can detect biological variation
Our goal was to develop a low cost, high-throughput platform for monitoring PH through time that would allow us to look at how genotypes develop and respond to different environmental conditions. A storm with high wind speeds occurred on July 1st of the 2018 growing season in St. Paul, providing a prime opportunity to evaluate the utility of this method to track plant responsiveness to environmental conditions. An example of this responsiveness can be seen in the height profiles for LH82 x DK3IIH6 plots, where plants at high density in the earlier planting date showed a reduction in height 45 days after sowing ( Figure 5a). Normally, a reduction in height would be unexpected and may reveal issues with the height estimation platform, but this height reduction was due to the lodging that occurred as a result of the storm. We were not able to track the lodging event in the PH R measurements (Figure 5b), as these measurements were not taken until several days after the wind storm and most plots recovered from the lodging prior to this measurement.
UAV imagery provides an opportunity to collect data rapidly in response to events such as this storm event and track plants at regular intervals that is much more difficult to do with hand measurements.
We expected to note differences in height within genotypes and planting dates based on planting density. Plants grown at higher density tend to grow faster and taller when water and nutrient resources are not limited as they are competing with more plants for sunlight (Tetio-Kagho & Gardner, 1988). Planting density did contribute to differences observed for lodging across genotypes in the early planting date. Plots planted at higher densities experienced a more severe degree of lodging when exposed to strong winds compared to plots planted at a lower density (Figure 5c). However, our data suggested that the realized height remained pretty stable across densities and that in many cases, the differences in height due to genotypic variation were reproducible across planting dates and densities ( Figure S2a-c).
The PH UAV estimates were also able to document differences in growth rates throughout the growing season between genotypes.
For example, LH82 × LH145 planted at high density in the late planting had a faster rate of growth compared to LH82 × DK3IIH6 planted in the same conditions as soon as 20 days after sowing ( Figure S2c).
LH82 × LH145 was able to maintain an advantage in terms of height throughout development and achieve a higher realized height at the end of the season. The ability to track these differences in growth patterns among genotypes that may or may not impact the end of season heights that are typically collected provides additional information to breeders in the process of selecting superior parents and hybrid combinations. Overall, visualization of the PH profiles for each genotype × condition suggested a fairly robust ability to monitor biological variation for height and responsiveness to environmental conditions under different growth conditions.  Figure S3) and mean plot PH R (2) (adjusted r-squared value of .60; Figure 6c). When looking at the different treatments separately, we can see that the ability to predict height across planting dates increases for the late planting treatment and the magnitude of this difference increases with higher planting densities. This is likely due to the challenges in estimating mean plot height for plots that suffered lodging and the higher density treatments in the early planting date experienced the most lodging early in the season.

| CON CLUS IONS
High-throughput collection of PH measurements throughout the growing season can allow a better understanding of how different environments influence maize growth and can reveal important factors affecting crop productivity. We have developed a novel approach for extracting plot PH UAV from individual plots that involves estimating the ground height individually for each plot and using that to calculate mean plot height. Overall, PH R (2) and PH UAV measurements showed a high correlation (adj. r-square value of .96); however, correlations within PH estimates collected on a single date are much lower as they are impacted by the low variability in height expressed on a single date. The reproducibility and quality of PH R (2) estimates were found to vary throughout development with lower correlations and higher error rates observed in early vegetative stages, and improved correlations and lower error rates observed in later developmental stages.
When assessing the error in generalizing mean PH R for a plot by measuring only a subset of plants in a row, PH R measurements were found to be greatly variable and impacted by the number and selection of plants sampled. Plot mean PH R (2) measurements obtained from different sets of plants within a row show a correlation and error across most timepoints comparable to those obtained when correlating plot mean PH UAV measurements to plot mean PH R (2) measurements.
Implementing this pipeline to visualize height across different genotypes, planting densities, and planting date treatments showed a strong ability to monitor biological variation for height across development in response to the environment. Differences in rate of growth between genotypes and planting dates as well as in lodging responses between planting densities were documented. Moreover, we saw that growth rates generated by PH measurements collected at multiple timepoints early in development can be useful in improving predictions of PH at the end of the season. Together, these findings show the utility of using UAVs for collecting PH data at high temporal resolution to track differences in growth and to predict an end-season trait early in development.

ACK N OWLED G M ENTS
We would like to thank Amanda Gilbert and Pete Hermanson for technical support for this experiment. We would also like to thank

CO N FLI C T O F I NTE R E S T
The authors do not have any conflict of interest to declare.

AUTH O R CO NTR I B UTI O N S
SBT, CNH, and NMS conceived the experiments; SBT conducted the experiments and analyses; SBT, CNH, and NMS wrote the manuscript.

DATA AVA I L A B I L I T Y S TAT E M E N T
The scripts and processes used to perform the image analyses and trait extraction are available at https://github.com/SBTir ado/ UAV_PH.git.