Effect of the Solar Zenith Angles at Different Latitudes on Estimated Crop Vegetation Indices

: Normalization of anisotropic solar reﬂectance is an essential factor that needs to be consid-ered for ﬁeld-based phenotyping applications to ensure reliability, consistency, and interpretability of time-series multispectral data acquired using an unmanned aerial vehicle (UAV). Different models have been developed to characterize the bidirectional reﬂectance distribution function. However, the substantial variation in crop breeding trials, in terms of vegetation structure conﬁguration, creates challenges to such modeling approaches. This study evaluated the variation in standard vegetation indices and its relationship with ground-reference data (measured crop traits such as seed/grain yield) in multiple crop breeding trials as a function of solar zenith angles (SZA). UAV-based multispectral images were acquired and utilized to extract vegetation indices at SZA across two different latitudes. The pea and chickpea breeding materials were evaluated in a high latitude (46 ◦ 36 (cid:48) 39.92 (cid:48)(cid:48) N) zone, whereas the rice lines were assessed in a low latitude (3 ◦ 29 (cid:48) 42.43 (cid:48)(cid:48) N) zone. In general, several of the vegetation index data were affected by SZA (e.g., normalized difference vegetation index, green normalized difference vegetation index, normalized difference red-edge index, etc.) in both latitudes. Nevertheless, the simple ratio index (SR) showed less variability across SZA in both latitude zones amongst these indices. In addition, it was interesting to note that the correlation between vegetation indices and ground-reference data remained stable across SZA in both latitude zones. In summary, SR was found to have a minimum anisotropic reﬂectance effect in both zones, and the other vegetation indices can be utilized to evaluate relative differences in crop performances, although the absolute data would be affected by SZA.


Introduction
Remote sensing applications utilizing light interactions (from visible to shortwave infrared) with crop/plant leaves can be used to monitor crop responses. The biochemical and biophysical factors within a leaf, including cellular organization, determine its optical properties, such as reflectance [1]. However, multi-angular remote sensing can affect the reflectance data due to the anisotropic reflectance property, which affects both the radiance intensity and spectral distribution captured by the sensor [2]. Many approaches to represent anisotropic reflectance using the bidirectional reflectance distribution function (BRDF) have been developed, and it has progressively gained recognition in the remote sensing community [3]. BRDF is crucial for radiometric adjustments, especially in terms of removing disruptions that could limit image interpretation skills [4].
Different mathematical methods based on BRDF have been devised to correct the anisotropic reflectance property [5]. For example, Latifovic et al. [6] compared four BRDF models (modified Walthall, semi-empirical Roujean, Ross-Li, and a new nonlinear temporal angular models) using satellite data. Renhorn et al. [7] reported the polarimetric BRDF model utilization for the multispectral data, which was composed of four parameters that account for various surface structure types by applying a generalized Gaussian distribution. Wierzbicki et al. [4] developed the dark-object subtraction technique based on the modified Walthall's model to compensate radiometric disturbances in near-infrared images acquired from an unmanned aerial vehicle (UAV). Nonetheless, quantification and use of BRDF models in remote sensing applications can be challenging due to data accuracy [4], parameter selection [8], and high measurement costs [7].
UAV applications in high-throughput phenotyping have become increasingly prominent and common in recent years due to their accessibility and affordability [9], combined with the handiness, flexibility, and ability to collect high-quality non-destructive data [10]. In lentil, Ahmed et al. [11] utilized UAV/multispectral imaging to detect and assess breeding plots using normalized difference vegetative index (NDVI) images. Similarly, Quirós et al. [12] found that NDVI created from UAV-based multispectral imaging consistently correlated with yield potential in pea across different breeding trials. Zhang et al. [13] demonstrated that the approach of multispectral imagery from UAV could be applied to monitor yield traits in breeding trials of pea and chickpea using vegetation indices (VI) such as NDVI and soil-adjusted vegetation index (SAVI). In pea, Zhang et al. [14] reported that the flowering features have a strong correlation with digital crop traits from multispectral images. In rice, Kang et al. [15] integrated UAV multispectral imagery with artificial neural network (ANN) analysis and indicated that medium resolution imaging spectrometer terrestrial chlorophyll index (MTCI) was the most important predictor for yield and protein content. Ogawa et al. [16] reported four quantitative trait loci (QTL) for vegetation fraction using RGB (red, green, and blue)/UAV imagery.
Several studies have highlighted the capability of UAV-based multispectral imagery for crop/plant phenotyping schemes in field conditions. Nevertheless, the influence of solar zenith angles (SZA) on digital crop traits extracted from UAV data has not been studied. Therefore, this reported work serves as a case study, where we assessed the stability of digital traits such as vegetation indices across different solar zenith angles in two latitude zones, utilizing breeding trials. The major goal was to assess the role of SZA on the variability in the vegetation indices based on its relationship to ground reference data, such as seed/grain yield acquired from multiple breeding trails at two latitudes. The specific objectives were to: (1) evaluate and identify key vegetation indices associated with the yield potential in the breeding trials; and (2) assess the stability of the vegetation indices across different solar zenith angles at two latitude zones with distinct climatic environments. The SZA is a function of both location and time of the year/season; nevertheless, the SZA can be computed based on the location and season, and the findings reported herein may be relevant. We anticipate that these findings will be critical for utilizing such data in phenotyping and precision agriculture applications.

