Assessing the Effect of Drought on Winter Wheat Growth Using Unmanned Aerial System (UAS)-Based Phenotyping

: Drought signiﬁcantly limits wheat productivity across the temporal and spatial domains. Unmanned Aerial Systems (UAS) has become an indispensable tool to collect reﬁned spatial and high temporal resolution imagery data. A 2-year ﬁeld study was conducted in 2018 and 2019 to determine the temporal effects of drought on canopy growth of winter wheat. Weekly UAS data were collected using red, green, and blue (RGB) and multispectral (MS) sensors over a yield trial consisting of 22 winter wheat cultivars in both irrigated and dryland environments. Raw-images were processed to compute canopy features such as canopy cover (CC) and canopy height (CH), and vegetation indices (VIs) such as Normalized Difference Vegetation Index (NDVI), Excess Green Index (ExG), and Normalized Difference Red-edge Index (NDRE). The drought was more severe in 2018 than in 2019 and the effects of growth differences across years and irrigation levels were visible in the UAS measurements. CC, CH, and VIs, measured during grain ﬁlling, were positively correlated with grain yield (r = 0.4–0.7, p < 0.05) in the dryland in both years. Yield was positively correlated with VIs in 2018 (r = 0.45–0.55, p < 0.05) in the irrigated environment, but the correlations were non-signiﬁcant in 2019 (r = 0.1 to − 0.4), except for CH. The study shows that high-throughput UAS data can be used to monitor the drought effects on wheat growth and productivity across the temporal and spatial domains. validation, UAS can be an important alternative tool to obtain CH measurements of wheat research plots. UAS-based canopy traits can be used to capture the response of various physiological attributes that cumulatively affect yield under drought. Integration of UAS-based remote sensing, agronomy, wheat breeding, and data analytics can potentially develop digital tools that can be used to select genotypes for drought tolerance and improve the genetic gain of wheat breeding programs.


Introduction
Wheat (Triticum aestivum L.) is the second most important crop in Texas, after cotton (Gossypium spp.), where it is produced under dryland and irrigated conditions for grain and forage. In 2020, it was planted on 1.9 million hectares, of which 0.97 million hectares were harvested for grain with a total production of 16 million Mt [1]. The frequent occurrence of drought poses a threat to winter wheat production in this region and reduces the yield under both the dryland and irrigated conditions [2]. The historic drought of 2011 and 2012 reduced wheat yield in Texas to 1.2 Mt compared to the five-year average of 2.5 Mt per hectare [3]. Another drought episode followed in 2018 when precipitation dropped below 150 mm. The frequent and unpredictable magnitude of drought observed this past decade reinforces the need to prepare for future climate threats.
Development of drought-tolerant cultivars along with better water management strategies can help reduce the impact of drought [4]. However, drought tolerance is a complex tion, stunted growth [41], and ultimately reduced yield. This study is a part of an ongoing project to integrate a UAS-based HTP system into the Texas A&M's winter wheat breeding program. The major objectives of this study are to (i) assess the impact of drought on winter wheat growth based on multi-temporal UAS-based canopy features, (ii) outline the procedure of UAS data collection, processing, and extraction, and (iii) determine the relationships between grain yield and UAS-based canopy features.  (Table 1). Agronomic practices such as nutrient and weed management were conducted when needed following the regional practices. The irrigated field received 343 mm and 95 mm irrigation in 2018 and 2019, respectively. Irrigation was applied with a three-span, lateralmove sprinkler system (Model 6000, Valmont Irrigation, Valley, NE, USA). The 2018 season was extremely dry and relatively hot compared to 2019 (Figures 2 and 3). Growing season precipitation in 2019 was 700 mm and 150 mm in 2018 (40 percentage of climate normal: calculated using the data from www.ncdc.noaa.gov, accessed on 1 January 2021).

UAS Data Collection
We obtained multispectral (MS) images by flying a DJI Matrice 100 quadcopter (SZ DJI Technology Co., Ltd., Shenzhen, China) equipped with an MS sensor: SlantRange 3P (SlantRange, Inc., California, USA). The SlantRange 3p has four 1.2-megapixel resolution sensors with central wavelengths of 550 nm (green), 650 nm (red), 720 nm (red-edge), and 850 nm (NIR). The spectral range covered by the green, red, red-edge, and NIR bands were 545-555 nm, 640-660 nm, 710-720 nm, and 840-860 nm, respectively. The SlantRange 3P sensor has a global shutter, and it has an onboard ambient illumination sensor that can be used to perform radiometric calibration in the post-processing stage. Therefore, data captured by the SlantRange 3P sensor throughout a single flight and from flight-to-flight during the growing season are comparable and consistent concerning changing weather conditions. Another platform, DJI Phantom 4 Pro equipped with a red, green, blue (RGB) sensor (SZ DJI Technology Co., Ltd., Shenzhen, China), was flown to collect RGB images. Phantom 4 Pro was flown at 15-20 m with 80-85% forward and side overlap to obtain 0.3-0.5 cm/pixel GSD (Ground Sampling Distance) when processed to generate orthomosaic images. Matrice 100 was flown at 30-35 m with a 70-75% overlap to get 1.2-1.7 cm/pixel GSD when processed to generate orthomosaic images. Nine semi-permanent ground control points (GCPs) were placed across the field after planting, and precise coordinates of the GCPs were surveyed using Trimble R8s Global Positioning System (GPS) (Trimble Inc., California, USA). Coordinates of the GCPs were used in the structure from motion (SfM) process to ensure accurate georeferencing of the UAS-based geospatial data products. Data were collected at 7 to 15 days intervals throughout the growing season ( Figure   4).

