PhenoCams for Field Phenotyping: Using Very High Temporal Resolution Digital Repeated Photography to Investigate Interactions of Growth, Phenology, and Harvest Traits

Understanding the interaction of plant growth with environmental conditions is crucial to increase the resilience of current cropping systems to a changing climate. Here, we investigate PhenoCams as a high-throughput approach for field phenotyping experiments to assess growth dynamics of many different genotypes simultaneously in high temporal (daily) resolution. First, we develop a method that extracts a daily phenological signal that is normalized for the different viewing geometries of the pixels within the images. Second, we investigate the extraction of the in season traits of early vigor, leaf area index (LAI), and senescence dynamic from images of a soybean (Glycine max) field phenotyping experiment and show that it is possible to rate early vigor, senescence dynamics, and track the LAI development between LAI 1 and 4.5. Third, we identify the start of green up, green peak, senescence peak, and end of senescence in the phenological signal. Fourth, we extract the timing of these points and show how this information can be used to assess the impact of phenology on harvest traits (yield, thousand kernel weight, and oil content). The results demonstrate that PhenoCams can track growth dynamics and fill the gap of high temporal monitoring in field phenotyping experiments.

Understanding the interaction of plant growth with environmental conditions is crucial to increase the resilience of current cropping systems to a changing climate. Here, we investigate PhenoCams as a high-throughput approach for field phenotyping experiments to assess growth dynamics of many different genotypes simultaneously in high temporal (daily) resolution. First, we develop a method that extracts a daily phenological signal that is normalized for the different viewing geometries of the pixels within the images. Second, we investigate the extraction of the in season traits of early vigor, leaf area index (LAI), and senescence dynamic from images of a soybean (Glycine max) field phenotyping experiment and show that it is possible to rate early vigor, senescence dynamics, and track the LAI development between LAI 1 and 4.5. Third, we identify the start of green up, green peak, senescence peak, and end of senescence in the phenological signal. Fourth, we extract the timing of these points and show how this information can be used to assess the impact of phenology on harvest traits (yield, thousand kernel weight, and oil content). The results demonstrate that PhenoCams can track growth dynamics and fill the gap of high temporal monitoring in field phenotyping experiments.

