Modeling of Diurnal Changing Patterns in Airborne Crop Remote Sensing Images

: Airborne remote sensing technologies have been widely applied in ﬁeld crop phenotyping. However, the quality of current remote sensing data suffers from signiﬁcant diurnal variances. The severity of the diurnal issue has been reported in various plant phenotyping studies over the last four decades, but there are limited studies on the modeling of the diurnal changing patterns that allow people to precisely predict the level of diurnal impacts. In order to comprehensively investigate the diurnal variability, it is necessary to collect time series ﬁeld images with very high sampling frequencies, which has been difﬁcult. In 2019, Purdue agricultural (Ag) engineers deployed their ﬁrst ﬁeld visible to near infrared (VNIR) hyperspectral gantry platform, which is capable of repetitively imaging the same ﬁeld plots every 2.5 min. A total of 8631 hyperspectral images of the same ﬁeld were collected for two genotypes of corn plants from the vegetative stage V4 to the reproductive stage R1 in the 2019 growing season. The analysis of these images showed that although the diurnal variability is very signiﬁcant for almost all the image-derived phenotyping features, the diurnal changes follow stable patterns. This makes it possible to predict the imaging drifts by modeling the changing patterns. This paper reports detailed diurnal changing patterns for several selected plant phenotyping features such as Normalized Difference Vegetation Index (NDVI), Relative Water Content (RWC), and single spectrum bands. For example, NDVI showed a repeatable V-shaped diurnal pattern, which linearly drops by 0.012 per hour before the highest sun angle and increases thereafter by 0.010 per hour. The different diurnal changing patterns in different nitrogen stress treatments, genotypes and leaf stages were also compared and discussed. With the modeling results of this work, Ag remote sensing users will be able to more precisely estimate the deviation/change of crop feature predictions caused by the speciﬁc imaging time of the day. This will help people to more conﬁdently decide on the acceptable imaging time window during a day. It can also be used to calibrate/compensate the remote sensing result against the time effect.


