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

Abstract

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.

Keywords: 2D/3D trait analysis; aerial phenotyping; genetic mapping; predictive modelling; rice; static and dynamic traits.

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).

Fig. 2

Algorithmic steps for processing unmanned aerial vehicle (UAV) collected 3D point clouds to generate aligned canopy height model (CHM) within region of interest (ROI) together with plot segmentation for plot‐based trait analysis. (a) A 3D point cloud file produced from pre‐processing (in LAS format). (b) Outliers (red) removed in the point clouds using the Statistical Outlier Removal (SOR) algorithm. (c, d) The Cloth Simulation Filter (CSF) algorithm applied to differentiate ground‐level and aboveground 3D points, resulting in a digital elevation model (DEM) and a digital surface model (DSM). (e) 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.

Fig. 3

Algorithmic steps for quantifying plot‐based morphological traits such as rice seedling number and canopy‐related traits using both spatial and spectral signals. (a) Plot masks rescaled (scale = 0.9) to segment a canopy height model (CHM) image collected at 93 d after sowing (DAS); the segmented CHM (right) pseudocoloured according to the unified height scale bar (0–150+ cm). (b) 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). (c) A field‐level CHM (69 DAS) used to compute canopy coverage index (CCI; left); top 10% of 3D points (H _90th) in a plot used to create plot‐based canopy masks (right). (d) After overlapping the canopy masks with 2D orthomosaic (69 DAS; left), edges of the canopy removed using the scaling function (scale = 0.7; middle), resulting in refined plot‐based canopy regions for computing canopy coverage and canopy ExG indices.

Fig. 4

Algorithmic steps for quantifying dynamic phenotypes of an example trait, canopy height growth rate; three types of rice landraces are shown to illustrate the procedure and the capability of estimating growth‐related traits. (a) Eight canopy height values (red dots) recorded between sowing and grain‐filling for a given japonica landrace, which were relatively evenly distanced during key growth stages, between 10 and 115 d after sowing (DAS). The Gaussian function applied to produce a growth curve of canopy height ($f{\left(x\right)}_{\text{height}}$; green colour), based on which a growth‐difference curve ($f{\left(x\right)}_{\text{diff}}$; black dash curve) was created. $f{\left(x\right)}_{\text{diff}}$ measures value changes on $f{\left(x\right)}_{\text{height}}$, indicating the rate of canopy height change during the season. Turning points (i.e. knee points, KPs; red crosses) on $f{\left(x\right)}_{\text{diff}}$ located to represent the rapid growth phase of canopy height (${\mathrm{RGP}}_{\text{height}}$; in days; red shading area), indicating the most rapid period of stem elongation. Within the ${\mathrm{RGP}}_{\text{height}}$, the fastest growth rate (${\mathit{\text{FGR}}}_{\mathit{\text{height}}}$; the light‐green cross) located by computing the first derivative of $f{\left(x\right)}_{\text{height}}$ within the ${\mathit{\text{RGP}}}_{\mathit{\text{height}}}$ period. (b) The same algorithmic steps followed to analyze dynamic phenotypes for two indica and intermediate landraces. (c) The maximum canopy height (${\mathit{Max}}_{\mathit{\text{height}}}$; in cm), its associated DAS, and key growth stages such as the beginning of ripening estimated using maximum and minimum curvature values on a normalized‐curvature curve $f{\left(x\right)}_{\mathrm{cuv}}$ (dotted blue) derived from the $f{\left(x\right)}_{\text{height}}$ for three types of rice landraces.

Fig. 5

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.

Fig. 6

A series of 2D orthomosaics, pseudocoloured height maps and 3D point clouds collected by low‐cost unmanned aerial vehicles (UAVs) in the 2019 and 2020 seasons from 254 rice landraces in Shanghai. (a–c) 3D point clouds (from a 60° perspective) and overhead 2D orthomosaic of 254 landraces generated from a series of UAV phenotyping conducted over the 2019 season in Shanghai (to the left). Pseudocoloured height maps (to the right), showing plot‐based canopy plant height values for all the plots in the field. (d) Quantification of growth curves using AirMeasurer‐measured canopy height values for three types of rice landraces (i.e. indica, japonica and intermediary) over the 2019 season. Coloured shading areas denote confidence intervals (15^th–85^th percentiles). The three coloured dashed arrows indicate when the average maximum height values of the three types of landraces were reached (in days after sowing, DAS). (e–h) Experiments of the same 254 landraces repeated in the 2020 season, producing 3D point clouds, 2D orthomosaics, the height maps and derived growth curves. The unified height scale bar (0–150+ cm) for the subfigures is shown.

Fig. 7

Matrix generated to provide a comprehensive overview of the performance of 254 rice landraces in the 2019 season, through which dynamic analysis of normalized canopy coverage index (CCI), excess green (ExG) and visible atmospherically resistant index (VARI) were performed. (a) Eight 2D orthomosaics collected between 20 July and 8 October 2019 used to generate the performance matrix, where each cell was an example canopy image of a rice genotype, such that genotypes were columns and UAV phenotyping timepoints were rows. In the performance matrix, the 254 rice landraces rearranged according to three domestic types, i.e. indica (blue), japonica (red) and intermediary (green). (b–d) Using the matrix, dynamic analysis was performed to study traits such as CCI, ExG and VARI, demonstrating their different growth patterns and the time points when their maximum values were reached, e.g. ${\mathit{\text{Max}}}_{\mathit{\text{CCI}}}$ (85–95 d after sowing, DAS), ${\mathit{\text{Max}}}_{\mathit{\text{ExG}}}$ (73–79 DAS) and ${\mathit{\text{Max}}}_{\mathit{\text{VARI}}}$ (74–78 DAS).