INTRODUCTION
Understanding the interaction of growth dynamics and phenology with environmental conditions is crucial to improve the resilience of our cropping systems to a changing climate (Pretty et al., 2010;Asseng et al., 2013Asseng et al., , 2015Siebert and Ewert, 2014). One approach to this is to characterize the growth of many different genotypes in different environments. Phenomics or phenotyping aims to describe the appearance of plants as a function of the interaction of its genetic background with environmental conditions and is currently one of the most rapidly developing disciplines in crop science (Fiorani and Schurr, 2013;Walter et al., 2015). New platforms and imaging techniques (Li et al., 2014;Cendrero-Mateo et al., 2017;Hund et al., 2019) aim to facilitate high-throughput analyses of plant traits. For field phenotyping applications, where plants need to be monitored in outdoor environments, the development of standardized methods for plant trait mapping is crucially needed (White et al., 2012;Busemeyer et al., 2013;Deery et al., 2014;Kirchgessner et al., 2017;Virlet et al., 2017). This setting, referred to as field phenotyping, has been identified as one of the largest challenges in plant phenotyping (Pieruschka, 2016;Araus et al., 2018).
Today, the most common technique to assess the growth of many different varieties under field conditions is still the visual rating by trained experts. In recent years, large infrastructures have been put into practice that can screen crop traits in a (semi-) automatic fashion Virlet et al., 2017). Also unmanned aerial systems (Bendig et al., 2015;Liebisch et al., 2015;Zaman-Allah et al., 2015;Gómez-Candón et al., 2016;Haghighattalab et al., 2016;Kefauver et al., 2017;Madec et al., 2017;Yang et al., 2017;Aasen and Bolten, 2018;Hu et al., 2018;Roth et al., 2018) and manual and (semi-) automated ground-based approaches (Gnyp et al., 2013;Tilly et al., 2015;Madec et al., 2017;Jimenez-Berni et al., 2018) are now used for this purpose. All of these approaches can deliver very high spatial resolution data. However, dynamic processes such as phenology might only differ by a couple of days between different genotypes of the same crop species or under different environmental conditions (e.g., Anderegg et al., 2020). Thus, methods that can deliver data at very high temporal resolution (at a daily rate) of many experimental plots at the same time are needed in the context of plant breeding and variety testing.
In the field of ecology, the technique of digital repeat photography-also called PhenoCams-that continuously capture images of a given area with an RGB or near-infraredenabled cameras has been used for more than a decade to estimate phenology (Richardson et al., 2007(Richardson et al., , 2009Ahrends et al., 2009;Penuelas et al., 2009;Graham et al., 2010;Ide and Oguma, 2010;Kurc and Benton, 2010;Migliavacca et al., 2011;Keenan and Richardson, 2015;Hufkens et al., 2016). Usually, the cameras are mounted on towers or lookouts and view the canopy horizontally or obliquely to record the objects within their field of view several times a day. Within these images, the brightness information of the RGB channels can then be used to track changes in phenology. Based on the digital numbers (DN) of the RGB channels, different color indices (CIs) such as the green or red chromatic coordinates or the excess green index were evaluated for that purpose. These CIs normalize the scene illumination or have an improved capability to distinguish vegetation . Based on CIs derived with Phenocams, many studies have investigated the timing of phenological events [also called phenological transition dates (Richardson, 2019)], such as budburst and leaf-out (Klosterman et al., 2014) or the start and end of the season, as well as the peak of redness and end of redness during senescence (Xie et al., 2018). Different software packages have been published to process PhenoCam data. PHENOR in combination with the PHENOCAM and DAMETR R package allows extraction and modeling of phenological information from images. Hufkens et al. (2018a) showed how this data could be used to evaluate phenology models. The Phenopix R package automates different fitting techniques as well as the extraction of phenological events . Since different terminologies are used throughout the literature, we will use the term "PhenoTimePoint" as a point in time that refers to an event within the phenological signal and the term "PhenoPhase" as a temporal period during the phenological development.
Still, not too many studies have been published for agricultural crops. Sakamoto et al. (2012) used two compact digital cameras to capture visible and near-infrared images to observe seasonal changes in crop growth of maize and soybean (Glycine max) during one year. They found that the camera-derived CIs were closely related to VIs calculated using a multispectral sensor (SKYE) and MODIS satellite reflectance. Additionally, they found that CIs were related to the change of green LAI, total LAI, and above-ground dry biomass of stalks and leaves during the season. Zhu et al. (2016) used a tower-based system to detect the wheat heading stage based on ear feature detection within the images and Brocks and Bareth (2018) used stereo images from a tower based system to estimate biomass in barley.
To the authors' knowledge, the PhenoCam approach has not yet been used (i) in the typical field phenotyping setting, where the growth of multiple genotypes needs to be characterized simultaneously, and (ii) to investigate the impact of the timing and duration of PhenoPhases on harvest parameters. The aim of this paper is to investigate the PhenoCam approach within a field phenotyping setting, in which the performance of several genotypes is compared. First, we will introduce a data processing workflow that takes into account the special needs of crop monitoring in the context of field phenotyping. Then, we will follow the phenological signal throughout a season and evaluate the data for a soybean variety testing trial. We show how the phenological signal can be used to estimate in season traits such as early vigor, dynamic traits such as LAI development and senescence dynamics, which are very important traits in the context of plant breeding (Hund et al., 2019). Additionally, we demonstrate how information on phenological timing can be used to assess the impact of phenology on harvest traits such as thousand kernel weight (TKW), yield, and oil content.

METHODOLOGY
In the following, the dataset and methodology to extract in season and harvest traits from the phenological signal are introduced. Figure 1 gives an overview of the workflow.

Camera System and Dataset Description
We use images of a LUPUSNET HD -LE971 camera 1 , which provides a resolution of 1,920 × 1,080 pixels. It was mounted on a pole of the ETHZ field phenotyping platform  at 24.5 m height above the corner of our research field at the ETH Research Station for Plant Sciences Lindau-Eschikon, Switzerland (47.449 N,8.682 E; 556 m above sea level). Viewing direction was toward north, north-west. In 2015, we used 44,000 images recorded with a frequency of 4 to 60 images per hour. Figure 2 shows an example image of the soybean field from 2015.
The soil of the experimental field was slightly acidic (pH = 6.2) cambisol with sandy loam soil structure and low humus content. The soil analysis indicated sufficient availability for phosphate, potassium, and magnesium. 10 soybean cultivars (Glycine max [L.] Merrill) differing in the time needed for seed maturation and in nutritional food composition of the beans were selected ( Table 1). All of the cultivars were commercially available in Switzerland, and eight of them were on the national list of Swiss soybean varieties (Schwärzel and Hiltbrunner, 2015). 88 experimental plots of 1.5 m by 6.5 m containing seven rows of a particular soybean cultivar were sown on April 10, 2015. Each cultivar was planted in eight replications in a randomized complete block design, whereas four were used for ground validation measurements. During the season between the sowing date (10.04.2015) and the harvest (10.09.2015), crop cultivation measures were done according to best management practices for breeding trials. The genotype "Amphor" was excluded from analysis due to low viability resulting from old seed age (plots 1 https://www.lupus-electronics.de with bad growth in Figure 2A). The details on the set of genotypes are given in Table 1 (see Supplementary Table S2 for details of field management).

Extraction of Plot Values
To extract a representative value for each plot, we created a plot mask with a margin sufficient to keep the mask inside the plot's extent during the complete growing period. This is important, since the plot's position in the image slightly changes during the growing period due to the interaction of viewing angle and crop height. These masks are used to extract the RGB values of the plots. Depending on the position in the image, the masks of each plot contain between 1,000 and 9,400 pixels. The green chromatic coordinate (gcc; Gillespie et al., 1987) allows estimating the "greenness" of objects Sonnentag et al., 2012;Brown et al., 2016;Richardson, 2019). It divides the green pixel value by the sum of the pixel values of the green, red, and blue channels (Eq. 1). Previous studies showed that gcc compensates changing illumination conditions and can be used with non-calibrated cameras (e.g., Sonnentag et al., 2012). To exclude the influences of obstacles such as cables of the field phenotyping platform, we averaged all gcc values within the 25th to 75th percentiles of the pixels in each mask. We apply this methodology to all day images of the dataset. As a result, we get gcc values for each plot with a temporal resolution of the capturing frequency of the image. gcc = green/(green + red + blue) (1)

Extraction of Daily Signal
Different methods have been proposed to get a representative signal for each day. Often, the 90th percentile of all values per day, along with a 3-day moving average is used to minimize dayto-day variation mainly caused by changing weather conditions influencing the illumination Sonnentag et al., 2012). Other studies used a sigmoid function or logistic functions to normalize for these effects (Klosterman et al., 2014). After trying the median, mean, and different percentiles, we use the median per day to get a representative value per day, since we found that it is robust toward illumination changes (c.f. Figure 3). Nevertheless, the data is not fully free of day-to-day noise, which complicates the following analysis steps. Thus, we applied a Savitzky-Golay filter (Savitzky and Golay, 1964) with a 3rd order polynomial and frame length of 7, which is used widely to investigate time series data (Chen et al., 2004;Jönsson and Eklundh, 2004;Cao et al., 2018). The result is a daily phenological signal for each plot.