UAS Data Collection
We obtained multispectral (MS) images by flying a DJI Matrice 100 quadcopter (SZ DJI Technology Co., Ltd., Shenzhen, China) equipped with an MS sensor: SlantRange 3P (SlantRange, Inc., California, USA). The SlantRange 3p has four 1.2-megapixel resolution sensors with central wavelengths of 550 nm (green), 650 nm (red), 720 nm (red-edge), and 850 nm (NIR). The spectral range covered by the green, red, red-edge, and NIR bands were 545-555 nm, 640-660 nm, 710-720 nm, and 840-860 nm, respectively. The SlantRange 3P sensor has a global shutter, and it has an onboard ambient illumination sensor that can be used to perform radiometric calibration in the post-processing stage. Therefore, data captured by the SlantRange 3P sensor throughout a single flight and from flight-to-flight during the growing season are comparable and consistent concerning changing weather conditions. Another platform, DJI Phantom 4 Pro equipped with a red, green, blue (RGB) sensor (SZ DJI Technology Co., Ltd., Shenzhen, China), was flown to collect RGB images. Phantom 4 Pro was flown at 15-20 m with 80-85% forward and side overlap to obtain 0.3-0.5 cm/pixel GSD (Ground Sampling Distance) when processed to generate orthomosaic images. Matrice 100 was flown at 30-35 m with a 70-75% overlap to get 1.2-1.7 cm/pixel GSD when processed to generate orthomosaic images. Nine semi-permanent ground control points (GCPs) were placed across the field after planting, and precise coordinates of the GCPs were surveyed using Trimble R8s Global Positioning System (GPS) (Trimble Inc., California, USA). Coordinates of the GCPs were used in the structure from motion (SfM) process to ensure accurate georeferencing of the UAS-based geospatial data products. Data were collected at 7 to 15 days intervals throughout the growing season ( Figure   4).

UAS Data Collection
We obtained multispectral (MS) images by flying a DJI Matrice 100 quadcopter (SZ DJI Technology Co., Ltd., Shenzhen, China) equipped with an MS sensor: SlantRange 3P (SlantRange, Inc., San Diego, CA, USA). The SlantRange 3p has four 1.2-megapixel resolution sensors with central wavelengths of 550 nm (green), 650 nm (red), 720 nm (red-edge), and 850 nm (NIR). The spectral range covered by the green, red, red-edge, and NIR bands were 545-555 nm, 640-660 nm, 710-720 nm, and 840-860 nm, respectively. The SlantRange 3P sensor has a global shutter, and it has an onboard ambient illumination sensor that can be used to perform radiometric calibration in the post-processing stage. Therefore, data captured by the SlantRange 3P sensor throughout a single flight and from flight-to-flight during the growing season are comparable and consistent concerning changing weather conditions. Another platform, DJI Phantom 4 Pro equipped with a red, green, blue (RGB) sensor (SZ DJI Technology Co., Ltd., Shenzhen, China), was flown to collect RGB images. Phantom 4 Pro was flown at 15-20 m with 80-85% forward and side overlap to obtain 0.3-0.5 cm/pixel GSD (Ground Sampling Distance) when processed to generate orthomosaic images. Matrice 100 was flown at 30-35 m with a 70-75% overlap to get 1.2-1.7 cm/pixel GSD when processed to generate orthomosaic images. Nine semipermanent ground control points (GCPs) were placed across the field after planting, and precise coordinates of the GCPs were surveyed using Trimble R8s Global Positioning System (GPS) (Trimble Inc., Sunnyvale, CA, USA). Coordinates of the GCPs were used in the structure from motion (SfM) process to ensure accurate georeferencing of the UASbased geospatial data products. Data were collected at 7 to 15 days intervals throughout the growing season ( Figure 4).

UAS Data Processing
Raw images obtained from the MS sensor were first processed using SlantView software (SlantRange, Inc., San Diego, CA, USA) to acquire radiometrically calibrated MS images. These images, along with raw RGB images, were processed using Agisoft metashape software (Agisoft LLC, St. Petersburg, Russia, 191144) to generate geospatial data products such as a 3D point cloud, Digital Surface Model (DSM), and orthomosaic images (Figure 5). CC, ExG, NDVI, and NDRE were computed from orthomosaic images. Digital Terrain Model (DTM) of bare ground surface was generated from the images obtained after planning. Canopy Height Model (CHM) was obtained by subtracting DTM from DSM generated after every flight during the crop growth ( Figure 6). A plot boundary polygon layer was overlaid on CHM and measurement in 95th percentile height value within each polygon was extracted as the representative canopy height for each plot.

UAS Data Processing
Raw images obtained from the MS sensor were first processed using SlantView software (SlantRange, Inc., San Diego, CA, USA) to acquire radiometrically calibrated MS images. These images, along with raw RGB images, were processed using Agisoft metashape software (Agisoft LLC, St. Petersburg, Russia, 191144) to generate geospatial data products such as a 3D point cloud, Digital Surface Model (DSM), and orthomosaic images ( Figure 5). CC, ExG, NDVI, and NDRE were computed from orthomosaic images. Digital Terrain Model (DTM) of bare ground surface was generated from the images obtained after planning. Canopy Height Model (CHM) was obtained by subtracting DTM from DSM generated after every flight during the crop growth ( Figure 6). A plot boundary polygon layer was overlaid on CHM and measurement in 95th percentile height value within each polygon was extracted as the representative canopy height for each plot.