Introduction
Modern plant breeding efforts depend on phenotypic data to select high-yield and stress-tolerant plants quickly and efficiently [1,2]. Ag remote sensing technologies have been developing rapidly for many years. Currently, field phenotyping activities are being conducted with satellites, airborne platforms (manned and unmanned), and ground-based vehicles [3,4]. Various sensors such as RGB (Red-Green-Blue), hyperspectral and thermal cameras are carried by these platforms to take images of the crop field. These technologies have been proven effective through various Ag remote sensing projects [3,5]. However, the quality of Ag remote sensing data is still limited by various noise sources such as the changes of day light, wind speed, temperature, sun angle, etc. [6][7][8][9]. Among these noises, the diurnal variability is one of the major factors causing significant quality issues in Ag remote sensing.
The diurnal impact to plants' reflectance characteristics is a complicated process. It introduces strong noises in plant phenotyping result [8,10,11]. For example, the reflectance characteristics of the same plants at noon can be very different from those in the afternoon. These variabilities are caused by interactions between camera's sensitivity, camera's view angle, canopy geometry, solar zenith angle, solar azimuth angle, and shadows [12,13]. Meanwhile, the plant itself is endogenously sensitive to the environmental conditions with the complicated interactions between the genetic backgrounds, the external environments and the treatments (G (Genotype) * E (Environment) * t (Treatment)). All these affect the final reflectance characteristics of plants throughout the day, causing the diurnal variabilities.
The diurnal variance on phenotyping data has been reported as unignorable impacts in many plant studies. Gardener [14] stated that it was a major unresolved noise issue in using reflectance measurements for estimating leaf area, plant biomass, or phenology as it was affected by diurnal changes. In the last decades, diurnal variability has been documented in various plant phenotyping studies for corn [15], soybean, wheat [11] with both passive [7,16] and active sensors [15]. These variances are retained in the captured images, weakening the signal power of the data. Sometimes, the generated variances are even larger than the plant differences caused by biotic or abiotic stresses, which severely limits the accuracy in phenotyping data [11,16]. However, for most of current remote sensing studies, people rarely consider the impacts from diurnal variabilities, which introduces severe noises into the final analysis results. For example, the NDVI values had over 10% difference at different time points on the same plant from the raw remote sensing measurements [11,16]. Similar diurnal variances are also found in many other plant features such as plant temperature, spectra features, chlorophyll content, etc. [9,17]. Therefore, it is extremely important to reduce the diurnal variabilities in the data for improved remote sensing quality.
With the development of plant phenotyping, researchers are better aware of the diurnal variabilities while there are limited studies on how to deal with diurnal variability in Ag remote sensing. To reduce impacts of diurnal variability, remote sensing technologies, such as unmanned aerial vehicles (UAV), usually have strict rules on the sampling time window and weather condition requirement [18][19][20]. According to Bellvert et al. [21], the proper time of the day to acquire thermal and multispectral images should be around noon, because there is almost complete absence of shadow effects. Meanwhile, plant physiology changes on a cyclical diurnal nature based on photosynthetic activity and processes dependent on incident solar radiation [22]. Consequently, UAV data collections are required to operate during a certain range of time for accurate monitoring crop physiology [20].
It is important to quantitatively model the diurnal variability of crop phenotyping features for improved Ag remote sensing quality. The combined PROSPECT leaf optical properties model and SAIL canopy bidirectional reflectance model (PROSAIL) has been widely used to predict the change of plant canopy spectral reflectance influenced by the changing environmental conditions such as solar angle [23,24]. However, the model does not meet the accuracy requirement in plant phenotyping remote sensing. For example, Berger et al. [25] compared PROSAIL's simulation result with the spectra data collected from the field and found that PROSAIL's predicted spectra showed severe drifts from the field measurements especially at the green and red edge (~700 nm) wavelengths [25]. This was also confirmed in the field remote sensing experiment in 2019 at Purdue University, in which the PROSAIL model only predicted less than 5% of the diurnal variance of the observed NDVI from a hyperspectral camera. Another potential solution arises from the use of different regression methods to predict and compensate the diurnal effects. For example, a correction model with the polynomial regression method was developed to predict the crop reflectance as a function of solar zenith angle, time of day, and ICI (instantaneous clearness index). The capability of the model in reducing the diurnal variation with Green Normalized Difference Vegetation Index (GNDVI) and some individual bands [7] Remote Sens. 2021, 13, 1719 3 of 19 was tested. However, the plant data was collected on a small portion of the leaf by a handheld radiometer with four bands, which may not properly simulate airborne remote sensing platforms carrying hyperspectral or multispectral cameras. The simple polynomial regression model could successfully describe the changes in data over three consecutive days. However, it may fail to represent the general pattern on other days when the plants are at different stages of their growth cycle. Therefore, a more accurate diurnal variability model is still urgently needed.
Modeling diurnal variability in plant phenotyping data is difficult because of the complicated interactions between plants and real-time environment conditions such as cloud coverage, wind and so on. It is necessary to collect hundreds or thousands of time series images of the same field in order to model the diurnal changes, but this had been difficult with the current airborne platforms. For instance, the UAV platforms were rarely reported to image the same field for more than one time a day due to various environmental and logistical limitations [26]. In order to solve this problem, the Ag engineers at Purdue University deployed a fixed field VNIR hyperspectral gantry platform as a "mock drone" system in Purdue's research farm in 2019, which is similar to the phenotyping system called Terra-Ref [27]. With this field gantry imaging system, hyperspectral images of the same field can be continuously collected under various weather conditions throughout the daytime in minutely base. A local weather station and distributed soil sensors in the field provide synchronized environmental condition data with the images.
Recent advances in remote sensing technologies enable the measurement of various image-based phenotypic features such as vegetation indices (e.g., NDVI), biochemical constituents (e.g., chlorophylls), water, or raw spectral features, that dramatically boosts modern crop studies [3,28]. However, as described above, the problem of diurnal variations in these phenotypic features is still unsolved due to the limits from both hardware and software aspects. In order to close this gap, in this paper, we performed a comprehensive investigation on the diurnal variability in over 8000 repeated hyperspectral images of the same corn field taken by the field imaging gantry at Purdue University over the 2019 growing season. Images of corn plants from vegetative stage V4 to reproductive stage R1 were captured. Hyperspectral image processing algorithms were applied to calculate the reflective canopy spectra and predict plant physiological features such as NDVI, RWC, and so on. The time series decomposition method was applied to obtain diurnal curves for these plant phenotyping features. After combining the imaging results over 31 days, it clearly showed repeatable diurnal changing patterns for NDVI and other phenotyping features. Regression models were successfully developed to precisely describe these diurnal changing patterns. With the result of this work, Ag remote sensing users will be able to more precisely understand the deviation/change of crop feature predictions caused by the specific imaging time of the day. It'll also be easier to decide on the acceptable imaging time range and calibrate the diurnal factors with the developed models.