PhenoTimePoint and PhenoPhase Extraction
To extract the PhenoTimePoints, we scaled the uncorrected phenological signal for each plot from 0 to 1. Then, we extracted significant points from the phenological signal by identifying minima and maxima and slope-based features. At the beginning of the green up phase, we defined the start of green up point (SOG) as the point in time at which the difference between two consecutive days increased above 0.015, which corresponds to a slope of 0.015, since the values were scaled between 0 and 1. The first maximum of the signal after the green-up was defined as green peak (GP).
The last peak was defined as senescence peak (SP). The first minimum after the SP was defined as end of senescence/season (EOS). Table 2 summarizes the extracted features. For each plot, we extracted the days after sowing (DAS) for the identified PhenoTimePoint and corrected them for spatial effects (c.f. section "Spatial Correction"). Then, we calculated The cultivars "Falbala" and "Tiguan" were not recommended according to the Swiss cultivars catalog (Schwärzel and Hiltbrunner, 2015). 1 Reference cultivar Maple Arrow (0 days), 2 www.dsp-delley.ch, 05.03.2020, 3 Within the maturity group.
the duration of the PhenoPhases as the difference in DAS between two PhenoTimePoints. Figure 4 shows the concept for data of one plot.

Spatial Correction
Within breeding trials, it is common to correct for field heterogeneity resulting from differences in the soil or from slight differences in the management. Multiple approaches exist to correct for spatial effects (Gilmour et al., 1997;Piepho and Williams, 2010). Velazco et al. (2017) compared different mixed model approaches to correct for spatial trends. They found that the 'Spatial Analysis of field Trials with Splines' (SpATS) approach (Rodríguez-Álvarez et al., 2018), which uses twodimensional P-splines, performed comparably to more elaborate and trial-specific spatial models, with the advantage of being FIGURE 3 | (A) Evolution of the phenological signal of plot 4 denoted by the green chromatic coordinate (gcc) during the season. The colored dots represent the individual data points extracted from the images, while their color corresponds to the measurement time (around DAS 80 the measurement frequency of the camera was increased for a short period). The black dots and triangles are the daily averages derived from the median and the 90th percentile, respectively. The solid line represents the Savitzky-Golay smoothed median values. On the right axis and denoted with the gray ribbon is the corresponding heritability of the signal across the genotypes measured by means of heritability (H 2 ) across the season. (B) Images taken by the digital repeat photography device at six different days, denoted in days after sowing (DAS), during the growing season of soybean.
flexible and user-friendly. It calculates the field heterogeneity into the following model: where Y is the measured phenotypic value and f (r, c) is a smoothed bivariate surface defined over row (r) and range (c) positions across the field. The vector c g is the random coefficient of the genotypes associated with the design matrix Z g of the experiment, c r is the random coefficient of the rows associated with design matrix Z r , and ε is the random error vector (Rodríguez-Álvarez et al., 2018). Besides the field heterogeneity, our data is additionally influenced by the different viewing geometries within the image (c.f. Figure 2). Since the viewing geometry changes continuously across the image and consequently across the field, we consider it as a continuous component of spatial heterogeneity that adds to the spatial effects resulting from field heterogeneity and thus becomes part of the f (r, c) term. With SpATS, we calculated the Best Linear Unbiased Estimator (BLUE) for each genotype; Best, meaning they have the lowest variance, Linear, meaning they are linear functions of the data, and Unbiased, means the expected value of a mean estimate for an individual equals its true value (Molenaar et al., 2018). While best linear unbiased predictions (BLUPs) have some advances compared to BLUEs in many cases (Piepho et al., 2008;Molenaar et al., 2018), in our case the limited number of genotypes and the lack of information on relationships to a base population did not allow to accurately estimate the genetic variance components. Thus, we set genotype effects as fixed and used BLUEs to normalize the field heterogeneity and compare between genotypes.
Heritability measures the proportion of the overall observed phenotypic variance that can be explained by the genetic variance. Heritability of the spatially corrected traits was calculated according to (Rodríguez-Álvarez et al., 2018) based on the genetically effective dimensions provided by SpATS with genotypes set as random effect: where ED g is the effective dimension for the genotypes and m g is the total number of genotypes evaluated. The denominator (m g -1) reflects the upper bound for the effective dimension (see Rodríguez-Álvarez et al., 2018 for further details). We calculated genotypic BLUEs and heritability for the phenotypic gcc values of every measurement date, the PhenoTimePoints before calculating the PhenoPhases and the field measured crop traits. Throughout the manuscript, the spatial corrected genotypic BLUEs are referred to as genotypic values, whereas the observed values are referred to as phenotypic values.