UAS Data Processing
Raw images obtained from the MS sensor were first processed using SlantView software (SlantRange, Inc., San Diego, CA, USA) to acquire radiometrically calibrated MS images. These images, along with raw RGB images, were processed using Agisoft metashape software (Agisoft LLC, St. Petersburg, Russia, 191144) to generate geospatial data products such as a 3D point cloud, Digital Surface Model (DSM), and orthomosaic images (Figure 5). CC, ExG, NDVI, and NDRE were computed from orthomosaic images. Digital Terrain Model (DTM) of bare ground surface was generated from the images obtained after planning. Canopy Height Model (CHM) was obtained by subtracting DTM from DSM generated after every flight during the crop growth ( Figure 6). A plot boundary polygon layer was overlaid on CHM and measurement in 95th percentile height value within each polygon was extracted as the representative canopy height for each plot.   CC was obtained by applying the Canopeo algorithm developed by Patrignani and Ochsner [42] to the RGB orthomosaics. In this approach, a certain threshold was used to classify the image into the canopy and non-canopy classes, and a binary image was produced by applying equation 1 during image processing.

, 2 , 2 3
(1) where R, G, and B are Digital Number (DN) values for the red, green, and blue bands, respectively. P1 and P2 are parameters to classify pixels in the green band and P3 sets the minimum value for ExG (Table 2) to select green vegetation. In this study, parameter (P1, P2, and P3) values were set to default as explained in the Canopeo algorithm (P1 = 0.95, P2 = 0.95, and P3 = 20) [42]. A plot boundary polygon layer was overlaid on the classified image and the percentage of canopy cover was calculated using Equation (2) (Figure 7).

Canopy Cover CC 100
(2) CC was obtained by applying the Canopeo algorithm developed by Patrignani and Ochsner [42] to the RGB orthomosaics. In this approach, a certain threshold was used to classify the image into the canopy and non-canopy classes, and a binary image was produced by applying equation 1 during image processing.
where R, G, and B are Digital Number (DN) values for the red, green, and blue bands, respectively. P1 and P2 are parameters to classify pixels in the green band and P3 sets the minimum value for ExG (Table 2) to select green vegetation. In this study, parameter (P1, P2, and P3) values were set to default as explained in the Canopeo algorithm (P1 = 0.95, P2 = 0.95, and P3 = 20) [42]. A plot boundary polygon layer was overlaid on the classified image and the percentage of canopy cover was calculated using Equation (2) (Figure 7).
Canopy Cover (CC) = Total number of canopy pixels Total number of pixels × 100 Three VIs were obtained, of which two indices, NDVI and NDRE, were obtained from the MS sensor, whereas ExG was obtained from the RGB sensor (Table 2). A plot boundary polygon layer was overlaid on the VIs map and the plot average value was extracted. Black denotes non-canopy pixels; white denotes canopy pixels. Plot boundaries denoted by red was overlaid on the classified image to calculate percentage of CC and extract data from each plot.
Three VIs were obtained, of which two indices, NDVI and NDRE, were obtained from the MS sensor, whereas ExG was obtained from the RGB sensor (Table 2). A plot boundary polygon layer was overlaid on the VIs map and the plot average value was extracted. Table 2. List of spectral vegetation indices. Reflectance obtained in red, green, blue, red-edge, and near-infrared region of the spectrum is denoted by R, G, B, RE, and NIR, respectively.

Ground Measurements
Manual plant height measurements were collected before harvest. A single measurement was taken per plot from the base of the plant to the tip of the spike. The heading date was recorded when 50% of spikes in a plot were fully exposed. In 2018, genotypes reached the heading date at 113-123 DOY and 119-127 DOY in the dryland and irrigated environments, respectively. Heading was much earlier than normal in 2018 due to drought stress. In 2019, it was in line with the long-term average; 122-129 DOY in dryland and 124-135 DOY in the irrigated environment. Each plot was harvested using a small plot combine (Classic Plus, Wintersteiger AG, Germany) to determine the final grain yield. Black denotes non-canopy pixels; white denotes canopy pixels. Plot boundaries denoted by red was overlaid on the classified image to calculate percentage of CC and extract data from each plot. Table 2. List of spectral vegetation indices. Reflectance obtained in red, green, blue, red-edge, and near-infrared region of the spectrum is denoted by R, G, B, RE, and NIR, respectively.

Vegetation Indices Formula References
Excess Green Index (ExG)

Ground Measurements
Manual plant height measurements were collected before harvest. A single measurement was taken per plot from the base of the plant to the tip of the spike. The heading date was recorded when 50% of spikes in a plot were fully exposed. In 2018, genotypes reached the heading date at 113-123 DOY and 119-127 DOY in the dryland and irrigated environments, respectively. Heading was much earlier than normal in 2018 due to drought stress. In 2019, it was in line with the long-term average; 122-129 DOY in dryland and 124-135 DOY in the irrigated environment. Each plot was harvested using a small plot combine (Classic Plus, Wintersteiger AG, Germany) to determine the final grain yield.

Data Analysis
Obtained UAS measurements were averaged across plots for each measurement date to compare the seasonal dynamics of canopy traits between the dryland and irrigated environment. Obtained measurements were plotted against the day of the year (DOY) to develop graphs for seasonal trends. DOY is the number representing the day between 1 Remote Sens. 2021, 13, 1144 8 of 20 and 365 of a year (DOY of January 1 is 1). Standard deviation was calculated across environments to show the variation in the UAS measurements. Data analysis was conducted using SAS version 9.4 (Statistical Analysis System Institute, Cary, NC, USA). PROC CORR procedure in SAS was used to determine Pearson correlation coefficients (r) between grain yield and UAS-based canopy traits. The coefficient of determination (R 2 ) was calculated to assess the relationship between ground-based plant height and UAS-based canopy height measurements.