High-Throughput Field Imaging Acquisition System
The field VNIR hyperspectral gantry platform at Purdue University's Agronomy Center for Research and Education (ACRE) was used to collect imaging data in this study. A weatherproof VNIR push-broom hyperspectral camera (MSV-101-W, Middleton Spectral Vision, Middleton, WI, USA, Table 1) was carried by the 7 m high gantry platform to scan a 50 m by 5 m strip field under a wide range of weather conditions. The VNIR images had spectral range from 376 to 1044 nm with spectral resolution of 1.22 nm, and spatial resolution of 0.5 cm/pixel ground sample distance (GSD). The system could be configured to automatically scan the crops in the field repeatedly. It takes 6.5 min to scan the 250 square meters field, and the scanning frequency can be higher if a sub-portion of the field needed to be scanned ( Figure 1). This single sided imaging gantry stands by the north side of the field and the length of the camera structure on the top was restricted by a certain ratio of the gantry's height. This unique design kept 100% of the shadow away from the crops. This shadow-free feature enabled the system to more realistically simulate drone remote sensing in the field. The hyperspectral camera utilizes sunlight as the lighting source, so the gantry system can function any time after sunrise until sunset on each day. A white reference panel was installed 0.5 m underneath the hyperspectral camera and moves together with it. A local mini weather station and distributed soil sensors were installed in the same field to collect real-time environmental condition data such as temperature, solar irradiation, wind speed and soil moisture when each image was taken.

Experiment Design and Data Collection
To study the diurnal variation in phenotyping data, two genotypes of corn plants, including genotype B73 × Mo17 and P1105AM, were grown in the field underneath the gantry in 2019. Each genotype was treated with three different nitrogen solutions: high nitrogen (HN) with 56 kg/ha (32 mL 28-0-0 in 1 L water), medium nitrogen (MN) with 28 kg/ha (16 mL 28-0-0 in 1 L water), and low nitrogen (LN) with 0 kg/ha (water). Each of the genotype by nitrogen treatment combination is repeated in 5 2-rows-by-3-m mini-plot replicates so there are 30 plots in total in the field ( Figure 2).
In order to capture the instant impacts to the images from the changes of cloud coverage, wind speed and other environmental conditions, the team decided to have an imaging frequency of every 2.5 min. This allowed us to scan up to six different plots considering the extra time needed for data transfer, real-time image processing, and homing the gantry cart. The six plots were selected to cover all three nitrogen levels and two genotypes.
The continuous imaging started when the corn plants reached the V4 stage and lasted until 31 days later when the plants are at R1 stage on average ( Figure 3). On each of the 31 days, imaging started at 8:00 a.m. and ended at 7:30 p.m. Around 280 hyperspectral images were collected with the repetitive imaging frequency of every 2.5 min. The gantry was only turned off during extreme weather conditions such as thunderstorms to protect the equipment's safety. By the end of the 31 days, a total of 8631 hyperspectral images were collected.  (h) binary image after segmentation; (i) layout of the plots with three nitrogen treatments and two genotypes. The green boxes were for Genotype P1105AM, and the blue boxes were for Genotype B73 × Mo17. Nitrogen treatments of HN, MN and LN were also labeled; (j) NDVI heatmap; (k) the spectra from different genotypes and nitrogen treatments.

types.
The continuous imaging started when the corn plants reached the V4 stage and lasted until 31 days later when the plants are at R1 stage on average ( Figure 3). On each of the 31 days, imaging started at 8:00 a.m. and ended at 7:30 p.m. Around 280 hyperspectral images were collected with the repetitive imaging frequency of every 2.5 min. The gantry was only turned off during extreme weather conditions such as thunderstorms to protect the equipment's safety. By the end of the 31 days, a total of 8631 hyperspectral images were collected. In the middle of the project when the plants were at V9, we collected ground truth measurements such as nitrogen content, RWC and plant fresh weight. Two plants were randomly sampled from each plot. The plant shoot was cut to measure the fresh weight. A small section (2.5 × 5.0 cm) of the top collared leaf was taken to measure the RWC using the Equation (1) [29,30]. The remaining part of top collared leaf was sent to the great lakes A and L laboratories (A and L Great Lakes Laboratories, Inc., Fort Wayne, IN, USA) for measuring the nitrogen percentage.
where FW is fresh weight, DW is dry weight, and TW is the turgid weight.