In-Season Trait Estimation
To evaluate and validate the feasibility and relevance of the extracted PhenoTimePoints and PhenoPhases in relation to crop growth and development in a breeding trial context, we selected early vigor, leaf area index (LAI) and senescence. All these traits are known to be an important selection criteria for soybean breeding in Switzerland and similar environments.
Early vigor was rated as an indicator for plant health and growth development at 40 DAS. The plots were rated by two persons, independently of the germination rate, and the values were averaged. The used scale ranged from one to eight. One means that the plant is small and depicts a low vitality, whereas eight stays for very vigorous plants showing a large and vital plant habitus. The rating mainly integrates size and color information (Peter et al., 2009). The LAI was measured in regular intervals during homogeneous weather conditions such as cloud-free sky or constant cloud cover with an LAI-2200 Plant Canopy Analyzer (LI-COR, Inc., Lincoln, NE, United States). The LAI measurement combined a first measure above the canopy, followed by ten representative measurements diagonally spread through the plot canopy followed by a second above canopy measurement. With a 270 • view-restricting cap, the potential negative influence by the measuring person was reduced. The senescence was estimated as the percentage of yellow and brownish colored leaves as well as fallen leaves according to (Munger et al., 1997). The senescence rating was done six times for each plot at 109, 117, 124, 132, 137, and 150 DAS, respectively. For each trait, we calculated the genotypic BLUEs and correlated both the (spatially not corrected) phenotypic and genotypic BLUEs with the uncorrected and spatially corrected phenological signal.