Study Locations
Four breeding trials were utilized to evaluate the effect of solar zenith angles on vegetation indices extracted from individual plots at each latitude zone. Two pulse trials were located at Genesee, Idaho, United States (representing high latitude, Figure 1a)  In these trials, 20 lines of rice (Oriza sativa L.), 12 lines of IR64 and eight lines of cv 'Ciherang', were planted on 16 May 2020, transplanted on 10 June 2020, and harvested on 25 September 2020. Rice trials were arranged in randomized block designs with four replications and plots of 2 m × 2 m. Transplanting was carried out with 20 cm between plants, 25 cm between rows, and 50 cm between plots. Agronomic management (irrigation, fertilization, and pest control) was performed in accordance with standard paddy rice production procedures.

UAV Systems and Sensors
In pulse trials, images were acquired at the early pod development stage. Our previous work [12,13] reported a strong and significant correlation between seed yield and VIs at flowering/early pod development stages. Multispectral images were acquired using a Quad Sensor (Sentera, LLC., Minneapolis, MN, USA) attached to a quadcopter UAV (DJI Inspire 1, DJI Inc., Shenzhen, China). Another quadcopter UAV (DJI Phantom 4 Pro, DJI Inc., Shenzhen, China) was utilized to capture RGB images (Figure 1c In rice, images were acquired at the late milky stage. Our previous research [17] highlighted the milky stage as a promising phenological stage for rice grain estimation in different rice trials. Multispectral images were acquired using an Altum sensor (Altum, MicaSense Inc., Seattle, WA, USA) attached to an octocopter UAV (DJI Matrice 600, DJI Inc., Shenzhen, China) using SkyPort adapter, which allowed the operator/pilot to trigger and control the camera through the flight mission. The RGB images (12.35 MP) were captured with a quadcopter UAV (DJI Mavic Pro Platinum, DJI Inc., Shenzhen, China) ( Figure 1d). The Altum sensor included five spectral channels blue, green, red, red edge, and nearinfrared (with center and bandwidth of 475 ± 32, 560 ± 27, 668 ± 16, 717 ± 12, and 842 ± 57 respectively). Both multispectral and RGB images were captured on 5 September 2020, at the following times and corresponding solar zenith angles: 08:30 (53.  (Figure 2c). UAV flight missions were created using Ground Station software (SZ DJI Technology Co. Ltd., China) in a single pattern, speed of 3 m s −1 , 80% of horizontal and 75% vertical overlaps, and at 45 m above ground level resulting in spatial resolution of 1.05 cm per pixel for multispectral sensor and 0.60 cm per pixel for RGB sensor. A MicaSense reflectance panel (MicaSense Inc., Seattle, WA, USA) was used to convert digital numbers into reflectance values for radiometric calibration. For all trials, light intensity data (2-4 data points per hour) and photosynthetically active radiation were collected using a quantum sensor (MQ-100, Apogee Instruments, Inc. USA) or an optical light meter (Lux-1330B, China) at different SZA to capture solar intensity at different time points (Figure 2b,d).

Ground-Truth Measurements
Seed/grain yield (GY) of both rice and pulse trials served as the primary groundreference data. In addition, plant height (PH), tiller number (TN), and above-ground dry biomass data (DB) were acquired for rice trails. For rice, 28 plants per genotype were randomly selected (seven plants per replication), wrapped in paper bags, and taken to the seed lab. The grain and vegetative tissues (stems and leaves = biomass) were separated into independent paper bags; grain was dried to 14% of moisture, and vegetative tissues were dried for 82 h until a stable weight was obtained. Finally, grain weight, number of tillers, and weight of the vegetative tissues were recorded per plant. For pea and chickpea, the entire plot was harvested using a small plot combine to acquire seed weight. The summary statistics of agronomic traits for both pulse and rice trials are presented in Table 1.

Image Features Equation
Normalized Difference Vegetation Index In rice and pulse trials, RGB images were processed in Pix4Dmapper to generate digital surface models (DSM), which reported elevation data in meters above mean sea level. The topography of the terrain as Digital Terrain Model (DTM) data were created based on the interpolation of elevation data over bare soil points on DSM using open-source software Quantum Geographic Information System (QGIS) (QGIS.org, 2021, version 3.10.4).
The DSM and DTM as .tif files were imported into MATLAB to calculate and extract the CH (Table 2).

Data Analysis
RStudio software [28] (Version 1.4.1106, 2021, Boston, MA) was used to evaluate the linear relationship between agronomic traits and vegetation indices using Pearson's correlation analysis. The VI data were summarized using boxplots (Q3, median, Q1, and suspected outliers) to visualize the variability across SZA and locations. Transformationestimating functions were used for data normalization using the bestNormalize package. Analysis of variance (ANOVA) linear model was conducted to assess the effect of SZA on VI at each location; for this model, the dependent variable score (Y ij ) is expected to be the result of the added effects of the factors treatment (τ i ), block (β j ), and error term ( ij ); the effects are added to the global mean (µ): Based on the results, post-hoc pairwise comparisons (α = 0.05) of group means were performed using Tukey Honest Significant Difference (HSD) to test significant differences among SZA in each location. Table 3 shows a higher data dispersion in rice trials (33 to 70%) than pulse trials (23 to 30%). In both rice and pulse trials, a significant effect of SZA was detected (p-value < 0.0001 for all bands in all trails, except green in Ciherang, where p-value was 0.001).
In rice trials, there was a statistically significant correlation between SR, NDRE, and GNDVI and grain yield with minor differences across SZA (Figure 4a-f). The GNDVI displayed moderate to high significant correlations with GY in the IR64 (r = 0.52 to 0.59) and Ciherang lines (r = 0.64 to 0.70) across SZA as presented in Figure 4c,f. In addition, the matrix inter-correlations between SR (e.g. SR at one angle compared to SR from another angles) and NDVI for pulse lines (Figure 3a,d), and SR and GNDVI for rice trails (Figure 4c,f) at multiple SZA, were between 0.74 to 0.99, which indicated stability across SZA. Figure 5 shows correlation plots between SR or GNDVI with GY for chickpea, pea, and rice plots across SZA. The correlation with seed/grain yield and other agronomic traits are summarized in Supplementary materials, Figures S1-S4.

Correlation Analysis with Dry Biomass and Tiller Number
Biomass (the vegetative tissues) is a reliable indicator of rice production, thus, a strong correlation (r = 0.85 to 0.86, p ≤ 0.001) between rice grain yield and dry biomass was observed. Pearson correlation demonstrated that SR, NDRE, and GNDVI showed stable and highly significant correlations (p-values ≤ 0.001) in both trials with dry biomass across SZA (Figure 4a-f). For the IR64 lines, VIs displayed a good correlation (Figure 4c) with biomass, while for the Ciherang lines, a higher correlation (Figure 4d) with biomass was observed, especially with SR. The relationship among SZA between tiller number and VI indicated a significantly moderate correlation with SR, NDRE, and GNDVI in the IR64 lines (Figure 4a-c). In contrast, TVI showed the best correlation (r = 0.46 to 0.55) with the tiller numbers for the Ciherang lines (Supplementary materials, Figure S3d). Table 4 presents the correlations between UAV-based plant height estimates and ground-truth reference data and VIs across SZA. In pulses, NDVI, GNDVI, and TVI showed a significant correlation with plant height data across most SZA. The highly significant correlations (p-values ≤ 0.001) for plant height were found in pulses at 09:00 (63.82 • ), 11:00 (43.80 • ), and 13:00 h (28.00 • ) with NDVI. Moderate to a high correlation between measured and estimated plant height was observed for the IR64 and Ciherang trials across SZA. For the rice trials, except for Ciherang at noon (3.11 • ), NDRE and SR demonstrated a significant correlation with plant height across SZA. Based on ANOVA analysis, there was the SZA effect (p-values ≤ 0.0001); however, the post-hoc test indicated reasonable stability in data acquired at most SZA (Figure 6a-c).

Influence of SZA Based on Analysis of Variance and Post-Hoc Test
The VI data that were substantially correlated with grain yield were selected to perform the ANOVA across SZA. The mean differences of SR, OSAVI, NDRE, and GNDVI across SZA and locations (p-values ≤ 0.0001) strongly indicated that the mean of at least one group significantly differed from other groups (Figures 7-9). This indicates that the Lambertian property may not be present on the VI acquired from pulse and rice crop surfaces. Tukey's multiple comparison test showed the different groups (Figures 7-9). SR and OSAVI were used for pulse trials for this test, while NDRE, GNDVI, and SR were used for rice trials. In pulse trails, in comparison to the SR, OSAVI showed more variation across SZA (Figure 7); indicating a more anisotropic effect on OSAVI associated with solar view angles, especially when the images were captured at the earlier and later times (Figure 7b, d); however, OSAVI showed low mean square error (0.00035 to 0.00077). This comparison indicates that SR can be measured over a broader solar zenith angle range, 08:00 h (73.98 • ) to 14:00 h (25.79 • ) (Figure 7a,c), Tukey test analysis on NDRE, GNDVI, and SR of the IR64 lines also showed statistically different groups among solar view angles; in general, noontime had low VI values across SZA. GNDVI had the highest values at early and late time points (Figure 8c,d). We observed that trends of NDRE and SR indices at most SZA were similar, except at 12:00 (3.11 • ) and 14:00 (29.07 • ) (Figure 9a,e). In contrast, NDRE and SR over the Ciherang lines revealed few groups compared to the IR64 lines. This result indicates that the VI values were statistically similar across SZAs (Figure 9a,e). Based on this, we can infer the minimum anisotropic effect.

Discussion
The models developed to estimate BDRF are found to offer robust and effective compensation of reflectance anisotropic [4,8]. Nevertheless, such measurements can still be challenging due to data accuracy [4] and parameter selection [8], in addition to variability within the crop. In this study, we evaluated the performance of some of the VI widely used for plant breeding to identify VI that may not need anisotropic correction. As expected, we found the effects of the SZA on SR, OSAVI, NDRE, and GNDVI at both high and low latitudes (Figures 8 and 9). In the tallgrass prairie, Middleton [29] reported significant dependence of SR and NDVI on the SZA. In this study, the variation caused by the SZA was found to be the least for SR and the greatest for other VI in high and low latitudes. This result is supported by Hashimoto et al. [30], who simulated the reflectance on the rice paddy field, indicating the least saturation of SR compared with NDVI. Ma et al. [31] found that the NDVI data extracted from savanna were more sensitive than enhanced vegetation index (EVI) to the SZA effects. Since the absolute VI data is affected by the SZA, the data need to be standardized prior to the use of such data in crop models and similar applications. Moreover, consistency during data acquisition should be maintained for proper interpretation of the data. The VI data that is independent of the anisotropic reflectance effect may also be used to create stable datasets. Nevertheless, further validation is needed across crops, phenological stages, and latitudes to explore these solutions. Middleton [29] and Ishihara [32] found varied results when assessing VIs data extracted from imagery from prairie and croplands.
In addition to a stable correlation between SR and GY across SZA, a very strong correlation between SR data at different SZA was found in this study. Naito et al. [16] reported a similarly strong correlation between SR and GY in paddy rice as well, based on data collected around noon. Even though SR, OSAVI, NDRE, and GNDVI data were statistically different across SZA, they showed a stable correlation with grain yield. In general, SR data exhibited low variability (mean square error = 0.013 to 0.37) and better consistency. Similar results were found with NDRE data, especially for Ciherang rice lines. Based on the correlation analysis with seed/grain yield, the results indicate that for phenotyping applications in crop breeding programs, where relative differences between crop performance across different varieties are evaluated, multiple vegetation indices can be utilized, even if there are differences in the SZA during data acquisition. Although the results are encouraging, further validation across breeding programs is needed. Moreover, standardization of data acquisition is necessary to develop crop growth and development curves using remote sensing data.

Conclusions
In this study, we found that the evaluated vegetation indices are significantly affected by different SZAs. There were two major findings: (1) The evaluated VIs maintained a consistent correlation with seed/grain yield in all four trials at both latitudes (low and high). The evaluated vegetation indices can be applicable for phenotyping applications in breeding programs to study the relative differences between crop varieties; and (2) based on the analysis to find the best VI with minimum bidirectional effect (low variance) and high correlation with grain yield across SZAs, SR was identified. Thus, SR has the potential as a yield indicator or predictor across solar zenith angles for a wide range of remote sensing applications without anisotropic reflectance correction.
Nevertheless, further work is needed, using different crops, phenological stages, latitudes, and different sensor angles. Such research will be useful to further validate the findings from this study and identify similar/different vegetation indices with minimum anisotropic reflectance effect, allowing image acquisition across broader solar zenith angles, which is essential for the interpretation of images in precision agriculture and crop modeling applications.