Image Segmentation and Feature Extraction
After data collection, standard imaging processing protocols were performed to extract the interested plant phenotyping features [5,31,32]. The raw hyperspectral images were firstly calibrated with the real-time white reference. It ensured each scanning line from this push-broom sensor was calibrated with the reference under the same lighting conditions. The image calibration was performed with the following equation: (2) where is the calibrated image, is the raw hyperspectral image, is dark reference image and is the hyperspectral image of the white reference. The calibrated images were then processed using a segmentation procedure with convolution methodology [33,34]. A vector of sequential integers from −20 to 20 was multiplied by the reflectance intensity vector from the red-edge region (680-720 nm). By choosing threshold 7 as the boundary, the plant tissue was successfully segmented from the background (Figure 1h). In some of images, there were some weeds which were ignorable in size compared to the corn plants (See the red boxed in Figure 1h).
The average reflectance spectrum from each plot was calculated. In total, 51,786 spectra (8631 images * 6 plots/image) were calculated for the plots with different genotypes and nitrogen treatments. These spectra were used to calculate the crop remote sensing results such as NDVI and predicted RWC. The formula below was used for calculating In the middle of the project when the plants were at V9, we collected ground truth measurements such as nitrogen content, RWC and plant fresh weight. Two plants were randomly sampled from each plot. The plant shoot was cut to measure the fresh weight. A small section (2.5 × 5.0 cm 2 ) of the top collared leaf was taken to measure the RWC using the Equation (1) [29,30]. The remaining part of top collared leaf was sent to the great lakes A and L laboratories (A and L Great Lakes Laboratories, Inc., Fort Wayne, IN, USA) for measuring the nitrogen percentage.
where FW is fresh weight, DW is dry weight, and TW is the turgid weight.

Image Segmentation and Feature Extraction
After data collection, standard imaging processing protocols were performed to extract the interested plant phenotyping features [5,31,32]. The raw hyperspectral images were firstly calibrated with the real-time white reference. It ensured each scanning line from this push-broom sensor was calibrated with the reference under the same lighting conditions. The image calibration was performed with the following equation: where R cali is the calibrated image, R raw is the raw hyperspectral image, R dark is dark reference image and R white is the hyperspectral image of the white reference. The calibrated images were then processed using a segmentation procedure with convolution methodology [33,34]. A vector of sequential integers from −20 to 20 was multiplied by the reflectance intensity vector from the red-edge region (680-720 nm). By choosing threshold 7 as the boundary, the plant tissue was successfully segmented from the background (Figure 1h). In some of images, there were some weeds which were ignorable in size compared to the corn plants (See the red boxed in Figure 1h). The average reflectance spectrum from each plot was calculated. In total, 51,786 spectra (8631 images * 6 plots/image) were calculated for the plots with different genotypes and nitrogen treatments. These spectra were used to calculate the crop remote sensing results such as NDVI and predicted RWC. The formula below was used for calculating NDVI.
where R 800nm and R 650nm are the reflectance values of wavelength 800 nm and 650 nm respectively [35,36]. The partial least squares regression (PLSR) model was used to predict RWC from the spectra. In order to avoid prediction drifts between facilities [37][38][39], instead of using an existing RWC prediction model developed from the other facilities, the team decided to build a new PLSR model with the spectra and RWC ground truth data collected in the same project. The new model predicted RWC with the cross-validation coefficient of determination (R 2 ) of 0.722 and root mean square error (RMSE) of 6.22% ( Figure 4). This RWC model was then applied to all the 51,786 spectra to predict RWC in those plots over the 31 diurnal cycles.

Data Quality Check
After calculating the remote sensing result data such as spectra, NDVI and RWC, the data quality was checked, and outlier data was removed. For each feature from each day, the measurements between the upper inner fence (Q3 + 1.5IQR) and lower inner fence (Q1 − 1.5IQR) were kept [40]. IQR is the interquartile range, being equal to the difference between 75th (Q3) and 25th (Q1) percentiles. This quality filtering removed the blunders/gross errors in the image acquisition setup/process. For example, some of the removed outlier images were taken under very-high wind-speed conditions, when the plant

Data Quality Check
After calculating the remote sensing result data such as spectra, NDVI and RWC, the data quality was checked, and outlier data was removed. For each feature from each day, the measurements between the upper inner fence (Q3 + 1.5IQR) and lower inner fence (Q1 − 1.5IQR) were kept [40]. IQR is the interquartile range, being equal to the difference between 75th (Q3) and 25th (Q1) percentiles. This quality filtering removed the blunders/gross errors in the image acquisition setup/process. For example, some of the removed outlier images were taken under very-high wind-speed conditions, when the plant stems were bent severely, making the images look significantly different. The data before 10:00 a.m. and after 5:30 p.m. also showed severe variances and noises. This could be caused by the dews on the surface of leaves and dim lighting conditions [41]. Only the data collected during 10 a.m.-5:30 p.m. were used in the diurnal pattern analysis as very few airborne remote sensing activities happen outside this time range. We ended up using this selected data pool (Table 2)