Temporal Dynamics of CH
In general, CH began to increase from the first week of March (60 DOY) in 2018 ( Figure 8A) and the second week of March in 2019 (74 DOY) ( Figure 8B) with the increase in air temperature. Dryland CH was relatively low in both years once the growth started. The difference in CH between the two water levels was highly significant throughout the season in 2018. Extreme drought in 2018 resulted in reduced growth and the plants did not grow much more than 0.5 m tall. However, under irrigation, CH increase was linear and reached a maximum of 1.2 m by the second week of May (134 DOY). In 2019, under both environments, CH followed a similar trend as in the 2018 irrigated environment. CH reached maximum 10 days earlier (140 DOY) in dryland compared to irrigated environment (150 DOY). Average CH was 1.1 m and 0.75 m in irrigated and dryland environments, respectively ( Figure 8B). Significant differences in CH measurements between irrigated and dryland environments were observed after this point. CH measurements in the dryland environments were low in both years compared to the irrigated environments. In general, there was a sudden decline in CH after reaching the plateau except in 2018 dryland where the plants dried down early amid severe drought stress.
Obtained UAS measurements were averaged across plots for each measurement date to compare the seasonal dynamics of canopy traits between the dryland and irrigated environment. Obtained measurements were plotted against the day of the year (DOY) to develop graphs for seasonal trends. DOY is the number representing the day between 1 and 365 of a year (DOY of January 1 is 1). Standard deviation was calculated across environments to show the variation in the UAS measurements. Data analysis was conducted using SAS version 9.4 (Statistical Analysis System Institute, Cary, NC, USA). PROC CORR procedure in SAS was used to determine Pearson correlation coefficients (r) between grain yield and UAS-based canopy traits. The coefficient of determination (R 2 ) was calculated to assess the relationship between ground-based plant height and UAS-based canopy height measurements.

Temporal Dynamics of CH
In general, CH began to increase from the first week of March (60 DOY Figure 8B). Significant differences in CH measurements between irrigated and dryland environments were observed after this point. CH measurements in the dryland environments were low in both years compared to the irrigated environments. In general, there was a sudden decline in CH after reaching the plateau except in 2018 dryland where the plants dried down early amid severe drought stress. CH measured from UAS were validated with ground-based manual measurements ( Figure 9). A linear relationship was found between UAS-obtained CH measurements and ground-based measurements. The R 2 values ranged from 0.5 to 0.7 depending on the environment. Hassan et al. [36] found a strong correlation (R 2 = 0.8-0.9) between UAS and ground-based plant height data measured at the booting and mid-grainfilling stages. Relationships were relatively strong under irrigated condition in both years compared to the CH measured from UAS were validated with ground-based manual measurements ( Figure 9). A linear relationship was found between UAS-obtained CH measurements and ground-based measurements. The R 2 values ranged from 0.5 to 0.7 depending on the environment. Hassan et al. [36] found a strong correlation (R 2 = 0.8-0.9) between UAS and ground-based plant height data measured at the booting and mid-grainfilling stages. Relationships were relatively strong under irrigated condition in both years compared to the dryland. Lower correlations are due to the time at which manual measurements were taken. Manual measurements were taken from a single plant per plot after crop maturity when there are high chances of lodging and disturbance in the canopy structure. Underestimation of CH by UAS as seen in some previous studies [45] was found in our study across all environments. Considering the seasonal trend of CH and accuracy of data dryland. Lower correlations are due to the time at which manual measurements w taken. Manual measurements were taken from a single plant per plot after crop matur when there are high chances of lodging and disturbance in the canopy structure. Und estimation of CH by UAS as seen in some previous studies [45] was found in our stu across all environments. Considering the seasonal trend of CH and accuracy of data va dation, UAS can be an important alternative tool to obtain CH measurements of wh research plots.