Fig. 8

Coefficient of determination (R ²) computed to evaluate correlations between AirMeasurer‐derived and manually scored maximum plant height, normalized canopy coverage index (CCI) and normalized excess green index (ExG). The correlations between AirMeasurer‐derived, Gaussian‐fitted and manually measured canopy height values also are provided. (a) Plot‐based correlation between the maximum height measured by AirMeasurer and manual scoring using 177 rice landraces in the 2019 season. (b) Correlation between the maximum height trait measured by AirMeasurer and manual scoring using 254 landraces in the 2020 season. (c) Correlation between AirMeasurer‐derived canopy height values (based on calibrated 3D point clouds) and manual scoring of point cloud data to derive canopy height values using 177 landraces measured at eight time points (35–115 d after sowing (DAS)) across the 2019 season, 1416 plots in total; and (d) correlation between AirMeasurer‐derived canopy height values and Gaussian‐fitted values. (e) Correlation between normalized canopy coverage index (0–1) measured by AirMeasurer and the manually scored canopy area of plot images (in pixels) using 177 plots in 2019. (f) Correlation between normalized ExG index (0–1) measured by AirMeasurer and manually measured green values (0–255) using 177 plot images (35–115 DAS) in 2019.

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 ${\mathit{\text{Max}}}_{\mathit{\text{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 ${\mathit{\text{AGR}}}_{\mathit{\text{height}}}$ trait, between 0 d after sowing (DAS) and the ${\mathit{\text{Max}}}_{\mathit{\text{height}}}$ day, in 2020. (d) QTL for the ${\mathit{\text{AGR}}}_{\mathit{\text{height}}}$ trait (0 DAS – the ${\mathit{\text{FGR}}}_{\mathit{\text{height}}}$ day) in 2021. (e) Four loci associated with the ${\mathit{\text{RGP}}}_{\mathit{\text{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 ${\mathit{\text{AGR}}}_{\mathit{\text{CCI}}}$ in 2021, determined over the period between 0 DAS and the ${\mathit{\text{FGR}}}_{\mathit{\text{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 ${\mathit{\text{Max}}}_{\mathit{\text{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 ${\mathit{\text{AGR}}}_{\mathit{\text{ExG}}}$ trait for the interval 0 DAS – the ${\mathit{\text{Max}}}_{\mathit{\text{ExG}}}$ day. The major QTL co‐locates with SLB1/SLB2 and D61. Table 2 summarizes the QTLs associated with the above growth‐related traits. Abbreviations: maximum canopy height (${\mathit{\text{Max}}}_{\mathit{\text{height}}}$; cm), average growth rate for a target trait (${\mathrm{AGR}}_{\text{trait}}$; %), the fastest growth rate of canopy height (${\mathrm{FGR}}_{\text{height}}$; %), the rapid growth phase (${\mathit{\text{RGP}}}_{\mathit{\text{height}}}$; days), canopy coverage index (CCI), excess green (ExG), maximum CCI (${\mathit{\text{Max}}}_{\mathit{\text{CCI}}}$), maximum ExG (${\mathit{\text{Max}}}_{\mathrm{ExG}}$) and the fastest growth rate of CCI (${\mathit{\text{FGR}}}_{\mathit{\text{CCI}}}$; %).

Fig. 10

Manhattan plots and quantile‐quantile (QQ) plots for AirMeasurer‐derived traits subjected to a genome‐wide association study (GWAS) of 254 rice landraces trialled in 2019 and 2020. The significance threshold is shown by the horizontal grey dotted line. Known genes that co‐locate with significant loci are indicated by blue arrows. See Fig. 9 legend for trait abbreviations. (a) Manhattan plot and a QQ plot for the AirMeasurer‐derived ${\mathit{\text{Max}}}_{\mathit{\text{height}}}$ trait measured in 2019. The strongest signal on chromosome 1 was close to the sd1 gene and a strong signal on chromosome 3 was close to the OsHox32 gene. (b) Manhattan plot for the manually scored maximum plant height trait collected in the 2019 and 2020 seasons. (c) The 2019 ${\mathit{\text{AGR}}}_{\mathit{\text{height}}}$ (0 DAS – the ${\mathit{\text{Max}}}_{\mathit{\text{height}}}$ day) was used to identify four significant SNPs, co‐locating with known genes: sd1, OsHox32 (chromosome 3), NOG1 (chromosome 1) and OsGSK2 (chromosome 5). Analysis repeated using the same trait collected in 2020 and reproduced two SNPs, close to sd1 and OsGSK2. (d) GWAS performed with the trait ${\mathit{\text{AGR}}}_{\mathit{\text{height}}}$ (0 DAS – the ${\mathit{\text{FGR}}}_{\mathit{\text{height}}}$ day). Similar results were obtained in both seasons. (e) Plots for the 2020 ${\mathit{\text{AGR}}}_{\mathit{\text{CCI}}}$ (0 DAS – the ${\mathit{\text{Max}}}_{\mathit{\text{CCI}}}$ day) trait. Two signals were identified, one close to the Pit gene on chromosome 1 and the other near the PFPβ gene on chromosome 6. (f) In the analysis of the 2019 ${\mathit{\text{AGR}}}_{\mathit{\text{ExG}}}$ trait (0 DAS – the ${\mathit{\text{FGR}}}_{\mathit{\text{ExG}}}$ day), the strongest signal co‐located with the CCP1 gene on chromosome 1. Table 3 lists all of the significant association signals of the above growth‐related traits. DAS, days after sowing.