Evaluating the Impacts from Treatments, Stages and Genotypes to Diurnal Changing Patterns
Before modeling the diurnal changing patterns, the diurnal patterns between different nitrogen treatments, growth stages and genotypes were compared to decide if any of these factors had significant impact on the diurnal pattern. If not, the data from different treatments, stages and genotypes could be combined for the diurnal pattern modeling. Otherwise, the modeling should be done separately for each different case.
The comparison of the changing patterns was done by applying the dynamic time warping (DTW) method to calculate the similarity between the relative different ratio (RDR) curves from the different plants plots [42]. An RDR curve was calculated to describe the diurnal changing pattern on each day as the % of the change of the phenotyping feature value relative to the feature's value at the reference time point (Equation (4)). For the reference time, we selected the solar noon time as this is the center point of the daytime, and this is when the lowest NDVI value was observed every day [11].

RDR =
Feature t − Feature solar noon Feature solar noon The DTW method was selected as it is an algorithm more commonly used to measure the similarity between two time-series data [42]. More specifically, DTW is a time series alignment algorithm developed to align two sequences of feature vectors by warping the time axis iteratively until an optimal match [43]. Thus, a distance score is generated during the process of alignment, which can be used as the difference between two curves. For example, small distance score means higher similarity between two curves. DTW also allowed non-linear mapping which was appropriate for the purpose of pattern matching. With the distance scores from DTW, the similarity of RDR curves of different plots were quantitively compared and discussed.

Diurnal Patterns Calculation by Time Series Signal Decomposition
Inspired by the idea of time series decomposition method [44], we decomposed the changing signal of each feature into two major parts: day-to-day trend (T t ) and diurnal pattern (D t ). T t is calculated with LOESS (locally estimated scatterplot smoothing) method [45]. By fitting a non-parametric regression curve on the scattered plot of the data, the day-to-day changing trend can be clearly extracted from the raw signal [44]. This trend is majorly reflecting the changes of plant growth stage and general weather conditions over the 31 days of imaging, which would be bothering if it is included in the diurnal pattern modeling. The diurnal component (D t ) was calculated by subtracting the day-to-day trend (T t ) from the raw signal. D t is also called the detrended data. D t contains the higher frequency variance components majorly caused by the short-term impacts such as plant's circadian behavior, sun angle, solar irradiation and temperature changes during the day.

Raw Time Series Signal
The mean curve of D t for 31 days was calculated for the diurnal pattern fitting. Meanwhile, the 95% confidence interval was also calculated to evaluate the consistency of the diurnal patterns across the days.

Diurnal Pattern Fitting
Based on the observation for the diurnal changing pattern of NDVI, a piecewise linear regression function was selected to fit the pattern. By calculating the 1st order derivative of the mean curve, the critical lowest point t a (which perfectly matched the time with the highest sun angle) was derived.
where t is the time offset (in hours) from solar noon. (i.e., t at 12:15 pm is -1.5, on a day and at a location where solar noon is 13:45.). Non-linear diurnal changing patterns were observed for the other features such as RWC prediction and a few single spectral bands. We selected the 2nd order regression model (Equation (7)) for fitting the pattern of those features.

Model Performance Evaluation
The performance of the developed diurnal model was evaluated and compared with the R 2 and RMSE between the prediction results and the desired diurnal changes. Moreover, to further assess the impacts from nitrogen treatments and genotypes, the general diurnal model was tested on each of the six plant plots, and the R 2 and RMSE were obtained for each plot, respectively.

Diurnal Models' Applications
These regression models quantitively summaries the diurnal impacts. For one example, the models can be used to calculate the exact imaging window based on any specific diurnal variance requirement. Besides, the models can also be used to remove the diurnal impact. For example, the NDVI measured at any other time point of the day can be converted to the NDVI at the highest sun angle time with the fitted model.

Results
In this study, the diurnal changing patterns models were built for various plant phenotyping features (including NDVI, RWC, Band670 (Red) and Band760 (NIR)). The detailed procedures and corresponding results for modeling NDVI's diurnal changing pattern was demonstrated and discussed first, as NDVI is one of the most common plant features in remote sensing [46]. The summarized modeling results are reported for the other plant phenotyping features including RWC, Red and NIR, respectively.