Temporal Dynamics of CC
In 2018 canopy growth started early in the season ( Figure 10A). The CC was alrea 40-50% in the first week of January and remained in the same range until the end of Ja uary (31 DOY). Low temperature caused a dip in CC values in the first week of Februa (38 DOY). CC values gradually increased after reduced winter growth in both enviro ments and reached 90-100% in the irrigated environment and 45-55% in the dryland. Ca opy growth started late in 2019 ( Figure 10B). CC was less than 5% until the last week February (57 DOY) when it started to grow rapidly and reached a maximum of 100% a 88 % in the first week of April (98 DOY) under irrigated and dryland environments, spectively. CC remained the same until the third week of May (139 DOY) and started decrease gradually and reached 1% in the first week of June (168 DOY) under irrigati In the dryland environment, the reduction in CC values started a week earlier (133 DO and reached 1% by 161 DOY. CC measurements clearly depicted the environmental co ditions. In the severe drought year of 2018, peak CC measured in the irrigated enviro ment was double that measured in the dryland environment. Significant differences in C were also found in 2019, but the difference was considerably less than in 2018. Cano growth followed a sigmoidal growth function in 2019 with slow early growth in the ginning followed by a linear growth phase, and then a steady stage. The decline in C after the steady stage is attributed to senescence. The changes in the color of the cano across different environments are shown in Figure 4. The trend of CC aligns well with

Temporal Dynamics of CC
In 2018 canopy growth started early in the season ( Figure 10A). The CC was already 40-50% in the first week of January and remained in the same range until the end of January (31 DOY). Low temperature caused a dip in CC values in the first week of February (38 DOY). CC values gradually increased after reduced winter growth in both environments and reached 90-100% in the irrigated environment and 45-55% in the dryland. Canopy growth started late in 2019 ( Figure 10B). CC was less than 5% until the last week of February (57 DOY) when it started to grow rapidly and reached a maximum of 100% and 88 % in the first week of April (98 DOY) under irrigated and dryland environments, respectively. CC remained the same until the third week of May (139 DOY) and started to decrease gradually and reached 1% in the first week of June (168 DOY) under irrigation. In the dryland environment, the reduction in CC values started a week earlier (133 DOY) and reached 1% by 161 DOY. CC measurements clearly depicted the environmental conditions. In the severe drought year of 2018, peak CC measured in the irrigated environment was double that measured in the dryland environment. Significant differences in CC were also found in 2019, but the difference was considerably less than in 2018. Canopy growth followed a sigmoidal growth function in 2019 with slow early growth in the beginning followed by a linear growth phase, and then a steady stage. The decline in CC after the steady stage is attributed to senescence. The changes in the color of the canopy across different environments are shown in Figure 4. The trend of CC aligns well with the seasonal change in ground cover and reached maximum when the ground was fully covered by the canopy. This shows the usefulness of estimating CC using UAS to calculate ground cover and assess the effect of drought on vegetation growth.
Remote Sens. 2021, 13, x FOR PEER REVIEW 10 of 21 seasonal change in ground cover and reached maximum when the ground was fully covered by the canopy. This shows the usefulness of estimating CC using UAS to calculate ground cover and assess the effect of drought on vegetation growth.

Temporal Dynamics of VIs
The temporal dynamics of NDVI, NDRE, and ExG obtained from UAS during the wheat growing season is shown in Figure 11, Figure 12, and Figure 13, respectively. General growth pattern of these indices was similar to CC across all the environments. The pattern can be divided into four stages based on the greenness of the canopy; germination and stand establishment, accelerated vegetative growth accompanied by tillering, steady stage when the ground was fully covered by the canopy, and a declining stage during physiological maturity.
NDVI was obtained from MS sensor and uses the visible and NIR light. This index has been widely used to monitor vegetation dynamics and crop health. Drought significantly affected NDVI measurements in our study. The effect was proportionate to the intensity of drought as seasonal NDVI values were low in dryland condition compared to irrigated condition in 2018 ( Figure 11A). As the experimental site received frequent precipitation in 2019, dryland NDVI measurements were similar to irrigated measurements (

Temporal Dynamics of VIs
The temporal dynamics of NDVI, NDRE, and ExG obtained from UAS during the wheat growing season is shown in Figure 11, Figure 12, and Figure 13, respectively. General growth pattern of these indices was similar to CC across all the environments. The pattern can be divided into four stages based on the greenness of the canopy; germination and stand establishment, accelerated vegetative growth accompanied by tillering, steady stage when the ground was fully covered by the canopy, and a declining stage during physiological maturity.

Temporal Dynamics of VIs
The temporal dynamics of NDVI, NDRE, and ExG obtained from UAS during the wheat growing season is shown in Figure 11, Figure 12, and Figure 13, respectively. General growth pattern of these indices was similar to CC across all the environments. The pattern can be divided into four stages based on the greenness of the canopy; germination and stand establishment, accelerated vegetative growth accompanied by tillering, steady stage when the ground was fully covered by the canopy, and a declining stage during physiological maturity.
NDVI was obtained from MS sensor and uses the visible and NIR light. This index has been widely used to monitor vegetation dynamics and crop health. Drought significantly affected NDVI measurements in our study. The effect was proportionate to the intensity of drought as seasonal NDVI values were low in dryland condition compared to irrigated condition in 2018 ( Figure 11A). As the experimental site received frequent precipitation in 2019, dryland NDVI measurements were similar to irrigated measurements (   NDVI was obtained from MS sensor and uses the visible and NIR light. This index has been widely used to monitor vegetation dynamics and crop health. Drought significantly affected NDVI measurements in our study. The effect was proportionate to the intensity of drought as seasonal NDVI values were low in dryland condition compared to irrigated condition in 2018 ( Figure 11A). As the experimental site received frequent precipitation in 2019, dryland NDVI measurements were similar to irrigated measurements ( Figure 11B) until heading (135 DOY). The difference was significant once the heading was completed, and plants began to senescence. The decline in NDVI during senescence was slow under  NDRE was also calculated from radiometrically calibrated MS orthomosaics. It uses the NIR and red-edge light which is the transition region between the visual red and NIR light. NDVI is prone to early saturation in dense vegetation [46] and NDRE can serve as an alternative to NDVI to assess spatial and temporal variabilities [47,48]. In our study, both NDVI and NDRE curves closely resemble each other. NDRE values ranged from 0.1 to 0.8 in all environments except in 2018 dryland in which the values ranged from 0.1 to 0.5. The effect of drought was severe in 2018 compared to 2019 as NDRE measurements were relatively low throughout the growing season ( Figure 12A). Dryland and irrigated NDRE measurements were about the same in 2019 and differences were not significant until after heading (135 DOY) ( Figure 12B). ExG was calculated utilizing the broad band lights from RGB sensor. Compared to MS-based indices, the trends were similar with early slow increase followed by linear growth, and decline at later growth stage, but there were higher oscillations across all environments ( Figure 13). In 2018, under irrigated conditions, ExG reached a maximum at 99 DOY and had a dip on 105 and 113 DOY ( Figure 13A). It again reached a plateau on 124 DOY and then started to decline. The dryland ExG never exceeded 0.1 in the dryland. In 2019, ExG values ranged from 0.1 to 0.45 in irrigated condition and the values in dryland ranged from 0.1 to 0.35. The drop in ExG values were highly proportionate to the intensity of drought. In 2019, ExG values increased consistently but after attaining the plateau inconsistent patterns were noted. In some DOYs, values were high while in others there were dips ( Figure 13B). The fluctuations happened when CC, NDVI, and NDRE remained steady during the 90-100% ground cover. ExG was highly sensitive to weather variations and change in the color of the canopy (Figure 4).  ExG was calculated utilizing the broad band lights from RGB sensor. Compared to MS-based indices, the trends were similar with early slow increase followed by linear growth, and decline at later growth stage, but there were higher oscillations across all environments ( Figure 13 Figure 13B). The fluctuations happened when CC, NDVI, and NDRE remained steady during the 90-100% ground cover. ExG was highly sensitive to weather variations and change in the color of the canopy (Figure 4). NDRE was also calculated from radiometrically calibrated MS orthomosaics. It uses the NIR and red-edge light which is the transition region between the visual red and NIR light. NDVI is prone to early saturation in dense vegetation [46] and NDRE can serve as an alternative to NDVI to assess spatial and temporal variabilities [47,48]. In our study, both NDVI and NDRE curves closely resemble each other. NDRE values ranged from 0.1 to 0.8 in all environments except in 2018 dryland in which the values ranged from 0.1 to 0.5. The effect of drought was severe in 2018 compared to 2019 as NDRE measurements were relatively low throughout the growing season ( Figure 12A). Dryland and irrigated NDRE measurements were about the same in 2019 and differences were not significant until after heading (135 DOY) ( Figure 12B).
ExG was calculated utilizing the broad band lights from RGB sensor. Compared to MS-based indices, the trends were similar with early slow increase followed by linear growth, and decline at later growth stage, but there were higher oscillations across all environments ( Figure 13). In 2018, under irrigated conditions, ExG reached a maximum at 99 DOY and had a dip on 105 and 113 DOY ( Figure 13A). It again reached a plateau on 124 DOY and then started to decline. The dryland ExG never exceeded 0.1 in the dryland. In 2019, ExG values ranged from 0.1 to 0.45 in irrigated condition and the values in dryland ranged from 0.1 to 0.35. The drop in ExG values were highly proportionate to the intensity of drought. In 2019, ExG values increased consistently but after attaining the plateau inconsistent patterns were noted. In some DOYs, values were high while in others there were dips ( Figure 13B). The fluctuations happened when CC, NDVI, and NDRE remained steady during the 90-100% ground cover. ExG was highly sensitive to weather variations and change in the color of the canopy (Figure 4).

Correlations Between Grain Yield and UAS-Based Parameters
There was a significant (p < 0.05) difference in grain yield between dryland and irrigated environments in both years. The Yield gap between the two environments was

Correlations Between Grain Yield and UAS-Based Parameters
There was a significant (p < 0.05) difference in grain yield between dryland and irrigated environments in both years. The Yield gap between the two environments was  In general, the correlation of grain yield with the UAS-based canopy traits obtained early in the growing season was weak compared to traits obtained after heading through maturity (Table 3). In 2018 dryland, CC, CH, and EXG had lower correlation coefficients (r) compared to NDRE and NDVI most of the times during the growing season. Significant (p < 0.05) r values greater than 0.3 for CC and ExG were on 87-97 DOY (boot stage) and 117-128 DOY (heading, anthesis and grain filling). Similarly, CH had r values greater than 0.3 during 117 DOY. NDRE and NDVI had r significant and greater than 0.3 during 99 DOY to 134 DOY (boot, heading, anthesis, and grain filling stage). The highest r value found in this environment was 0.58 (p < 0.05) from NDRE during grain filling period. CC, CH, and ExG were obtained from RGB images while NDVI and NDRE were obtained from radiometrically calibrated MS images. Plots in 2018 had early ground cover compared to 2019 ( Figure 10A) but this early biomass accumulation was not related to final grain yield. Lack of seasonal precipitation along with the reduction of available soil water during heading, anthesis, and grain filling might be the reason for not having higher correlation between yield and UAS-based canopy features measured early in the season. The 2018 irrigated environment was better with higher r values of canopy traits with grain yield. CC and ExG performed better than NDRE and NDVI with higher r values. CC had r values greater than 0.3 most of the times except 23-38 DOY (early tillering stage) and 99 DOY. CH had higher r values at 134-150 DOY (after heading). Higher r values of ExG (r > 0.45, p < 0.05) was during 105 DOY to 150 DOY (pre-heading, heading, anthesis, and grain filling stage). NDRE and NDVI also had higher r values multiple times during the growing season. The highest r value found in 2018 irrigated environment was 0.57 (p < 0.05) from ExG measured during the grain filling period (134 DOY).
The 2019 dryland environment was favored by frequent rainfall events. Additionally, early season canopy coverage was low compared to 2018. CC, CH, ExG, NDRE, and NDVI In general, the correlation of grain yield with the UAS-based canopy traits obtained early in the growing season was weak compared to traits obtained after heading through maturity (Table 3). In 2018 dryland, CC, CH, and EXG had lower correlation coefficients (r) compared to NDRE and NDVI most of the times during the growing season. Significant (p < 0.05) r values greater than 0.3 for CC and ExG were on 87-97 DOY (boot stage) and 117-128 DOY (heading, anthesis and grain filling). Similarly, CH had r values greater than 0.3 during 117 DOY. NDRE and NDVI had r significant and greater than 0.3 during 99 DOY to 134 DOY (boot, heading, anthesis, and grain filling stage). The highest r value found in this environment was 0.58 (p < 0.05) from NDRE during grain filling period. CC, CH, and ExG were obtained from RGB images while NDVI and NDRE were obtained from radiometrically calibrated MS images. Plots in 2018 had early ground cover compared to 2019 ( Figure 10A) but this early biomass accumulation was not related to final grain yield. Lack of seasonal precipitation along with the reduction of available soil water during heading, anthesis, and grain filling might be the reason for not having higher correlation between yield and UAS-based canopy features measured early in the season. The 2018 irrigated environment was better with higher r values of canopy traits with grain yield. CC and ExG performed better than NDRE and NDVI with higher r values. CC had r values greater than 0.3 most of the times except 23-38 DOY (early tillering stage) and 99 DOY. CH had higher r values at 134-150 DOY (after heading). Higher r values of ExG (r > 0.45, p < 0.05) was during 105 DOY to 150 DOY (pre-heading, heading, anthesis, and grain filling stage). NDRE and NDVI also had higher r values multiple times during the growing season. The highest r value found in 2018 irrigated environment was 0.57 (p < 0.05) from ExG measured during the grain filling period (134 DOY). The 2019 dryland environment was favored by frequent rainfall events. Additionally, early season canopy coverage was low compared to 2018. CC, CH, ExG, NDRE, and NDVI measured until 126 DOY had non-significant correlation (p > 0.05) with yield (Table 3) period. The 2019 irrigated season showed a different phenomenon regarding the correlation between VIs and grain yield. Except CH, correlation was either negative (r = 0.14-0.40) or non-significant (p > 0.05). ExG had highest r values in 2018 irrigated condition but in 2019 the correlation was negative. A similar pattern was shown by CC, NDRE, and NDVI. Plots attained maximum ground cover in 94 DOY and maintained 100 percent canopy coverage until 144 DOY. Senescence in 2019 started late and was slower than in 2018. The variability among the plots was low with respect to change in canopy cover which might have contributed to lower correlation values in the 2019 irrigated environment. CH measured at 121-133 DOY, however, had better r values (r > 0.3, p < 0.05).

Growth Dynamics based on UAS-Based Canopy Traits
All the canopy traits: CC, NDVI, NDRE, and ExG, followed general growth characteristics of winter wheat. These traits started with slow growth at the beginning of the season (germination and early stand establishment), and as the season progressed (tillering), showed linear growth (stem elongation along with canopy growth), and reached a plateau of maximum ground cover. The plateau lasted until heading which was again followed by a slow decline until maturity. The decrease in VIs and CC after heading is attributed to the change in color of the canopy because of senescence and physiological maturity [49]. A similar growth pattern of canopy attributes and vegetation indices was observed in other studies [36,[50][51][52][53][54]. Although the growth trend was similar to these studies, the time of event differed in our study. For example, Lin et al. [55] recorded highest NDVI value during the heading stage of wheat but maximum CC, NDVI, and NDRE were achieved 10-15 days before heading in our study. This was mainly because high tillering and increased vegetation growth led to early canopy saturation. Similar results were observed by Song and Wang [56] in a phenology study conducted using Sentinel-1 data. In the case of plant height, except for the 2018 dryland, CH followed a sigmoidal growth pattern and reached maximum during heading followed by slight reduction once it reached its maximum. Anderson et al. [57] also found a similar trend in seasonal plant height measured by UAS in maize. There was some lodging in our study which reduced the average value in the late growing season. Although the nature of the trend was similar in 2018 and 2019, some canopy traits exhibited difference in the time and duration of phenological changes. For example, all the traits showed higher values between 0 to 40 DOY in 2018 compared to same period in 2019 which was due to early planting in 2018 complemented by moisture stored in the soil due to precipitation before planting in 2017.
Combining the seasonal trend of CC and CH obtained in our study might be useful to assess the phenology of winter wheat growth. The growth dynamics can be divided into four different stages by looking at the CC and CH curves. Slow early growth with an increase in ground cover during the tillering stage, the onset and duration of linear growth can be classified as stem elongation period. The point at which CH reached maximum and CC starts to decline (150 DOY) happened after heading and at the onset of flowering. The duration of the decline in CC can be classified as a ripening stage. Slow early growth can be used to assess early stand establishment and vigor [58]. Linear growth in CH and CC can be used to assess maximum canopy growth and duration to achieve and maintain maximum ground cover. The decline of CC can be used to assess senescence rate [50,52], and stay green [59,60] phenotypes. Another approach to differentiate phenological stages is by using VIs. Magney et al. [61] used spectral reflectance radiometers to collect daily NDVI measurements to monitor crop phenology of soft white spring wheat. Four phenological growth stages namely, tillering, stem elongation, heading, and ripening were derived by fitting a non-parametric regression curve on the daily measurements. Zheng et al. [62] used first order derivative of time series NDVI measurements to estimate the time of tillering, heading, and maturity in rice. Chu et al. [55] used MODIS time series data and fitted double Gaussian function to detect green-up, heading, and harvesting phases of wheat. Successful detection of phenological phases can be beneficial to wheat breeders to screen genotypes for respective traits. Early detection of these phases can be helpful to assess the future drought impact on canopy growth and can support farmers to make management decisions. Although it needs further research and validation, this study showed the possibility of using multi-temporal UAS data collected over the entire season to estimate the time of phenological phases in winter wheat.
Several physiological and biochemical traits [63][64][65][66][67] have been evaluated and proposed to screen genotypes for drought tolerance in wheat. However, in traditional breeding programs grain yield is the primary trait used to make selection decisions [68]. The spatial and temporal constraints of measuring physiological traits in many genotypes limits their use in plant breeding. Methodologies used in this study can be useful to obtain data of yield-associated traits that can be used to identify genotypes that can resist, tolerate, or adapt under water-limited conditions. The ability to track the seasonal dynamics of these traits can help to understand the time, intensity, and duration of the impact of drought on canopy growth. Additionally, integration of different traits can help to dissect the effect of drought according to key phenological stages [67]. For example, CC can be used to estimate early ground cover and NDRE can be used to estimate senescence rate and calculate stay-green duration, thus optimizing the utilization of available soil water [68].
MS sensor-based indices were more stable than ExG (obtained from RGB sensor). Weekly ExG measurements had a higher standard deviation on individual dates and were more sensitive to the changes in canopy color and weather conditions. A notable difference was seen during low winter temperatures and recovery in CC (60-60 DOY, Figure 10A). Radiometric calibration of the MS sensor can be attributed to its stable measurements [69] compared to the RGB sensor.
Drought significantly affected canopy growth. Differences in the growth curves of UAS-based canopy traits between the dryland and irrigated condition clearly illustrated the impact of drought on canopy growth, which other studies also found [36,[70][71][72]. Dryland CC values in 2018 started off strong but lagged as the season progressed due to the extreme drought in 2018. In 2019, the increase in CC started once irrigation was applied and CC in the irrigated environment remained higher throughout the season. The difference between dryland and irrigated environment was larger after heading which is indicative of moisture stress during the reproductive stage [60] as dryland plants could not maintain greenness. Drought significantly reduces leaf water potential, relative leaf water content [41], chlorophyll content [73], stomatal conductance, and photosynthesis [74] resulting in canopy wilting, stunted growth, and leaf area reduction [41]. Positive correlations exist between VIs and physiological traits such as leaf equivalent water thickness, canopy water content (CWC) [75] and chlorophyll content [76]. The difference in UAS derived canopy traits between dryland and the irrigated environment was clearly reflected in respective yield differences in both years. Thus, the seasonal patterns of UAS-based canopy features can enhance our understanding of the effect of drought and crop management practices on vegetation growth [77].

Association between UAS-Based Canopy Traits and Grain Yield
A consistent positive correlation between UAS-based canopy traits and grain yield was found under the dryland environment in both years. The relationship between spectral VIs and grain yield under dryland has been established in several previous studies [78][79][80][81]. Canopy features and vegetation indices determined by the amount and intensity of the green canopy has been used to estimate above-ground biomass in wheat [82][83][84] and the biomass is correlated with grain yield [85]. This leads to the association of canopy features with grain yield. Naser et al. [49] found a positive association between wheat grain yield and NDVI measured at anthesis, and during grain filling period under dryland environment. Thapa et al. [85] also found a significant correlation between yield and NDVI measured after the jointing stage of wheat. However, they cautioned about using NDVI for screening genotypes for yield under extreme drought conditions. Our 2018 dryland environment went through extreme drought with less than 150 mm precipitation during the entire growing season. Even with the extreme drought in 2018, a significant positive correlation was observed between yield and UAS-based canopy features measured during anthesis to maturity. Our results showed that VIs generated from MS sensor yielded higher correlation coefficients compared to ones obtained from the RGB sensor. In contrast to the 2018 experiment, 2019 was moderate in drought, and the correlation between grain yield and UAS-based features improved significantly.
A significant positive relationship between yield and each of CC and VIs was found in the 2018-irrigated environment. Like dryland environment, the UAS-based features obtained after heading had stronger correlations with grain yield. This result agrees with Thapa et al. [85] in which they showed a positive relationship between yield and Green Seeker-based NDVI data obtained in the same location using different sets of genotypes in the same year. Positive correlation between spectral indices and grain yield under irrigated environment was found in several other studies in wheat [36]. However, in 2019, the nature of this correlation was different. Except for CH, the relationships were either nonsignificant or negative. Genotypes in 2019 were not water-stressed as the available water was higher than normally required during the growing season. Early canopy saturation reduced the variability in UAS-based canopy features within the irrigated environment as plants attained maximum canopy size and maintained greenness which might be the reason for not obtaining a similar correlation trend as in other environments. The uniformity in the greenness of the canopy might be the reason for reduced ability of the canopy traits to differentiate genotypes with respect to yield. Naser et al. [49] also found a poor association of NDVI with wheat grain yield in irrigated environments and suggested that NDVI can perform a better job when the sensor is not saturated. The results obtained by Casadesus et al. [78] also support the findings of our study as the NDVI measured during anthesis of wheat genotypes was not able to provide better yield estimation because of complete soil coverage and high plant density. Four irrigation events on 66, 98, 107, and 136 DOY, along with higher precipitation, caused higher vegetative growth and excessive biomass accumulation in our study. It can be inferred that the CC and VIs used in this study are not effective to estimate yield in winter wheat with excessive biomass.
Finding optimal time for collecting UAS data to assess drought tolerance with respect to yield is important to plant breeders to reduce the cost of data collection and processing. In our study, we collected UAS data on a weekly basis, almost 20 times during the season. Even in a moderate drought year like 2019, RGB data collected after 100 DOY and MS data collected after 130 DOY were valuable to find the impact of drought in the canopy. Additionally, the yield was highly correlated with UAS-based canopy features measured after heading as in other studies. Hassan et al. [36] found lower correlation coefficients of yield with VIs measured during early growth stages like stem elongation, booting, and heading. Selecting cultivars based on the measurements obtained after heading especially, during anthesis and grain filling period [49,85,86], seem efficient, but our results suggest that taking multiple measurements across several growth stages can add significant benefits to the screening process.

Conclusions
This study demonstrated the collection, processing, and phenotypic feature extraction procedure of UAS. Experimental results showed that the multi-temporal UAS data is useful for monitoring winter wheat growth and assessing the impact of drought on canopy growth. UAS-based canopy traits were affected by the intensity of drought on wheat vegetation. Besides, the integration of season-long multi-temporal canopy features and vegetation indices showed the potential to identify the phenological growth stages and assess the impact of drought at these different stages. Traits obtained from the MS sensor were more stable than from the RGB sensor since radiometric calibration improved the stability of measurements across time. UAS-based traits obtained from heading to maturity correlated with grain yield. However, the relationship largely depends on the growing conditions and may not always predict yield, especially when the canopies are dense. Multi-temporal UAS-based canopy traits can be used to capture the response of various physiological attributes that cumulatively affect yield under drought. Integration of UAS-based remote sensing, agronomy, wheat breeding, and data analytics can potentially develop digital tools that can be used to select genotypes for drought tolerance and improve the genetic gain of wheat breeding programs.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author. The data are not publicly available due to privacy.