AirMeasurer: open‐source software to quantify static and dynamic traits derived from multiseason aerial phenotyping to empower genetic mapping studies in rice

Summary Low‐altitude aerial imaging, an approach that can collect large‐scale plant imagery, has grown in popularity recently. Amongst many phenotyping approaches, unmanned aerial vehicles (UAVs) possess unique advantages as a consequence of their mobility, flexibility and affordability. Nevertheless, how to extract biologically relevant information effectively has remained challenging. Here, we present AirMeasurer, an open‐source and expandable platform that combines automated image analysis, machine learning and original algorithms to perform trait analysis using 2D/3D aerial imagery acquired by low‐cost UAVs in rice (Oryza sativa) trials. We applied the platform to study hundreds of rice landraces and recombinant inbred lines at two sites, from 2019 to 2021. A range of static and dynamic traits were quantified, including crop height, canopy coverage, vegetative indices and their growth rates. After verifying the reliability of AirMeasurer‐derived traits, we identified genetic variants associated with selected growth‐related traits using genome‐wide association study and quantitative trait loci mapping. We found that the AirMeasurer‐derived traits had led to reliable loci, some matched with published work, and others helped us to explore new candidate genes. Hence, we believe that our work demonstrates valuable advances in aerial phenotyping and automated 2D/3D trait analysis, providing high‐quality phenotypic information to empower genetic mapping for crop improvement.