Harvest Trait Estimation
To evaluate the relevance of the extracted PhenoTimePoints and PhenoPhases for overall crop performance, we investigated their relationship to the harvest traits yield, thousand-kernel weight (TKW), and seed oil content. Yield in t ha −1 was estimated using the seeds of all hand-harvested pods of 1 m of plants dried for 3 days at 40 • C. Due to the experimental setup, machine harvest was only possible at the time of maturity of the latest ripening genotypes, when early genotypes were already affected by pod and kernel fall. The TKW (g) was calculated by counting and weighing the seeds of 10 randomly selected pods. Oil content (percentage of dry kernel weight) was measured in dry grains with a near-infrared spectrometer (NIRS, Infratec 1241 grain analyzer, Foss GmbH, Rellingen, Germany) using the company provided soybean program. Also, for the harvest traits, we calculated the genotypic BLUEs and correlated these with the BLUEs of the PhenoTimePoints and the duration of the PhenoPhases. Figure 3 shows the evolution of the gcc-based phenological signal of plot 4. After sowing, the phenological signal did not change much until around DAS 50. During this time, the vegetation growth was very limited due to the low temperatures. The fluctuations were mainly due to changes in soil moisture due to rain, although at DAS 23, the growth also started. At DAS 48, the gcc steeply started to increase. This increase was followed by a shoulder of high greenness around DAS 76, after which the greenness decreased before it rose again to a second peak around 125 DAS. After this second peak, the greenness steeply decreased until around DAS 144.

Phenological Signal
The black triangles and dots in Figure 3 represent two different approaches to extract a daily signal, namely the 90th percentile often found in literature and the median, respectively (c.f. subsection "Extraction of Daily Signal"). In direct comparison, the median values appeared less noisy than the 90th percentile, in particular when the environmental conditions differed, such as between DAS 50 and 60. Thus, we continued with the daily median to generate the phenological signal. Still, the median was also affected by changing illumination conditions, but the Savitzky-Golay filter was able to reduce these day-to-day fluctuations. The heritability of the phenological signal varied. Between DAS 15 and DAS 20, it increased to 0.95. At DAS 55, it slightly decreased to 0.925 at DAS 60, 0.85 at DAS 70, and 0.33 at DAS 80. Toward DAS 90, it increased again over 90 before it decreased again after DAS 140. Figure 5 shows the phenotypic and genotypic values of the different plots of the phenological signal at different DAS and for different PhenoTimePoints in relation to the viewing geometry. The influence of the viewing geometry on the phenological signal changed throughout the season. During the early development (DAS 40 and 50), the signal was positively related to viewing geometry. During the green up (DAS 61), only a weak trend was found. Around the GP (DAS 80) and SP (DAS 120), the relationship was negative. At the end of the season, only a weak positive trend was observed. Overall, the viewing geometry had only a weak relationship with the signal (R 2 < 0.08; Supplementary Table S3). For DAS 80 and 120, the viewing geometry had an R 2 of 0.55 and 0.26 with the gcc (Supplementary Table S3). In all cases, the spatial correction could mitigate the viewing geometry effect. In all cases, the timing of the PhenoTimePoints had a weak positive trend with the viewing geometry (R 2 < 0.09; Supplementary Table S3). Also here, the spatial correction could reduce this trend. The exact quantitative description of the relationships can be found in Supplementary Table S3.

Spatial Correction
The gray bars in Figure 6 show the heritability of the PhenoTimePoints and PhenoPhases. After spatial correction, all PhenoTimePoints showed a high heritability (H 2 ) above 0.76. Besides of the PhenoPhase EOS_GP (H 2 of 0.46), all PhenoPhases also showed high to very high heritability (0.68-0.97). For exact timings with their variability and heritability, see Supplementary Table S1.

In Season Trait Estimation of Leaf Area Index, Early Vigor, and Senescence
The gcc represents an integrated signal of canopy greenness. Until 70 DAS, the LAI increased to 6. Toward higher values, the scatterplots (Figure 7) show that the gcc signal saturated at an LAI of about 4.5. At low LAI (<1), the scatterplots indicate that there was no clear relationship between the LAI meter and gcc estimated LAI values. The LAI values of DAS 61 ranged from 0.35 to 2.93. In this span, the gcc corresponded well with the LAI on the phenotypic (R 2 of 0.8) and genotypic level (R 2 of 0.84). At DAS 40, a good relationship was found for plant vigor on the phenotypic level (R 2 of 0.64) and the genotypic level (R 2 of 0.78).
Visual senescence rating and gcc were not linearly related over time (Figure 8). While the visual rating constantly decreased over time, gcc showed a clear bimodal trend with a steady increase followed by a sharp decrease. Until a senescence rating of about 65%, there was a positive linear trend of gcc with the rating. After 65%, there was a strong negative trend (R 2 of 0.78 and 0.69 on the phenotypic and genotypic level, respectively) of gcc and the rating. This also showed that the SP approximately corresponds to a senescence rating of 65%. Figure 4 shows the relation of the phenological signal to the manual rating of senescence and LAI of one plot: The LAI rapidly increased after SOG toward the GP. Around the GP, the LAI increase slowed. The senescence steeply increased toward the SP and reached 100% at EOS. The SP approximately corresponded to a senescence rating of 65%. Figure 6 shows boxplots of the timings of the PhenoPhases and PhenoTimePoints. The timings of the PhenoTimePoint SOG only showed a small variability of 3 days. The GP, SP, and EOS showed greater variability (more than 7 days). The PhenoPhase SP_SOG showed a high range of durations (more than 14 days) being in the same order of magnitude as the maturation difference known for the selected genotypes (between 1 and 3 weeks). SP_SOG showed a high range of 14 days, while the difference in duration of the other PhenoPhases was between 6 and 7 days.