The NDVI Diurnal Fluctuations
The change of NDVI values was plotted to show the diurnal fluctuations from the raw imaging data ( Figure 5). Figure 5a is the result for all 31 days. Each gathering group represents one day, and the gaps between groups are from 17:30 until 10:00 next day without imaging data. In this figure, NDVI shows obvious diurnal variabilities with repeatable V-shaped patterns. Figure 5b is the zoomed in diurnal fluctuation in Day 1. It shows that NDVI keeps decreasing until the solar noon time at 1:39 p.m. on that day. This confirmed the similar findings were reported in previous research papers that the NDVI of wheat and soybean were higher at the beginning and end of the day [16,47]. This diurnal change could be due to the combination of the imaging lighting condition change and the plant's physiology changes on a cyclical diurnal nature, which is a central behavior to regulate the susceptibility and responses to biotic or abiotic stresses [48][49][50]. The significance of diurnal variance was assessed by calculating the difference ratio between highest and lowest NDVI during the same day (Table 3). NDVI had the highest diurnal variance of 15.71% for the HN and B73 × Mo17 genotype plot. This high diurnal variance is too severe to be accepted in most plant phenotyping studies [51,52]. This justifies the needs to model the diurnal pattern and calibrate the diurnal impacts.

The Impacts of Nutrient Treatments, Genotypes and Leaf Stages on Diurnal Variation
The relative difference ratio (RDR) curves were compared to study the impacts from treatments, genotypes and growth stages on diurnal variation patterns. For each plant The significance of diurnal variance was assessed by calculating the difference ratio between highest and lowest NDVI during the same day (Table 3). NDVI had the highest diurnal variance of 15.71% for the HN and B73 × Mo17 genotype plot. This high diurnal variance is too severe to be accepted in most plant phenotyping studies [51,52]. This justifies the needs to model the diurnal pattern and calibrate the diurnal impacts.

The Impacts of Nutrient Treatments, Genotypes and Leaf Stages on Diurnal Variation
The relative difference ratio (RDR) curves were compared to study the impacts from treatments, genotypes and growth stages on diurnal variation patterns. For each plant plot, the mean RDR curve was computed from all the 31 RDR curves (one RDR curve per day). In Figure 6, the mean RDR curves from different treatments and genotype plots showed very similar changing patterns. For a more quantitative comparison, the DTW distance scores were calculated to measures the similarity between RDR curves (Table 4). Small distance score means higher similarity between two curves. The distance scores between plots with the same genotype but different nitrogen treatments are between 1.37~3.75. Meanwhile, the distance scores between different genotype plots are 0.33, 0.20 and 0.30 for HN, MN and LN, respectively. Table 5 recorded the DTW distance scores between different leaf stages. Most of these distance scores are below 1. As a summary, comparatively, the nitrogen treatment had higher impact on the diurnal variation than genotype and plant stage. This confirmed the observation in previous study [53]. However, all these factors had minor impacts on the diurnal changing pattern. This allowed to combine the data for building one general diurnal changing model.  HN and B73 × Mo17 HN and P1105AM MN and B73 × Mo17 MN and P1105AM LN and B73 × Mo17 LN and P1105AM  For a more quantitative comparison, the DTW distance scores were calculated to measures the similarity between RDR curves (Table 4). Small distance score means higher similarity between two curves. The distance scores between plots with the same genotype but different nitrogen treatments are between 1.37~3.75. Meanwhile, the distance scores between different genotype plots are 0.33, 0.20 and 0.30 for HN, MN and LN, respectively. Table 5 recorded the DTW distance scores between different leaf stages. Most of these distance scores are below 1. As a summary, comparatively, the nitrogen treatment had higher impact on the diurnal variation than genotype and plant stage. This confirmed the observation in previous study [53]. However, all these factors had minor impacts on the diurnal changing pattern. This allowed to combine the data for building one general diurnal changing model.

Diurnal Changing Pattern
Based on the result in Section 3.2, the data of all six plant plots were combined for modeling the diurnal changing pattern because the nitrogen treatments, genotypes and leaf stages only had minor impacts on the diurnal variation and a general diurnal pattern model is preferred for ease of use.