Introduction
Rice (Oryza sativa) is one of the key staple foods, feeding > 50% of the global population (Muthayya et al., 2014). Breeding for rice improvements in yield potential, grain quality and resistance to stresses is vital for its climate-change adaptation and, thus, food security in many rice-consuming nations (Nakashima et al., 2007;Jagadish et al., 2012). This relies on selecting favourable phenotypes of agronomic traits from thousands of varieties in the field, which in turn heavily relies on specialists' visual assessment (Bevan et al., 2017;Roitsch et al., 2019). To help accelerate the selection procedure, many field-based phenotyping approaches have been introduced Yang et al., 2020).
Additionally, as agronomically important traits are controlled by the expression of multiple genes and modulated by the environment, phenotyping can assist researchers to understand underlying biological mechanisms that contribute to genetic gain (Hartung & Schiemann, 2014;Furbank et al., 2019). Through genome-wide association studies (GWAS), the genetic architecture of some agronomic traits in rice has been dissected Yang et al., 2014;Tang et al., 2019), laying the foundation of identifying functional diversity of alleles to discover valuable genes (Xing & Zhang, 2010). These contributions have led to advances in rice genetics and the development of new varieties with desired qualities, including high yield potential, resistance to stresses and increased resource-use efficiency (Barabaschi et al., 2015;Du et al., 2018;Li et al., 2018).
Certain traits such as plant height can be phenotyped at a specific time point; however, for growth-and yield-related traits that are genetically complex and influenced heavily by environmental factors, their phenotypes need to be examined dynamically (Naito et al., 2017;Mu et al., 2022). Nevertheless, to achieve this target, consistent data collection and trait analysis are required, which has posed significant challenges in developing reliable solutions for practical breeding programmes and field-based plant research (Shakoor et al., 2017;Pieruschka & Schurr, 2019). In essence, several problems need to be addressed, including: (1) scalability, trials are normally largescale and at multiple sites; (2) affordability, resources are usually limited and solutions need to be cost-effective; (3) accuracy and repeatability, analysis results should be consistent and reproducible in other trials; (4) processing cycle, the duration between breeding cycles or multiseason experiments is often brief, requiring data to be processed, analyzed and fed-back promptly to enable timely decisions (Großkinsky et al., 2015;Atkinson et al., 2018). Recently, several advances have been adopted by breeders and plant researchers, but many attempts remain at early stages (White et al., 2012;Juliana et al., 2019). New tools derived from some academic research have often worked at relatively small scale and with limited accessibility as a result of bespoke hardware, proprietary software and specialized packages, preventing them from being employed easily (Yang et al., 2020. Furthermore, to exploit genomic resources, traits of interest and genetic diversity need to be assessed across sites and seasons, demanding accessible data collection and analysis toolkits (Naito et al., 2017;Atkinson et al., 2018). Hence, methodological advances shall intend to address the above challenges, which is the emphasis of this study.
One of the most exciting advances recently was the rapid development of unmanned aerial vehicles (UAVs, also known as unmanned aircraft systems) and their applications in crop monitoring resulting from their mobility, throughput and affordability (Shi et al., 2016;Maimaitijiang et al., 2017;Jang et al., 2020). There are numerous examples in the literature reporting UAV-based phenotyping using sensors such as redgreen-blue (RGB) cameras, multi-and hyperspectral devices, Light Detection and Ranging (LiDAR), and thermal and infrared sensors (Kachamba et al., 2016;Gracia-Romero et al., 2017;Harkel et al., 2020;Hyyppä et al., 2020). Some work also reported quantitative trait loci (QTL) mapping of traits including plant height and vegetation fraction (Hassan et al., 2019;Wang et al., 2019;Ogawa et al., 2021). Nevertheless, most of these studies focused on estimating static traits collected at specific time points (Shakoor et al., 2017;Rodene et al., 2022), which often missed the dynamic nature of plant growth and development. Key agronomic traits (e.g. senescence and stem elongation) vary in time and space, which require new approaches to collect and analyse (Xu et al., 2018;. In fact, in a trial containing diverse genotypes, each line grows at a different pace, and thus dynamic analysis can provide meaningful comparisons between the genotypes (Hartung & Schiemann, 2014). Finally, changing behaviours of target traits, within or across seasons, can characterize the plant's complex responses to external stimuli, which are direct evidence to reveal spatial and temporal changes in the expression of genes and their regulators (Roitsch et al., 2019;Mu et al., 2022).
To extract meaningful information from UAV-collected imagery, many analytic solutions have been developed to measure traits related to yield, stress tolerance and growth patterns, using morphological, spectral and textural properties (Perez-Sanz et al., 2017;Jiang et al., 2021), most of which have focused on dryland crops. For example, EASY MPE (Tresch et al., 2019) combined excess green (ExG) and automatic thresholding to study soybean; AIRSURF (Bauer et al., 2019) employed deep learning to count and classify lettuces; GRID (Chen & Zhang, 2020;Tang et al., 2021) utilized pixel-wise Kmeans clustering to delineate irregular (e.g. zigzag) or regular (e.g. grid-based) trial layouts for wheat trials; R/UAS::PLOTSH-PCREATE  created polygon shapefiles using parameters (e.g. field direction and plot size) to study maize; FIELDIMAGER (Matias et al., 2020) incorporates manual inputs (e.g. row and column numbers) into the extraction of plot-based traits for potato.
Still, limited tools are available for nonexperts to examine multigenic traits and develop markers for paddy field crops (e.g. rice), which are complex as a consequence of changing water levels (e.g. resulting from rainfall and draining) and many voluntary plants (e.g. duckweed) compared with dryland crops (Ogawa et al., 2021). Moreover, few research groups have the resources to process large-scale aerial images, or to develop complex algorithms to address problems in automated trait analysis (Roitsch et al., 2019;Zhu et al., 2021). Hence, along with the development of open-source computer vision, machine learning and data science libraries (Howse, 2013;Virtanen et al., 2020), open solutions will be valuable to equip plant researchers with new toolkits to study complicated crops.
In order to address some of the challenges, we have developed AirMeasurer, an open-source platform that automates trait analysis for rice trials using 2D orthomosaics and 3D point clouds acquired by low-cost UAVs. First, we established tailored protocols for regular flight missions and data pre-processing. Secondly, varied 2D/3D analysis algorithms were integrated into the platform to quantify static traits such as seedling number, plant height, canopy coverage and vegetative indices, using morphological, spectral and textural signals. Thirdly, we developed an original algorithm to compute dynamic traits based on static traits, including growth rates of the target traits and their rapid growth phases, which were time-consuming or impossible to score previously. To ensure that our work could reach the broader research community, we created a graphical user interface (GUI) for nonexperts to use. Finally, to validate the platform and its utility in research, we applied the AIRMEASURER-derived traits collected from hundreds of rice landraces and recombinant inbred lines (RILs) in a multiseason case study (2019-2021) to genetic mapping studies (i.e. GWAS and QTL mapping) and identified reliable loci.

Plant materials and field experiments
In order to develop a UAV-based imaging protocol for multisite phenotyping, we established two experiments (2019-2021; Supporting Information Fig. S1a): (1) one focused on 254 landraces (Huang et al., 2012) in Shanghai (the 2019/2020 seasons), including 103 japonica, 40 intermedia and 111 indica types; (2) the other studied 191 homozygous RILs in Hainan (the 2020/2021 seasons), derived from the crossing parents Nipponbare (Oryza sativa ssp. Japonica) and 93-11 (Oryza sativa ssp. Indica), two popular varieties . In 2019, 177 RILs were used for manual assessment as a consequence of agronomic issues with some RILs during grain-filling. The sites were chosen owing to their geography and weather conditions. Crops at both sites were managed using standard husbandry and agronomic inputs according to local conditions. Landraces (Fig. S1b) and RILs (Fig. S1c) were sown in 2 × 1.1 m plots, 18 plants per plot. To maximize the efficient use of field space and facilitate initial selection (Payne, 2006), we did not introduce plot-level replicates; however, the same lines were repeatedly used in this multiseason case study. Details of the trial design, plant materials and geo-locations are provided in Notes S1.

Ground truthing
In-field ground truth measurements to validate AIRMEASURERderived traits were conducted by field workers. Maximum plant height was measured with a metre ruler in the late reproductive phase. After grain-filling, six plants in a plot were straightened and the distance from ground level to the top of rice spikes was measured. Heading date was scored manually, when there were five plants with panicles emerging 25 mm above the flag leaf sheath. To verify traits such as ExG and canopy coverage used for dynamic trait analysis, images of 29-30 randomly selected plots at six growth stages between early vegetation and early ripening (177 plots in total) were analyzed manually using the FIJI/IMAGEJ software (Schindelin et al., 2012), through which plot-based green-channel intensity values (0-255; measured from linear histogram) and canopy coverage (in pixels; using the autothresholding function) were obtained. To validate AIRMEASURERderived plant height at different growth stages, technicians manually measured calibrated 3D point clouds (with unwanted terrain features removed) to obtain plot-level canopy height at eight time points throughout the 2019 season (177 per point, 1416 in total).

Workflow of UAV-based phenotyping
When carrying out aerial phenotyping, we implemented a fourstep workflow (Fig. 1a): (1) experiment setupsincluding trial design (e.g. field layouts, target traits) and ground control points (GCPs; Figs 1b, S1d); (2) aerial imagingproviding guidelines to pilots to execute flight plans (Figs 1c, S1e); (3) data preprocessingproducing 2D field-level orthomosaic images (in TIFF) and 3D point cloud files (in LAS) from acquired aerial Fig. 1 A general workflow of unmanned aerial vehicle (UAV) based field phenotyping and phenotypic analysis established for collecting 2D/3D aerial images, processing 3D point clouds, and measuring plot-based morphological, spectral and textural traits. (a) A high-level workflow established to perform UAV-based field phenotyping and phenotypic analysis at multiple sites and over the course of multiple seasons. (b) Field experiments designed based on biological questions concerning plant varieties, target traits, treatments, trial layouts and in-field setups (e.g. ground control points, GCPs). (c) The selection of imaging protocols to collect aerial image series with 3D-or geo-referencing information. (d) Data pre-processing to produce 2D orthomosaics and 3D point clouds for the experimental field with plot-level plant resolution. (e) Automated trait analysis using a combination of 2D/3D image processing, spectral analysis, and machine learning techniques to perform plot segmentation and plot-based trait analysis using morphological, spectral, and textural signals (all the traits produced by AIRMEASURER are listed in Table 1

Aerial imaging using low-cost UAVs
At each site, in-field settings ( Fig. 1b) such as GCPs, height reference panels, spectral reflectance mats, or real-time kinematic positioning (RTK) were applied according to recommended practices published previously (Watanabe et al., 2017). To ensure that the imaging protocol could be adopted easily, we chose to use low-cost drones (e.g. Mavic 2 Pro; DJI, Shenzhen, China). Because smaller UAVs generated less downdraft, they could limit wind disruption of plant canopies during the low-altitude imaging. We designed two mission plans: (1) field-level imaging (25-35 m altitude), collecting RGB images speedily to limit colour distortion caused by natural illuminance (e.g. Fig. 1b left); (2) plot-level imaging, conducting flights with tailored flight parameters at low altitudes (10-15 m). Flights were normally carried out 10-12 times per season, among which eight flights were selected for time series measures (detailed mission plans, imaging protocols and guidelines are included in Notes S2).

3D point cloud processing and canopy height model
There can be unwanted noise in 3D point clouds generated by the Structure-from-Motion (SfM) algorithm (Singh & Frazier, 2018). To measure morphological features reliably ( Fig.  2a), we first denoised the SfM-generated 3D points (e.g. for a 0.1-ha field, low-density 3D reconstruction could produce > 30 million points). Second, we implemented the Statistical Outlier Removal (SOR) algorithm (Hodge & Austin, 2004) to remove outliers (red-coloured points, Fig. 2b). Third, a ground-level filter was developed based on the Cloth Simulation Filter (CSF) algorithm , classifying denoised 3D points into ground-level and aboveground groups. Because the CSF was designed for ultra-large land surveillance, we optimized the filter by reducing its grids and nodes (Fig. 2c). Finally, we removed unwanted terrain features (e.g. the field-level slope) using geo-or 3D-coordinates recorded from GCPs (saved in a shapefile, SHP). The procedure to correct geometric distortion is included in Notes S2 and S3, which shows the improved height measurements after removing field-level slopes. Next, we generated a digital surface model (DSM, i.e. ground-level points) and a digital elevation model (DEM, i.e. aboveground points) using the LidarTinGridding function (Lindsay, 2016) in WHITEBOXTOOLS (Fig. 2d). We defined region of interest (ROI) according to the SHP file (red markers, Fig. 2e). The DSM was subtracted from the DEM to retain aboveground plant information, resulting in a canopy height model (CHM) representing plant spatial signals with greyscale values (i.e. the brighter a pixel, the higher the point). Finally, we combined the CHM with spectral signals using the getPerspectiveTransform function (Mezirow, 1978) in OPENCV, which realigned the CHM (Fig. 2f, upper) with the field-level orthomosaic (Fig. 2f, lower). Spatial features were pseudocoloured (Fig. 2g,upper), ranging from 0 (dark blue, for bare plots) to 150+ cm (dark red, for tall plants). PYTHON-based software implementation for the above algorithms is given in Notes S4.

Plot segmentation
In order to acquire plot-based trait information routinely, plant plots should be identified consistently. Recent solutions such as EASY MPE, GRID, AIRSURF, R/UAS::PLOTSHPCREATE and FIELDIMAGER have been applied to segment plots or plant blocks for species such as soybean, wheat and maize, which were valuable advances for dryland crops. We trialled them in our paddy rice experiments and encountered segmentation issues as a consequence of unclear plot boundaries, changing water levels and overlapped rice plants during ripening (Notes S5).
Consequently, we developed an optimized plot segmentation algorithm (Notes S6): (1) applied an iterative self-organizing data thresholding (Irvin et al., 1997) to a field-level CHM and generated a global mask to represent plot edges (Fig. 2g, lower); (2) the Hough transform (Duda et al., 1972) was employed to seek horizontal and vertical lines in the mask, respectively (Fig. 2h, upper); (3) when some boundaries were undetectable, vertical and horizontal lines could be drawn manually to improve plot delineation via the GUI; (4) as most of the plots were not distanced evenly even with RTK-assisted seed drilling, we merged the adjacent lines; (5) after detecting plot boundaries, we assembled the remained lines to generate plot masks (Fig. 2h, lower), based on which all the plots were labelled according to the trial design for indexing and cross-referencing purposes (Fig. 2i, upper left); (6) to minimize edge effects and remove overlapped plants between neighbouring plots, a scaling function was designed to rescale the plot masks to measure different traits (e.g. scale = 0.25-0.3 for height measurements, depending on the degree of plant overlapping; Fig. 2i, upper right); and (7) finally, the refined masks were used to generate plot-level sampling regions (Fig. 2i, lower).

Automated trait analysis
Rice growth and development can be associated with stem elongation (i.e. changes in height) over time (Hosoi & Omasa, 2012). We utilized both spatial and spectral signals to analyze growthrelated traits. For different morphological traits, varied vertical levels of spatial signals were used. For example, we chose the top 10% height values (H 90th , i.e. top 10% of the 3D points; see reasoning in Notes S7) in the CHMs (scale = 0.25-0.3) to compute canopy plant height after grain-filling. For early establishment, as colour or textural signals were unreliable in identifying seedlings owing to weedy plants and changing water levels, we therefore first segment plant signals from CHMs (scale = 0.9; Fig. 3a); then, seedling masks were generated using H 95th as a result of short seedlings (Fig. 3b, left); after removing noisy objects (e.g. nongreen pixels) based on ExG values, we separated seedling objects from Comparably, we applied similar steps to measure canopy coverage before canopy closure: (1) using H 50th to represent plot canopy as canopy density was low during early vegetative phase (Fig. 3c, left); (2) applying the local adaptive thresholding (Singh et al., 2012) to generate a field-level mask (Fig. 3c, right); (3) overlapping the mask with 2D orthomosaic (Fig. 3d, left); (4) removing plot edges as some gaps between plots were unclear at canopy closure (scale = 0.7; Fig. 3d, middle); (5) using the Lab colour space (McLaren, 1976) to filter nongreen pixels (Fig. 3d, right); and (6) computing normalized canopy coverage index (CCI; 0 to 1, where 1 stands for 100% coverage; Notes S8).
According to a recent report (Svensgaard et al., 2021), RGB sensors can be applied to perform reliable spectral analysis without radiation calibration. Hence, we used RGB sensor to compute growth-related vegetation indices in the study. A series of vegetative indices and textural traits (e.g. canopy uniformity) were produced (Notes S9). All of the traits produced by AIRMEA-SURER are listed in Table 1.

Analysis of dynamic traits
Because dynamic or longitudinal phenotypes can be more informative in revealing plant-environment interactions (Campbell et al., 2019), we derived dynamic traits from static traits collected at different growth stages, rather than using values scored at arbitrary time points to represent growth patterns. Inspired by previous research (Anderson et al., 2019), we chose to measure dynamic phenotypes from the fitted curves even if some phenotyping points might be missing. The following section describes steps to compute dynamic phenotypes for an example trait, canopy height growth: (1) Eight height values were used between sowing and grainfilling for a given japonica landrace (red dots in Fig. 4a). The eight points were relatively evenly distanced between 10 and 115 d after sowing (DAS). Because the height of rice canopy tends to decrease during the later grain-filling period, we tested several fitting functions (e.g. stepwise regression) and chose the Gaussian function to fit plant height changes (green curve, Fig. 4a). Region of interest (ROI), denoted by four red markers recorded from ground control points (GCPs) with 3D-or geo-coordinates; then, DSM subtracted from DEM to generate a canopy height model (CHM), which uses greyscale values (0-255) to present plant height values. (f) A 2D perspective transformation applied to produce aligned red-green-blue (RGB) and CHM images using the ROI markers. (g) Pseudocolour applied to the aligned CHM according to a unified height scale bar (0-150+ cm; right); then, the iterative self-organizing data (ISODATA) thresholding algorithm employed to produce a field-level mask from the CHM. (h) The Hough transform algorithm used to detect horizontal and vertical lines separately, followed by the assembly of detected lines to produce initial plot masks. (i) All of the plots labelled based on the trial design; then, the scaling function applied to remove edge effects and overlapping plants among neighbouring plots, resulting in refined sampling regions for height (scale = 0.25-0.3) and colour-related measures in all the plots. (2022)

New Phytologist
(2) The Gaussian-fitted height curve f x ð Þ height then was used to generate a growth-difference curve f x ð Þ diff (black dash curve, Fig. 4a) through the KneedLocator function (Satopää et al., 2011), which measures value changes on f x ð Þ height , signifying the rate of plant height changes (i.e. increasing, decreasing or constant).
(3) Turning points (i.e. knee points, KPs; red crosses, Fig. 4a) were located on f x ð Þ diff , indicating height change phases. To locate theses KPs, we found the first We named the phase between KP1 and KP2 as rapid height growth phase (RGP height ; in days), denoting the period of rapid stem elongation.
(4) Within the RGP height , we found the first derivative f 0 x ð Þ height (green curve) to locate the day when canopy height was changing at the fastest growth rate (i.e. the FGR height day, in DAS; light-green cross; Fig. 4a) together with computing average growth rate (AGR height ; %), between 0 DAS and the FGR height day.
Then, we applied the above steps to analyze indica and intermediate landraces (e.g. GP014 and GP543; Fig. 4b). RGP height values for the genotypes were identified together with the FGR height days and AGR height . To assess phenotypic changes for other growth-related traits, we employed the algorithm to study variables such as ExG (i.e. RGP ExG , FGR ExG and AGR ExG ) and CCI (i.e. RGP CCI , FGR CCI and AGR CCI ). We also applied f x ð Þ height (green curves, Fig. 4c) to estimate Max height and other key growth stages (e.g. the beginning of ripening) using a normalized-curvature curve f x ð Þ cuv (dotted blue, Fig. 4c). The maximum curvature on f x ð Þ cuv was located to represent the Max height day (blue crosses, Fig. 4c), followed by the estimation of the beginning of ripening (purple crosses, Fig. 4c) using the minimum curvature. Moreover, AGR height , AGR ExG and AGR CCI (all in %) between 0 DAS and the Max height day also were quantified. To compute f x ð Þ cuv , we used the equation below: The rescaled plot masks applied to divide a field-level CHM acquired at early establishment (21 DAS) with a new height scale bar (0-20 cm; left), displaying height differences for short rice seedlings. Top 5% of 3D points (H 95th ) in a plot utilized to produce a plot-based seedling mask, followed by overlapping the mask with 2D orthomosaic (collected at 21 DAS; middle); finally, excess green index (ExG) computed to remove nonseedling objects, resulting in the quantification of seedling number per plot (right (f x i ð Þ height , Gaussian-fitted height curve; i is between 10 and 115 DAS).

PYTHON-based software implementation
A relative lack of open analytic solutions impedes researchers from exploiting newly introduced methods. Hence, we chose to develop the AIRMEASURER GUI using PYTHON programming language together with a modular design, so that each function or module in the platform could be accessed and modified independently. We used the TKINTER toolkit (Shipman, 2013) to develop a cross-platform GUI (in EXE). Open-source libraries such as SCIPY, OPENCV and SCIKIT-LEARN were employed to develop 2D/3D trait analysis algorithms and machine-learning based predictive modelling. The GUI software, executable JUPYTER notebooks, and user guides are provided for academic use (Availability and Requirements).

GWAS analysis and QTL mapping
The AirMeasurer-measured and manual-scored traits collected from landraces were used to perform GWAS analysis to find the associated-loci controlling phenotypes. The RIL population was used to verify the static and dynamic traits through QTL mapping. For GWAS analysis, an efficient mixed-model association eXpedited (EMMAx) was performed (Kang et al., 2010). Single-nucleotide polymorphisms (SNPs) with minor allele frequency (MAF) < 0.05 were excluded. Lines with missing phenotypes as a result of agronomic reasons were excluded. In total 254 landraces (2019/2020 seasons) were used to conduct GWAS. The matrix of pairwise genetic distances obtained by the simple SNP matching coefficients was employed to model the variance-covariance matrix of the random effect. Permutation tests were applied to help define the threshold of association signals (Churchill & Doerge, 1994). For each trait, we reshuffled the phenotypic data and performed association analysis using EMMAx with same parameters. To determine the significant threshold in GWAS, we accomplished 100 permutation analyses for each trait. Manhattan and quantile-quantile (QQ) plots were produced by using the PERL scripts (Mägi & Morris, 2010). For the 191 homozygous RILs (2020/2021 seasons), sequencing and genotyping were conducted using the published pipeline and SEG-MAP . Windows QTL CARTOGRAPHER v.2.5 (Wang, 2007) was employed for QTL analysis of composite interval mapping. LOD value was computed to indicate the possibility of QTLs based on likelihood ratio tests.

Collected 2D/3D aerial images
Using low-cost UAVs to monitor rice experiments between 2019 and 2021, many series of 2D/3D aerial images were generated. For example, eight flights conducted in the 0.1-ha trial in Shanghai generated 10 GB 2D/3D imagery from over 100 GB raw images in a season. For the 0.2-ha trial in Hainan, 13 GB 2D/ 3D imagery was created from eight flights (145 GB raw images).

New Phytologist
We uploaded a series of testing files to our GitHub repository (10 GB in total) for researchers to test and improve AIRMEASURER.

The GUI software
The initial GUI window of AIRMEASURER software consists of an input (red dash rectangle) section and a unified workspace (green dash rectangle; Fig. 5a). Users can select 2D/3D image series and a SHP file to begin the processing, including: (1) 'tab a' shows the central portion of the input orthomosaics within several seconds so that users can choose one image to proceed (Fig. 5b); (2) 'tab b' defines ROI and aligns the selected orthomosaic (Fig. 5c); (3) 'tab c' generates a CHM using the associated 3D point clouds and performs initial plot segmentation (Fig. 5d); (4) if some plot boundaries fail to delineate, users can use a mouse to draw horizontal and vertical lines to improve the plot delineation (yellow circles; Fig. 5e); (5) users can refine the plot masks using the scale function (0-1, where 1 stands for 100% of the original masks; Fig. 5d right); and (6) 'tab d' visualizes initial results and a button for batch processing with a progress bar and a checkbox for generating a performance matrix for genotypes (Fig. 5f). A different image can be reselected in 'tab a' to repeat the above procedure. The initial plot masks will be used to benchmark all the input images during batch processing. Finally, plot-based trait analysis and processed plot-level images can be downloaded.
Using an ordinary computer (Intel Core i7 CPU, 16 GB RAM with integrated graphics), 16 2D/3D images (10 GB) took 3 h to process. A detailed step-by-step user guide is provided in Notes S10 and Video S1.

Multiseason plant height analysis
We applied the AIRMEASURER system to process flights conducted in multiseason rice trials. For visual display, we selected three 2D orthomosaics to present overhead imagery and three 3D point clouds to exhibit field-level plant spatial features, from a 30°perspective, on 53, 73 and 103 DAS in the 2019 season, when landraces entered vegetative, reproductive and ripening phases ( Fig. 6a-c, left). Pseudocoloured height maps (Fig. 6a-c, right; with a unified scale bar) were created using AIRMEASURER-derived We applied the Gaussian-fitted curves and categorised the landraces into three groups according to their domestic types (i.e. indica, japonica and intermediary), with coloured shading areas denoting 15 th -85 th percentile confidence intervals (Fig. 6d). The fitted curves roughly followed a sigmoid pattern but with dissimilar developmental rates. For example, the intermediary group peaked on 91 DAS (Max height = 1.18 m), followed by the japonica group (93 DAS; Max height = 0.97 m) and the indica group (94 DAS; Max height = 0.82 m). For the 2020 height measures, three flights (on 52, 70 and 100 DAS) were selected to visualize plot-based height differences (Fig. 6e-g). The Gaussianfitted curves identified similar growth patterns (Fig. 6h); for example, the Max height days are between 95 and 110 DAS, and the 2020 lines were 5-10 cm shorter than the same lines studied in 2019. Complete height measurements for the two seasons (Datasets S1, S2), 2D orthomosaics and pseudocoloured height maps for the 191 RILs are provided (Fig. S3).

Performance matrix and dynamic trait analysis
In order to analyze dynamic traits effectively, we created a new function to help reorganize plant genotypes. For example, by extracting plot-level images from the 254 rice landraces and inserting them into a matrix according to their domestic groups, a 'performance matrix' was created using eight orthomosaics collected between 20 July and 8 October 2019 (Fig. 7a). In the matrix, each cell was an overhead image of a rice genotype, such that genotypes were columns and phenotyping timepoints were rows (see the entire 2019 performance matrix in Notes S11).
We used the matrix to examine different traits. The first row was utilized to quantify plot-based seedling number (Dataset S3; Fig. S2). To perform dynamic analysis of traits such as CCI, we used Gaussian-fitted curves to study the increase of CCI until 100% coverage was reached (Fig. 7b), based on which FGR CCI , AGR CCI and the Max CCI day were computed (Dataset S4). For spectral traits such as ExG and visible atmospherically resistant index (VARI), their associated dynamic traits were also quantified Graphic user interface (GUI) of AIRMEASURER developed for nonexpert users to readily use, which is capable of batch processing a series of 2D orthomosaics and 3D point clouds for 2D/3D trait analysis. (a) Initial GUI window of AIRMEASURER, consisting of input and analysis sections. A series of 2D orthomosaics, 3D point clouds, and 3D-or geo-coordinates (in SHP format) could be selected in the input section to initiate the initial analysis. (b) 'tab a' used to select an image with relatively clear gaps between plots from a list of input 2D orthomosaics. (c) 'tab b' used to define region of interest (ROI) of the field experiment using 3D-or geo-coordinates. (d and e) 'tab c' used to generate a field-level plant canopy height model (CHM) and plot masks. If the generated masks failed to delineate all the plot boundaries, an 'Optimize plot segmentation' button (coloured green) could be used to draw horizontal or vertical lines using a mouse (yellow circles); also, the 'Scale plot mask' input box could be used to scale down the plot masks (0-1, where 1 stands for 100% of the original mask), removing plot edges and overlapping plants. (f) 'tab d' used to visualise pre-processing results, a 'Batch processing' button to initiate automated trait analysis together with a progress bar and a checkbox for generating a performance matrix for all the rice genotypes. After the batch processing, trait analysis results (in comma-separated values, CSV), plot-based red-green-blue and CHM images (in JPG format) could be downloaded via the GUI. GUI-produced traits are listed in Table 1 S4). Noticeably, to estimate height-related traits (e.g. AGR height , RGP height and FGR height ; Dataset S9), both the matrix and associated CHMs were used.

Estimation of heading date
Based on dynamic traits, we explored the estimation of a complex trait, heading date. We predicted the trait using multiple dynamic traits (e.g. height, CCI and VARI) and machine-learning modelling (Notes S12), including: trait selection (Fig. S5a), dynamic trait analysis (Fig. S5b), feature engineering and selection (Fig. S5c), model training and selection (Fig. S5d), and model validation (Fig.  S5e). Among many models tested, support vector regression (SVR) obtained the best coefficient of determination (R 2 = 0.725; with P < 0.001 in linear regression; Notes S12).

Validation of AIRMEASURER-derived traits
AIRMEASURER-derived traits were validated by a range of ground truth data. The correlation between the 2019-AIRMEASURERderived Max height and the 2019-manual-scored Max height was visualized using 177 landraces (R 2 = 0.8848, P < 0.001; root mean square error, RMSE = 16.041; Fig. 8a), showing a very strong positive correlation. We repeated the validation with 254 landraces measured in the field in 2020. A similar R 2 was obtained (0.8926, P < 0.001, RMSE = 21.163; Fig. 8b).
Owing to different methods, AIRMEASURER-derived height values were consistently shorter than manual scoring, which was expected as plants were straightened in manual scoring. The result indicates that AIRMEASURER could soundly estimate maximum plant height among diverse rice genotypes (including landraces) with a high accuracy. We also verified the AIRMEASURER-measured height at eight time points (35-115 DAS) using 1416 plots (177 plots per time point) against plot-based canopy height that was measured manually from eight 3D point clouds (R 2 = 0.9651, P < 0.001, RMSE = 6.675; Fig. 8c), as well as against the Gaussian-fitted canopy height (n = 1416 plots; R 2 = 0.9649, P < 0.001, RMSE = 6.092; Fig. 8d). Significant positive correlations were obtained, indicating the reliability of the AIRMEASURER-estimated height trait throughout the season. For traits such as ExG and CCI that also were used for genetic mapping, 177 plots (at six growth stages, with 29-30 plots per stage) were used to compare the two traits obtained by manual and AIRMEASURER-based

QTL mapping using height-related traits
In order to evaluate the biological relevance of AirMeasurerderived traits in genetic mapping studies, we first used the overlapped 191 homozygous RILs in 2020 and 2021 season for genetic linkage analysis. The AIRMEASURER-derived Max height trait was used to map QTLs in the population, with the x-axis representing the genetic distance of 12 chromosomes and y-axis the LOD value. The threshold (red horizontal line) was set as 2.5 and known loci were indicated with red arrows. Two QTLs related to Max height were identified (Fig. 9a), among which one significant QTL (LOD = 15.7; chromosome 1) indicated a locus controlling rice plant height in the two seasons, consistent with the QTL mapping using the 2021 manual data and reported previously using the same population (Wang et al., 2011). In fact, there is a known gene sd1 that shortens rice stems (Sasaki et al., 2002), c. 100-210 kb from the locus. The second highest peak located using the Max height trait was on chromosome 7 (LOD = 7.3), c. 0 kb from Ghd7.1, a gene plays an important role in grain productivity and rice heading (Yan et al., 2013); however, the second peak identified with the manual scoring was c. 1 Mb from Ghd7.1 (LOD = 7.9; Fig. 9b).

QTL mapping using growth-related traits
Then, we mapped QTLs using other AIRMEASURER-derived growth traits such as CCI and ExG. The 2021 AGR CCI (0 DAS the FGR CCI day) was used to locate two QTLs (Fig. 9f), including Oshox4 (c. 0 kb; LOD = 7.9), overexpressing the gene leads to dwarfing and increased tillers, and thus the canopy size (Dai et al., 2008). One strong locus was identified on chromosome 9 (LOD = 9.6; Fig. 9g) using the 2021 AGR CCI trait (0 DASthe Max CCI day), indicating a locus linking to canopy expansion. Actually, TAC1 (Yu et al., 2007) is c. 10 kb away, which controls a spread-out or compact plant architecture. The  Table 2.

GWAS using height-related traits
Besides QTL mapping, we utilized AIRMEASURER-derived traits in GWAS analysis with the 254 rice landraces. We identified several significant SNPs associated with the two-season heightrelated traits and presented them in the Manhattan plot and QQ plot, with a grey dotted line indicating the threshold of the genome-wide significant P-value (Table S1) and a false detection rate (FDR) of 0.2. For example, using the 2019 Max height trait, the strongest signal on chromosome 1 (−log 10 (P) = 12, indicated with a blue arrow; Fig. 10a, left) was c. 208 kb from sd1. On chromosome 3, the strongest signal (−log 10 (P) = 6.42) identified was c. 10.7 kb from the OsHox32 gene, which is known for pleiotropic effects on plant architecture and leaf development . We repeated the analysis using the 2020 Max height trait and produced similar results on chromosome 1 (c. 208 kb from sd1; −log 10 (P) = 7.14; Fig. 10a, right). The findings were consistent with the GWAS analysis using the twoseason manual Max height scoring (c. 202-208 kb from sd1; −log 10 (P) = 6.91-6.63; Fig. 10b). Next, we chose the 2019 AGR height (0 DASthe Max height day) and identified four significant SNPs associated with the trait (Fig. 10c, left). Besides the strongest signal (−log 10 (P) = 10.62; c. 208 kb from sd1), the second strongest signal (−log10(P) = 6.98) was c. 11 kb from OsHox32 on chromosome 3, followed by a strong signal (−log10(P) = 6.6) c. 122 kb from NOG1 on chromosome 1, and the last one c. 373 kb from OsGSK2 on www.newphytologist.com chromosome 5. The NOG1 gene increases rice grain production (Huo et al., 2017), whereas OsGSK2 regulates the mesocotyl length (Sun et al., 2018), both genes relate to plant growth and development. We repeated the analysis using the 2020 data (Fig.   10c, right) and reproduced two SNPs that close to sd1 (−log10 (P) = 8.29; c. 188 kb) and OsGSK2 (−log10(P) = 8.07; c. 373 kb). To test other height-related dynamic traits, we used the twoseason AGR height (0 DASthe FGR height day) and obtained

GWAS using other growth-related traits
Finally, we chose the AIRMEASURER-derived CCI and ExG traits to perform GWAS and found three SNPs (Table S2). Using the 2020 AGR CCI trait (0 DASthe Max CCI day), two signals were identified (Fig. 10e): (1) one signal (−log 10 (P) = 7.4) was c. 329 kb from the Pit gene on chromosome 1, a disease resistance gene (Hayashi & Yoshida, 2009); (2) another (−log 10 (P) = 6.49) was c. 224 kb from the PFPβ gene on chromosome 6, which associates with carbon metabolism during grain-filling (Duan et al., 2016). Using the 2019 AGR ExG trait (0 DASthe FGR ExG day), the strongest signal (−log 10 (P) = 6.07) was c. 304 kb from the CCP1 gene on chromosome 1 (Fig. 10f), which functions  Fig. 9 Genetic linkage analysis of various AIRMEASURER-derived growth-related traits and manually scored maximum plant height, collected from 191 homozygous recombinant inbred lines (RILs) trialled in 2020 and 2021. For the significant single-nucleotide polymorphisms (SNPs) identified, known genes are indicated by red arrows. (a) Chromosomal location of significant quantitative trait locus (QTLs) identified using AIRMEASURER-derived Max height trait in 2020. The x-axis denotes the genetic distance of 12 chromosomes and y-axis the logarithm of the odds (LOD) value, with a significant threshold set at 2.5 (red horizontal line). The QTLs are close to the sd1 gene (chromosome 1) and the Ghd7.1 gene (chromosome 7). (b) Height QTLs identified using manually measured maximum plant height in the 2021 season; these also were located close to the sd1 and Ghd7.1 genes. (c) QTL for the AGR height trait, between 0 d after sowing (DAS) and the Max height day, in 2020. (d) QTL for the AGR height trait (0 DASthe FGR height day) in 2021. (e) Four loci associated with the RGP height trait collected in the 2020 season, including one located near SUI2 (chromosome 5), and another significant locus on chromosome 12 that is not associated with any known gene. (f) Two QTLs for the AGR CCI in 2021, determined over the period between 0 DAS and the FGR CCI day. The major QTL co-locates with Oshox4. (g) QTL for the average growth rate of CCI in 2021 determined over the period between 0 DAS and the Max CCI day. One strong locus on chromosome 9 (three peaks between 19.2 Mb and 21.6 Mb) co-locates with a known gene (TAC1) that controls canopy structure, and LGD1 that regulates vegetative growth in rice. (h) QTLs for the AGR ExG trait for the interval 0 DASthe Max ExG day. The major QTL co-locates with SLB1/SLB2 and D61.

Research Methods
New Phytologist palea development (Yan et al., 2015). Furthermore, GWAS was attempted with the 2019-heading-date trait estimated by both manual and AIRMEASURER approaches (Notes S13).

Discussion
In order to exploit available genomic resources to address climate change challenges, selected traits need to be assessed under field conditions across locations and years. Conventional phenotyping requires making many measurements of target traits, which is arduous and difficult to implement at busy periods of the season, resulting in newly developed methods (Pieruschka & Schurr, 2019;Jang et al., 2020). This study demonstrates that the use of low-cost UAVs can acquire larger and regular plant data from the field, based on which high-quality 2D/3D aerial imagery with field-and plot-level resolutions can be generated to enable automated analysis of static and dynamic traits that are biologically relevant. Furthermore, this approach is potentially valuable for assessing rates of genetic gain in larger trials, facilitating the calculation of heritability for agronomic traits and accurate genetic mapping for developing molecular markers. Nevertheless, metrics such as economic costs, scalability, analysis accuracy and throughput, or processing time need to be considered to evaluate the above research objectives, which is beyond the scope of this study but important for future studies.

Static and dynamic traits
The use of AIRMEASURER helped us analyze target traits coherently, which was achieved by removing unwanted field-level terrain features, transferring 3D to 2D signals for efficient processing, analyzing 3D points at different vertical levels and identifying plots consistently. These methodological advances were proved to be useful in examining static traits (e.g. height, CCI and ExG) at different growth stages for paddy rice (in particular the landraces), which was complex to study during the season. Inspired by previous work (Würschum et al., 2014;Anderson et al., 2019), we developed a bespoke approach to estimate dynamic traits derived from time series measures of static traits, enabling us to gain insights into dynamic features (e.g. growth rate and RGP) of target traits, without excessive phenotyping. Instead of using phenotypes measured at arbitrary time points, dynamic analysis helped us evaluate phenotypic variation reliably with hundreds of genotypes.
Additionally, through integrating dynamic traits into machine-learning modelling, we predicted a complex trait, heading date, which could lead to new estimates of QTL × Environment interactions. Recent studies (Lowry et al., 2019;Mu et al., 2022) have reported similar approaches that used multi-location traits in QTL mapping. Also, we demonstrated that AIRMEASURER-derived traits could be used for multiseason QTL discovery, which was confirmed by the results highlighting the locations of known QTLs. Although the main objective of this study was not to discover novel QTLs, nor to validate robustness of such QTLs across different germplasm sets and environments, it would only require simple adjustments to trial designs (e.g. more replicates) and greater repetition of trials across locations and years in order to produce reliable estimates of trait heritability and QTLs.

AIRMEASURER as a research tool
Genetic mapping of dynamic or longitudinal traits can be a powerful tool for developing novel molecular markers that cannot easily be revealed using static measurement, partly because of temporal regulation of gene expression (Harder et al., 2019). We explored the use of AIRMEASURER-derived traits to identify associated loci. For example, QTLs were mapped for traits such as Max height , AGR height during RGP height or at the FGR height day. If the QTLs are shown to be robust across years, locations and different germplasm sets, these then could be used to develop growth-related molecular markers. Some of the identified loci were co-located with known genes, as well as with other genes within the interval had unknown functions, which could lead to new candidate genes. QTLs located using traits derived from CCI and ExG (e.g. Max cov , FGR cov and Max ExG ) also indicate potentially useful loci. Likewise, we used AIRMEASURER-derived traits in GWAS analysis. Comparable loci were identified from rice landraces. Height-related trait (e.g. AGR height ) led to the consistent identification of signals such as the nearby genes (e.g. sd1, OsHox32, NOG1 and OsGSK2) relevant to plant height, architecture and growth regulation, indicating the value of dynamic traits in studying genetically diverse landrace populations. Moreover, using dynamic traits such as AGR height , Max cov and FGR ExG , we located some previously unknown strong signals, which may be valuable for identifying small effects of individual allelic differences (e.g. loci on Chromosome 2; Fig. 10c,d) that jointly contribute to the regulation of trait expressions. Finally, AIRMEASURERestimated heading date traits could bring a new perspective to GWAS analysis. Loci identified using a small number of indica landraces (n = 97) were just c. 17.14 kb from OsSOC1 gene on chromosome 3 (−log 10 (P) = 5.79) and c. 30.24 kb from the Hd3a gene on chromosome 6 (−log 10 (P) = 4.65), both of which are heading-date related (Notes S13).

Limitations of the platform
We have encountered problems that are not uncommon when applying drones in aerial phenotyping: (1) weather conditionssmall UAVs cannot be operated in unstable weather such as high or gusty wind (> 15 ms −1 ), rainfall, or heavy fog (Tmušić et al., 2020); (2) geo-referencing -GPS modules installed on lowcost drones had metre-level deviation and thus geo-referencing errors needed to be rectified; (3) nature illuminanceimage colour and contrast could vary noticeably with changing light conditions, and we mitigated this issue by conducting a field-level imaging; and (4) aviation regulationsthe change in aviation regulations casts uncertainty on aerial phenotyping, requiring regular communications with local civil air traffic control authorities; without official authorization, the payload capacity of a drone was restricted, indicating the advantage and practicality in using small drones for routine phenotyping. In the study, we used an RGB camera for growth-related spectral analysis as a low-cost alternative to more costly hyper-or multispectral sensors. It is worth noting that visible spectra are limited in early diseases detection and sensing abiotic stress responses as accurate spectral information is key to assess plant responses to certain external stimuli (Tmušić et al., 2020).
For data pre-processing, we used the proprietary PIX4DMAPPER software to generate 3D point clouds and 2D orthomosaics. We have tested several open-source software types (e.g. VISUALSFM, MESH-ROOM) for the same task and encountered technical problems such as prolonged computational time, incorrect geo-referencing, and mismatched 2D/3D patches. Another problem during the processing was to denoise large-scale 3D point clouds. We used the SOR for the task, which required 15-20 min to denoise 60+ million 3D points. Hence, algorithms such as local-outlier and cluster-based outlier detection (Kriegel et al., 2009), and deep-learning (DL) methods (Casajus et al., 2019) should be considered to speed up this task. Although the AIRMEASURER's plot segmentation algorithm could reliably be applied to field experiments with regular gridded plot layout designs (Notes S10), it cannot be extended to analyze irregular plot layouts (e.g. zigzag arrangements). DL approaches such as multilayer perceptron that can incorporate multidimensional ground/plant signals might be more useful for this mission.

Future applications
Further developments could include the analysis of high-density 3D point clouds. Rice flowering starts 1 d after the heading, during which anthers (1-1.5 mm in diameter) can be observed on different panicles. By flying smaller UAVs (e.g. DJI Mini) at a 4-m altitude, we could achieve a ground-sampling-distance (GSD) of 1-1.5 mm per pixel. Thus, it is feasible to measure anther extrusion using highdensity 3D points acquired by smaller drones, which could find applications in hybrid breeding programmes where the selection of male parents with certain flowering characteristics is crucial.
Low-cost UAVs and dynamic trait analysis also could be applied to examine traits such as grain-filling, which are challenging to quantify using conventional approaches. By conducting daily flights during ripening, fitted curves could enable the estimation and eventually the prediction of the initiation and duration of this key trait. However, it is expected that the Gaussian function might not be suitable for such growth patterns, and thus other fitting methods shall be explored.
Although we did not thoroughly test AIRMEASURER to analyze other crops, we have successfully applied the platform to examine wheat trials with limited parametric changes (Notes S10, S14), suggesting potential applications of AIRMEASURER for other plant species. As the modular-designed AIRMEASURER was developed in PYTHON, which is widely supported, we trust that this platform could be shared, extended and upgraded by the community relatively easily, providing open and readily accessible solutions for the broader research community.

Supporting Information
Additional Supporting Information may be found online in the Supporting Information section at the end of the article.

Dataset S4
The 2019 season dynamic analysis of canopy coverage index.

Dataset S5
The 2019 season dynamic analysis of the excess green index.

Dataset S9
The 2019 season dynamic analysis of canopy height.

Fig. S5
Combining phenotypic traits and supervised machine learning to predict a complex trait, heading date, with high confidence.
Notes S1 Trial design and plant materials.
Notes S2 UAV imaging protocol and in-field setups.
Notes S3 Different plant height measures before and after removing terrain features.
Notes S4 3D point clouds processing and canopy height model.

Notes S5 Previous published segmentation solutions trialed in rice field experiments.
Notes S6 Source code of the plot segmentation algorithm.
Notes S7 The reasoning behind choosing H 90th for height measurement.
Notes S8 Source code for computing canopy coverage and ExG indices.
Notes S9 Vegetative indices and texture-based traits.
Notes S10 A step-by-step user guide of the AIRMEASURER GUI.
Notes S11 Entire performance matrix for all plots monitored.
Notes S12 Estimation of a complex traitheading date.
Notes S13 GWAS using heading dates estimated by the SVR model.
Notes S14 Applying AIRMEASURER to examine wheat plots under different nitrogen treatments.
Video S1 GUI of AIRMEASURER in operation.

New
Phytologist