Interaction Between Phenological Timing and Harvest Traits
With the information on the timing and duration of the PhenoTimePoints and PhenoPhases, it is possible to correlate these with harvest traits. In the following, we show the potential of the approach by discussing (a) how the early vigor corresponded to the timing of the PhenoTimePoints of the different genotypes and (b) how the timing of PhenoPhases corresponded to the harvest traits dry thousand kernel weight, oil content and pod yield.  Figure 9A shows a scatterplot of the vigor rating at DAS 40 to the timing of the GP of the different genotypes. A high rating of vigor strongly corresponded (R 2 of 0.87) to an earlier GP. For the other PhenoTimePoints the relationship was moderate ( Figure 9B). All relationships were negative, which indicated that a high early vigor generally led to earlier development. The PhenoPhases had a moderate to low (R 2 < 0.4) relationship to the early vigor rating and were negative. Only the duration between SP and EOS had a moderate (R 2 of 0.36) positive correlation to early vigor, which indicated that a better early growth prolongs senescence. Early vigor had no significant (p of 0.01) correlation to the total duration of growth after the GP peak was reached. Figure 10 shows the correlation of the harvest traits yield, TKW, and oil content with the timings of the PhenoTimePoints and PhenoPhases. It revealed that varieties that had an earlier start of the green up phase and reached the GP earlier had a higher TKW. Varieties that reached the GP later had a higher oil content. Varieties that had a shorter green period (SP_GP) and an earlier SP had a higher yield, while varieties with a longer senescence period (EOS_SP) had a higher yield. Similarly, the oil content was higher when the duration between the GP and SP or EOS was longer. It has to be noted that these results might be specific for this experiment reflecting the genotype selection.

Generation of the Phenological Signal
PhenoCam data are generally produced by uncalibrated cameras and are influenced by many factors. In this study, we identified (i) changing measurement geometry during the day, other (ii) changes in illumination conditions due to varying cloud cover, and (iii) the viewing geometry within an image as the main influences on the phenological signal. The gcc is known to suppress the effects of changes in scene illumination . During initial tests, we could confirm that the gcc was more stable than other color indices (data not shown). To account for short term fluctuations of the signal (i, ii) we calculated the daily median to normalize for the diurnal differences and used the Savitzky-Golay filter to buffer day-to-day variations resulting from short term differences due to changing weather conditions. The approach showed to be more successful than other approaches (e.g., 90th percentile) to derive a smooth signal across the season. Still, two uncertainties remain: With low crop cover, the soil background affects the signal. Wet periods of several days reduce the gcc value as seen around DAS 25. Also, plant canopies appear differently in direct and diffuse illumination due to the varied shading within the canopy (Aasen and Bolten, 2018), which is still an open problem in remote sensing and field phenotyping (e.g., Yu et al., 2017) and hard to address in a general manner without actual information on canopy structure. Nevertheless, the daily mean in combination with the Savitzky-Golay filter seemed to normalize fluctuating illumination conditions quite well. The different pixel count per plot, which results from different distances to the camera, can be assumed to have negligible effect on the signal, since pixels covering a larger area represent the average signal of this area. Still, a sufficient resolution is required for cameras to ensure that pixels potentially contaminated from border effects can be excluded.
The systematic bias resulting from viewing geometry (iii) could be mitigated by the spatial correction approach along with potential in-field heterogeneity. The gcc at most DAS and most timings of the PhenoTimePoints were influenced by the oblique viewing geometry (Figure 5), as has also been reported for LAI in forest environments (Keenan et al., 2014). The main reason for this is that the zenith angle interacts with canopy structure, which results in occlusion effects (Aasen and Bolten, 2018;Roth et al., 2018). Considering the pattern caused by the viewing geometry ( Figure 2B) as observational bias during the spatial correction could remove the trend ( Figure 5) and most of the remaining variance could be attributed to the differences between the genotypes as demonstrated by the high heritability (Figure 6). The drawback of this approach is that it relies on an experimental design that is suited for spatial correction procedures, such as those implemented in SpATS or similar approaches. Still, the viewing geometries mainly influenced the absolute values of the gcc. In comparison, the timing of the PhenoTimePoints was barely influenced (Figure 5 and Supplementary Table S3). Thus, two cases can be distinguished: (a) when the gcc should be used as a proxy for a trait (e.g., LAI and vigor), a spatial correction seems to be necessary; and (b) when the gcc dynamics are used to estimate PhenoTimePoints it does not seem to be absolutely necessary to apply a spatial correction. Regarding the absolute values, it also has to be noted that these results will mainly depend on the differences between the viewing geometries of each plot. When a camera is mounted further away from the experiment, the viewing geometries for each plot will be more similar and, consequently, will have less impact on the results. Additionally, the absolute influence of the viewing geometry will always depend on the canopy structure and thus differ between crops. Besides, further studies could advance the spatial correction approach by separating the effects resulting from the viewing geometry from fieldheterogeneity effects.