Decomposition
The changing signal of phenotyping features were decomposed into two major parts: day-to-day trend (T t ) and diurnal pattern (D t ). Figure 7 shows the decomposition result for NDVI. From the trend T t chart, the overall NDVI values increased from the early leaf stage to the late stage. This confirmed with the expectation that as plants grew up, NDVI increases until the reproductive stage [54]. Meanwhile, there were two big dips along the timeline, which matched well with two severe temperature drops in the West Lafayette area. Therefore, the day-to-day trend T t well described the changes of plant stage and plant health conditions.
On the other hand, the diurnal pattern chart in Figure 7 showed clear periodical oscillation with repeated V-shapes. For most days, the diurnal patterns were substantially consistent that the value decreased until reaching the solar noon time. We were not able to collect the complete data for the whole day on 36, 37, 58, 59 and 60 DAP due to the extreme weather conditions. It is interesting to investigate how the real-time weather conditions affect phenotyping results at different leaf stages. In this paper, we only focus on modeling the diurnal changing patterns under "normal" weather conditions. The modeling of the weather condition impacts will be reported in another paper.
To model the general diurnal pattern, the patterns through the 31 days were combined for the mean curve and 95% confidence interval (Figure 8). The narrow 95% confidence interval also provided strong evidence that the diurnal pattern is consistent through the days. As discussed in the methods, because the confidence interval is much lower compared with the mean curve's change over the diurnal cycle, the mean curve (black line) is used for modeling the diurnal pattern. This also helps to "average out" the impacts from abnormal weather conditions in a few certain days. Figure 7. The NDVI of HN and Genotype B73 × Mo17 plot from V4 stage to the R1 stage. The raw NDVI plot was decomposed into the day-to-day trend ( ) and diurnal pattern ( ). The red boxes are the days with incomplete data measurements due to the extreme weather conditions, which are 36, 37, 58, 59 and 60 DAP.
To model the general diurnal pattern, the patterns through the 31 days were combined for the mean curve and 95% confidence interval (Figure 8). The narrow 95% confidence interval also provided strong evidence that the diurnal pattern is consistent through the days. As discussed in the methods, because the confidence interval is much lower compared with the mean curve's change over the diurnal cycle, the mean curve (black line) is used for modeling the diurnal pattern. This also helps to "average out" the impacts from abnormal weather conditions in a few certain days.  To model the general diurnal pattern, the patterns through the 31 days were combined for the mean curve and 95% confidence interval ( Figure 8). The narrow 95% confidence interval also provided strong evidence that the diurnal pattern is consistent through the days. As discussed in the methods, because the confidence interval is much lower compared with the mean curve's change over the diurnal cycle, the mean curve (black line) is used for modeling the diurnal pattern. This also helps to "average out" the impacts from abnormal weather conditions in a few certain days.

Pattern Fitting
Both 1st order and 2nd order piecewise models were tested to fit the diurnal pattern of NDVI. Both had excellent regression performances and there was no significant difference between the two models. Therefore, we adopted the simpler 1st order model, and we concluded that the NDVI diurnal changing pattern can be described by a linear model. The fitting equations in Equation (8) accurately predicted the observed diurnal pattern. The predicted NDVI values and the measured values had R 2 of 0.99 and RMSE of 0.0012. Figure 8 illustrates the detailed fitted result. To test the validation of the model for different genotypes and nitrogen treatments, the fitted equation was applied on each different field plot, respectively. The model maintained accurate performance in all cases (Table 6). It proves this general diurnal pattern model is valid for plants from different nutrient treatments and genotypes. It can be noted that the critical time in Equation (8) is 13:42, which is within the actual solar noon time ranging from 13:35 to 13:45. Therefore, we conclude that the solar angel is one of the major impacts to the diurnal variance, so Equation (8) should be shifted based on the actual solar noon time on each day (replacing "13:42" in Equation (8) with the actual solar noon time on the imaging day).
NDV I(t) = −0.012 × t + NDV I solar noon i f t ≤ 13 : 42 0.010 × t + NDV I solar noon i f t > 13 : 42 (8) where t is the time offset (in hours) from solar noon. (i.e., t at 12:15 pm is -1.5, on a day and at a location where solar noon is 13:45.) NDV I solar noon is the NDVI measured at solar noon.

Applications of the Model
The diurnal changing pattern model can be utilized to determine the proper imaging window based on any specific quality requirement. In Figure 9 and Table 7, the proper imaging windows were calculated based on the different thresholds tolerated diurnal variance. For example, in order to control the NDVI's diurnal error within ±0.03, the imaging needs to happen between 10:55 and 16:30.
Remote Sens. 2021, 13, x FOR PEER REVIEW 16 of 21 Figure 9. The intersections between the adjusted NDVI's diurnal changes and three different thresholds for proper imaging windows. The x-axis is the diurnal time in a unit of hour. The y-axis is the adjusted NDVI's diurnal changes when adjustment at solar noon to be 0.   Besides the imaging time window, the diurnal changing model can also be used to easily calibrate the diurnal variances. For example, the NDVI measured 2 h before the solar noon time should be decreased by 0.024 (0.012 * 2), so the calibrated NDVI is as if it is taken exactly at the solar noon time. It'll be important to run this calibration process, which will help to remove the significant diurnal variance and enable phenotyping researchers to do "apple-to-apple" comparison between the NDVIs of the different plots.