Extraction of Traits, PhenoTimePoints, and PhenoPhases
The analysis showed that the gcc could be used to estimate LAI between approximately 1 and 4.5. Higher LAI values were not detectable due to saturation of the gcc resulting from the mentioned occlusion effects in combination with almost full canopy cover. The inaccuracies of lower LAI values likely result from the inaccuracies from our reference method. It is known that the Plant Canopy Analyzer (LAI-2200) cannot capture low LAI (<1) values reliably in row crops (c.f. Roth et al., 2018). Since crop growth is most rapid for LAI 1 to 4.5 (e.g., Figure 4), capturing the development in this phase is already a benefit for many applications in breeding and variety testing. It is likely that FIGURE 8 | Relationship of senescence ratio to gcc on the phenotypic (A) and genotypic level (B) in days after sowing (DAS). The color corresponds to the senescence rating. The coefficient of determination (R 2 ) is given for the relationship of gcc and senescence for a rating above 65%. the exact relationships between gcc and LAI are crop and site specific (e.g., in brighter soils also lower LAI values will be visible) and will be influenced by the sowing density and row orientation. Nevertheless, the general concept is broadly applicable. Future studies should investigate the relationships for other crops and sowing patterns.
We could identify four points within the phenological signal that could be attributed to phenological events, namely the start of green up (SOG), green peak (GP), senescence peak (SP), and end of senescence/end of season (EOS). Different terminologies have been used in the ecological literature, and we tried to align with the existing terminology. In the case of SOG, we refrain from calling it start of season, since the plants did already grow before, but our model could not reliably resolve it (see below). The GP and SP are uncommon in the ecological literature. But, the analysis showed that these PhenoTimePoints can be attributed to quantitative ratings for LAI and senescence of soybean. Figure 7 shows that the maximum value of the gcc is around an LAI of 4.5, and consequently, the peak during green up (GP) can be attributed to this value. Figure 8 shows that the gcc reaches a maximum during senescence at a rating of 65%, and consequently, the peak of gcc during the senescence phase (SP) can be attributed to this value. As previously mentioned, it has to be noted that it is likely that these values are species specific because they are influenced by canopy structure. Also, it has to be mentioned that senescence in field crops often starts at leaf levels close to the ground. This is a shortcoming of the presented method since oblique images only capture the top of canopy. This also needs to be taken into account when setting up a camera, since in crops where the canopy is not 100% closed at the time of senescence-such as cereals-viewing geometries close to nadir will allow seeing into the canopy and thus, possibly capture processes earlier than in parts of an image with oblique viewing geometries. Thus, ideally, a camera should be placed in some distance to the experiment such that the viewing geometry differences are small. Finally, the method to extract the PhenoTimePoints relies on significant changes in the slope of the gcc. This turned out to be a shortcoming in extracting the point for start of senescence, which did not work reliably for all genotypes. We ascribe that to the fact that some plots already started a very slow yellowing process around DAS 99, which did not produce a pronounced slope change (Supplementary Figure S1).

Investigating Interactions of Growth, Harvest Traits and Phenology With PhenoCams
While there is a conceptual understanding on how phenology may be influenced by climate conditions, it is challenging to provide quantitative estimates of the magnitude of these shifts (Richardson et al., 2013;Herrera et al., 2018). Field phenotyping with PhenoCams is a low cost, fully automated approach to gather very high temporal resolution data of an area sufficient for experimental trials. We could show that such an approach is particularly useful to monitor (i) rapid growth dynamics that might only differ by or last a few days, (ii) investigate differences in phenology across genotypes, and (iii) analyze the interactions between crop development (e.g., early vigor, LAI, and green up dynamics), phenology (e.g., senescence dynamics) and harvest traits.
With the highly temporal-resolved data, we could investigate the influence of early vigor on subsequent PhenoTimePoints and duration of PhenoPhases. While the focus of this manuscript is on the methodology, the obtained results are plausible. Early vigor is positively correlated with GP showing the faster canopy closure of more vigorous genotypes. Similar observations were reported using canopy cover measures instead of GP (Kipp et al., 2014;Kirchgessner et al., 2017). Besides, we found that yield and yield-quality traits such as oil content and TKW are related to a longer green period, late and short senescence, which has previously been reported to correspond to higher yield performance under normal or stress situations (Kropff et al., 1993;Thomas and Smart, 1993;De Souza et al., 1997;Rebetzke et al., 2016). While these interactions are not new knowledge per se, so far it was very hard to quantify these relationships for multiple genotypes in the field. Because these insights are connected to differences between genotypes, analyzing datasets of multiple years would also allow assessing, e.g., the impact of climate-induced phenological shifts.

Comparison to Other Approaches
PhenoCam data allows monitoring multiple areas of interest at the same time (e.g., Julitta et al., 2014;Browning et al., 2017; this study) in very high temporal resolution. As shown, they can monitor various types of traits across the whole season with minimal effort. Still, the approach is limited to some basic traits that are observable with the spatial ground sampling distance of the cameras. Besides, the approach is not suitable to estimate biochemical traits due to the missing radiometric calibration and spectral characterization. While encouraging approaches exist, e.g., using very high-resolution gray card calibrated nadir images along with plant segmentation to estimate leaf nitrogen concentration levels (Wang et al., 2014), the bandwidth of usual RGB cameras are typically too wide to fit the absorption features of e.g., pigments. Continuous spectral measurement systems may be superior for that purpose (e.g., D'Odorico et al., 2015;Paul-Limoges et al., 2018;Wieneke et al., 2018). However, temporally continuous spectral measurements are mostly limited to point measurements (Aasen et al., 2019).
For a field phenotyping scenario compared to the one of this study, where a couple of hundred plots (depending on the plot size) on an area of approximately 2 ha are to be screened for basic phenological parameters, the PhenoCam approach has lower operational and setup cost, higher temporal resolution but lower precision compared to (more or less) mobile platforms that monitor the canopy from close range (e.g., Cendrero-Mateo et al., 2017;Kirchgessner et al., 2017;Hund et al., 2019). Flying mobile platforms (UAV) that monitor the canopy from a distance of a couple of decameters with RGB sensors are quite flexible, can be used to monitor phenology  and can even be used to track the green up based on plant segmentation (e.g., Roth et al., 2018). But while RGB flights are almost always possible with a high temporal resolution (Aasen and Bareth, 2018;Roth et al., 2018), each flight campaign is connected with some effort, because fully autonomous systems are currently prohibited by legislation in most countries. Generally, mobile platforms can cover a larger area. Close-range PhenoCams that are mounted right above the canopy could combine the benefits of both approaches but would be very costly to set up. Table 3 summarizes these considerations. Looking forward, one could imagine a network of multiple PhenoCams that cover larger trials, even across different environments and for multiple years. In the end, many traits need to be captured by field phenotyping approaches to complement each other in order to get the full picture. PhenoCams can become an important tool providing high temporal resolution information on phenological processes.

DATA AVAILABILITY STATEMENT
The datasets generated for this study are available on request to the corresponding author.

AUTHOR CONTRIBUTIONS
HA developed the idea, conceived the study together with FL, analyzed the data, and wrote large parts of the manuscript. NK operated the PhenoCam as part of the FIP system, provided the technical details, and supported the data analysis. FL organized and conducted the field experiment and supported the data analysis. HA, NK, AW, and FL contributed to writing of the manuscript.

FUNDING
The project was supported by Agroscope within the project: «Knowledge transfer for non-destructive measurement techniques for agricultural experiments» under the contract ID: 655017678.

ACKNOWLEDGMENTS
We thank Hansueli Zellweger for the preparation and management of the field, Brigitta Herzog for the preparation of the seeds and her support with the sample processing, Lukas Müller for doing most of the field measurements in the frame of his master thesis at ETH, Marius Hodel for generating the shapes of the plots, and Andreas Hund for helping with the spatial correction. We also thank the four reviewers for their constructive and detailed comments that helped to improve the quality of the manuscript.