Other Image-Derived Phenotyping Features
Besides NDVI, the diurnal changing patterns of other plant phenotyping features, including RWC prediction, Red and NIR in the spectra, were modeled. From Figure 10a,c,e, the result showed consistent diurnal patterns of these features with clear mean curve and narrow 95% CI.    Second order models were selected, because these other features showed non-linear diurnal changing patterns. The diurnal patterns were fitted and plotted in Figure 10b,d,f. Equations (9)-(11) are the fitted equations. All the models showed excellent description of the observed diurnal patterns. The R 2 and RMSE between the models' prediction and the measured features are summarized in Table 8. Similarly, these models have the potential to help improving the efficiency and data quality in phenotyping practices. Moreover, successful application of the same proposed diurnal pattern modelling method on the RWC, Red and NIR indicated that this method has potential to be generalized well on other phenotypic features which can be further explored. N IR(t) = −0.0027 × t 2 − 0.031 × t + N IR solar noon i f t ≤ 13 : 45 : 00 0.00020 × t 2 + 0.040 × t + N IR solar noon i f t > 13 : 45 : 00 (11)

Strengths
Benefited from the unprecedented high imaging frequency of the new field hyperspectral imaging gantry at Purdue, with over eight thousand repeated images over one growing season, the detailed diurnal patterns of plant phenotyping features, including NDVI, RWC, and two individual spectral bands (Red and NIR), were quantitatively described and modeled for the first time. Moreover, most of previous studies reported the diurnal changes directly from the raw phenotyping features (e.g., NDVI, chlorophyll content) [11,16,23]. However, these results are inconsistent due to the changes of plant growth stage and weather condition. To address this issue, our proposed approach incorporates a modified time series decomposition method to filter out the signal variance caused by the changing plant stage and weather conditions, so the later modeling process can exclusively focus on the higher-frequency diurnal changes.

Limitations and Future Work
In this experiment, all the imaging data was collected from the same field imaging gantry facility. Extra validation is needed for the other remote sensing platforms such as UAV. In addition, the diurnal patterns are reported based on corn plants, which may limit the scope of application. It is necessary to conduct more tests including more diverse plant species (e.g., soybean, wheat, and rice).
In the 2020 season, the developed diurnal models will be externally validated with the field gantry's images. The model will also be validated on a UAV hyperspectral imaging platform in the same field. Finally, the diurnal changing patterns on other plant species, such as soybean wheat, and rice will be also explored and compared.
Besides the diurnal effect, the changing environmental conditions in the field such as cloud coverage, temperature, wind speed and so on all showed strong impacts on the phenotyping results. The data has shown that these effects could be precisely modeled as well. The environmental calibration models for plant remote sensing will be reported in future papers.

Conclusions
In this study, diurnal changing patterns in crop remote sensing images were quantitatively investigated with the proposed novel modeling approach. The diurnal patterns for four selected plant phenotyping features, including NDVI, RWC, two individual spectral bands (Red and NIR), were successfully modeled. The results showed that the NDVI presents a repeatable V-shaped diurnal changing pattern: it linearly decreases by 0.012/h before the solar noon time and increases by 0.010/h thereafter. On the other hand, the RWC, and the two individual spectral bands shows non-linear diurnal changing patterns, where RWC and Red change in inversed V-shaped patterns, and NIR changes in a normal V-shaped pattern. The diurnal patterns lightly vary between the nitrogen treatments, genotypes, and plant stages, but the difference is minor and ignorable in the modeling. With the result of this work, Ag remote sensing users will be able to more precisely understand the deviation/change of crop feature predictions caused by the specific imaging time of the day. The reported diurnal models can also be used to calibrate the remote sensing result so to remove the diurnal impacts.
Author Contributions: D.M. carried out the data collection, data analysis, and drafted the manuscript; J.J. contributed the research idea, led the designing, fund-raising and construction of the field imaging facility, designed the experiments and revised the manuscript; H.M. and M.R.T. guided the experimental design; T.U.R. and L.Z. made substantial contributions to system construction and data collection. All authors have read and agreed to the published version of the manuscript.
Funding: This research was funded by Sumitomo Chemical, grant number 16121941 and The APC was funded by Purdue University.

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