Examination of appropriate observation time and correction of vegetation index for drone-based crop monitoring

Use of information and communication technologies, as well as robotics, routinely saves labor and refines agricultural tasks; thus, innovative “smart farming” to maintain and enhance the quality of crops can improve the sustainability of agriculture. When managing crop growth using remote-sensing drones, the normalized difference vegetation index ( NDVI ) ̶ used to assess growth ̶ typically changes depending on sunlight conditions. In this study we have attempted to develop an empirical correction to correct for differences in sunlight conditions in drone NDVI images of paddy rice. Based on observations using a field sensor installed in a paddy field, and considering the effects of morning dew, we determined that 10:00 AM is the most appropriate time for NDVI observations in paddy rice, when the morning dew has largely evaporated. This observation time differs from that used in the radiative transmission models described in previous studies. In the drone observations, sections with lower NDVI were more strongly affected by solar altitude, and thus by time of day. Therefore, we found that when correcting NDVI according to sunlight conditions, it is necessary to adjust the correction parameters depending on the NDVI values. Based on the aforementioned results, we corrected the drone-observed NDVI and succeeded in mitigating the decline in NDVI value associated with changes in sunlight conditions, in terms of both NDVI values and NDVI images, within plots established in the ( Inoue and Yokoyama, 2017; Mee et al. , 2017 ) . The use of satellite data for farming guidance, such as optimizing harvest timing and fertilizer application ( Sakaiya et al. , 2016 ) , is an example of the


Introduction
Changes over recent years in the global climate, typically warming, have significantly changed the environments in which agricultural crops are grown. These changes typically diminish crop production and are predicted to become more widespread in the future IPCC, 2014 . Improving agricultural sustainability is a global challenge, and the most important mission in different production regions will be to develop innovative methods of cultivation management to cope with changes in the growth environment.
In the 20 years leading up to 2015, the agricultural workforce in Japan shrank by approximately 50 , and the area of abandoned fields and rice paddies now exceeds 400 000 ha MAFF, 2016 . Furthermore, the remaining workforce is aging considerably, with 66.3 years being the current average age of agricultural workers MAFF, 2016 . Consequently, it is necessary not only to secure crop quality and yield, but also to protect the sustainability of Japan's agriculture. Preserving and managing future agriculture and farmland is of immediate importance, with implications for the sustainability of society as a whole.
Given that robotics and information/communication technologies are now routinely used to save labor and refine agricultural tasks, maintaining and improving the quality of crops with innovative "smart farming" could significantly enhance the sustainability of agriculture. Smart farming activities may potentially save labor and lighten agricultural workloads, as well as attract new agricultural workers and allow existing farmers to improve their technical abilities in cultivation Inoue and Yokoyama, 2017 . In this regard, the importance of remote sensing RS technology contactless sensing and image analysis that can measure and estimate farmland targets has become widely recognized MAFF, 2019 .
RS developed as an information-processing and analysis method for observational data from artificial satellites; numerous studies have employed RS data for vegetation monitoring. Among the potential applications of RS-based vegetation observations, the monitoring of agricultural croplands presents practical challenges regarding image analysis and information processing technology Saito et al., 2017 . With RS, visible to near-infrared spectrum information can be used to measure and estimate plant structures of above-ground vegetation, such as crop biomass and leaf area index LAI , as well as the nutrient content chlorophyll, nitrogen, and phosphorus of crop plants Inoue and Yokoyama, 2017;Mee et al., 2017 . The use of satellite data for farming guidance, such as optimizing harvest timing and fertilizer application Sakaiya et al., 2016 , is an example of the practical use of RS in agriculture. However, RS by satellites has certain practical limitations. For example, using optical satellite imagery, the ground surface is generally not observable in the presence of clouds, and because of observation frequency the specific number of days between images of the same location , it is difficult to observe crops at the right time. By contrast, low-altitude RS using drones drone RS delivers a highly maneuverable platform that can make observations at the right time Hama et al., 2020b . With regard to RS using drones, we have previously shown that the normalized difference vegetation index NDVI can be used as an explanatory variable to estimate such paddy rice growth parameters plant structures as plant height and LAI e.g., Hama et al., 2016;2018 . In recent years, other groups have similarly used the vegetation index obtained by drones to estimate crop yield e.g., Zhou et al., 2017;Herrmann et al., 2020 . Studies have also been conducted to assess the development of detailed fertilization standards by using vegetation metrics based on the relationship between these indices and crop growth e.g., Guan et al., 2019;Mochizuki et al., 2020 . There are, however, concerns that sunlight conditions may affect the results of vegetation indices obtained using drone RS, depending on weather conditions and timing. Accordingly, as drones are increasingly making low-altitude RS easier and more effective, the influence of sunlight conditions on the quality of collected data warrants further consideration. In the aforementioned studies Hama et al., 2016;Zhou et al., 2017 , the observation time was uniformly set between 10 : 00 and 14 : 00. However, in practical application, it is not always possible to set the observational time uniformly when monitoring agricultural sites. In practice, vegetation index values do not remain fixed throughout the day, given that as solar altitude the angle of the sun above the horizon increases, vegetation index values are shown to decrease Cogliati et al., 2015;Rahman et al., 2015 . Furthermore, this effect is also influenced by the weather, with solar altitude affecting the vegetation index more on clear days when a greater proportion of the direct components of solar radiation reach lower altitudes and a smaller effect on cloudy days when the dispersion of sunlight components is greater Ishihara et al., 2015; 2019 . In addition, vegetation density also has an effect, with areas of greater density being found to be affected to a lesser extent by solar altitude Ishihara et al., 2015 . The effect of sunlight conditions on the NDVI can be explained thus: the higher the solar altitude, the greater is the proportion of incoming light reaching the ground surface beneath the vegetation layer; under these circumstances, the ground has a greater influence on the reflected light being observed from above Ishihara et al., 2015 . The diffusion of incoming sunlight is greater on cloudy days, and the influence of sunlight conditions is reduced at higher vegetation densities. Moreover, these effects of sunlight conditions are apparent even when converting observational data to reflectance values.
As previously mentioned, it is essential to correct for these effects of sunlight conditions when using vegetation indices to estimate growth, and when using RS data as a basis for determining growth management tasks such as fertilization control. One previous study applied corrections based on estimating the LAI from machine-learned or multiple regression models that use the solar zenith angle, ratio of diffuse solar radiation, canopy reflectance, understory reflectance, or vegetation index as explanatory variables Hashimoto et al., 2020 . This method has obtained sound results; however, it does not mitigate the effects of sunlight conditions on the vegetation index itself.
Furthermore, when using drone RS, it is difficult to capture accurate images of the reflectance and to utilize spectral data; thus, it is preferable to use normalized indices, including NDVI, which apply the same wavelength bands in the numerator and denominator section 2.1 below . When using established vegetation indices such as the enhanced vegetation index EVI Liu and Huete, 1995 , there is a greater effect compared to NDVI attributable to irregularities in brightness due to plant shadows, which in turn increases the influence of factors that are unrelated to vegetation status Matsushita et al., 2007 . In practice, given that the use of a normalized vegetation index is assumed, methods for correcting the vegetation index itself are also considered, even when no method is used to correct for the individual wavelength bands.
Therefore, in this study, we focused primarily on correcting for sunlight conditions in the NDVI drone-captured images of paddy rice. Specifically, we sought to develop an empirical correction formula in which the amount of correction corresponds to variations in growth. With practical use in mind, we strived to produce a corrective method that is simple and requires as few parameters as possible to make the corrections.

Materials and Methods
To achieve the stated objectives, the present study was conducted in 2020 according to the following process.
1. Analyze the appropriate time during the day for observations and set standard values for making corrections see section 2.1 below .
2. Monitor paddy rice when close to the period of fertilization and when close to the heading stage; develop an empirical formula to correct NDVI changes contingent upon changes in sunlight conditions see section 2.2 .

Study sites and field observation equipment
For the purposes of the present study, we used a paddy field located in Sakado, Saitama Prefecture 35 58 N, 139 23 E as the experimental field Fig. 1 . A ground control point GCP was installed at each of the four corners of the field, and the position of each was determined through measurements made with a total station theodolite Ushikata Co. TEO-Ray . "Koshihikari" Oryza sativa L. cv. Koshihikari rice was transplanted on May 22 nd ; the panicle formation stage was on July 11 th , the heading stage was on August 5 th , and the crop was harvested on September 10 th . To observe differences in NDVI results due to variations in growth, we established two 1 m 1 m plots for rice exhibiting dense, moderate, and sparse growth, respectively totaling six plots. We also installed a pyranometer Climatec CHF-SR05-DA1 and leaf wetness sensor Climatec C-LWS in the experimental field. The observation interval was set at 10 min.
To determine an appropriate time for daily observations, we installed field-use spectral reflectance sensors Decagon Devices SRS-Nr Field Stops Sensor, SRS-Ni Hemispherical Sensor in the center of the field, adjusted to a height of 1.5 m from the ground surface Picture 1 . By attaching two sensors that can detect near-infrared 810 2 nm and red 650 2 nm wavelength bands, one facing upward and the other downward, the sensor was able to continuously detect the reflectance-based NDVI Formula 1 .
where NDVI field is the field sensor-based NDVI, and NIR field and Red field are the field-sensor-based near-infrared NIR band and red band reflectance, respectively. The observation interval was set at 10 min, and observations were made from the end of planting until harvest. As NDVI represents paddy rice growth parameters plant structures , large and variable day-to-day fluctuations are undesirable. Thus, the appropriate observation time will be when the NDVI time-series representation is smooth. To examine the appropriate observation time, a curve smoothed by locally weighted scatterplot smoothing LOWESS was created of NDVI field against observation time. Then, the smoothness of the time-series changes was evaluated as the coefficient of determination R 2 analysis of the LOWESS curve against the observed NDVI field .

Drone image acquisition
In this study, we used a DJI Phantom 4 Multispectral P4M drone, which can be used to detect five wavelengths: blue 450 16 nm , green 560 16 nm , red 650 16 nm , red-edge 730 16 nm , and near-infrared 840 26 nm . It is also equipped with a spectral sunlight sensor located at the top of the drone, which records solar irradiance. This solar irradiance was used in converting the digital number DN to signal values, as described below. Aerial photography was conducted using P4M at roughly one-hour intervals on July 22 nd and August 5 th , 2020. July 22 nd between the panicle formation stage and the heading stage was an important time for observing paddy rice growth and using the results to adjust fertilization regimes. August 5 th heading stage was important for estimating the yield or protein content of brown rice by using the drone-based observation data Hama et al., 2020a;2020b . We developed a correction formula by using the data obtained on July 22 nd , and verified the formula using the data from both July 22 nd and August 5 th . The weather on July 22 nd was clear with occasional clouds, and was generally clear on August 5 th . DJI GS PRO iOS software was used to fly the P4M at an altitude of 50 m above the ground, with the flight-path adjusted for 90 image-height overlap and 70 image-width overlap. The camera was set perpendicularly, and photographs were taken at two-second intervals. Figure 1 shows the flight path of the drone when the sun was hidden by clouds. When the sun was unobscured, the GS PRO automatically recommended a flight path. At 50 m altitude, the ground resolution of the captured images was approximately 2.6 cm.

Data and analysis 2.3.1 Preprocessing of drone imagery
Although each of the images captured using the P4M had a recorded DN, the calculated value of the vegetation index varied according to different exposure times and different target brightness. Thus, it was necessary to convert the DN values to spectral reflectance values or to image data that would consistently represent the reflectance. The P4 Multispectral Image Processing Guide DJI, 2020 for the P4M explains the process of converting DN data to consistent signal values that account for the effects of exposure time, gain, solar irradiance, and so on. This pixel-value conversion can be performed automatically using Structure-from-Motion/Multi-view Stereo SfM-MVS software such as DJI Terra and Pix4Dmapper Pix4D, Lausanne, Switzerland , enabling the creation of a single orthomosaic from multiple images captured by the drone.
However, as we used Agisoft Metashape Professional v1.6.0; Agisoft LLC, Russia in the present study, we were unable to automatically convert pixel values. Thus, following the P4 Multispectral Image Processing Guide, we converted the DN values to signal values see the P4 Multispectral Image Processing Guide for details on the conversion method .

Image processing
We used the images that had been converted into signal values to create an orthomosaic by using Agisoft Metashape Professional for each wavelength band. By subjecting the orthomosaic to analysis by using the open-source Geographic Information System GIS software QGIS www.qgis.org , we calculated the NDVI and the average values within the six evaluated plots. NDVI was calculated using the following formula: where NDVI P4M, NIR P4M, and R P4M are the P4M-based NDVI, NIR, and red signal values, respectively.

Correction of the vegetation index
We used the following process for NDVI correction, taking into account sunlight conditions. Initially, we separated the global solar radiation at the time of observation into direct and diffuse components, based on Inanuma and Takeda 2002 method Formulae 3 to 7 . Here, θ is the solar altitude at the time of observation radians , H D is the solar altitude at the time of observation , Kt is the Picture 1. The spectral reflectance sensor installed in the experimental field. SRS-Ni hemispherical sensor was facing upward and SRS-Nr field stops sensor was facing downward.
clearness index, I G is global solar radiation W/m² , I S diffuse solar radiation W/m² , I D is the direct solar radiation W/m² , Q is the global solar radiation outside the atmosphere W/m² , and C Ri is the cloud ratio. Next, we performed NDVI corrections according to the following procedure Formulae 8 to 11 .

Cf = H D H o I D /I G 8
NDVI std / NDVI obs = S Cf + 1 9 S = a NDVI std + b 10 NDVI std = NDVI obs 1 + b Cf / 1 a NDVI obs Cf 11 Here, Cf and S are the NDVI correction variables, H D is the solar altitude at the time of observation , H o is the assumed standard solar altitude , a is the slope, b is the intercept, NDVI std is the NDVI at the assumed standard time, and NDVI obs is the NDVI at the time of observation. The assumed standard solar altitude is the solar altitude at the appropriate time for observing NDVI. In Formula 8 , the correction parameter values increase as the solar altitude increases beyond the assumed standard solar altitude. The effects of weather are also considered in Formula 8 and structured such that the correction factor increases with increasing amounts of direct solar radiation. Next, using Formula 9 , we obtain the variable S for changing the amount of correction depending on the NDVI values for the six plots established within the experimental field. See section 3.3 below for further explanation of the variable S. This is based on regression analysis using the NDVI rate of variability of each observation time relative to the assumed standard time as the objective variable and Cf as the explanatory variable. Next, using Formula 10 , we determine the parameters a and b, based on regression analysis using S as the objective variable and NDVI obs as the explanatory variable. Finally, based on the correction variables sought in formulae 8 -10 , NDVI correction is performed as shown in Formula 11 . In the regression analyses, p-values ≤ 0.01 were considered significant. Figure 2 shows the changes in NDVI field detected by the field sensor on clear and cloudy days. On clear days, we found that NDVI field decreased around 12 : 00, when the solar altitude was high, and thereafter increased. By contrast, on cloudy days NDVI field time changes were smaller. These results are similar to those obtained by Ishihara et al. 2015 and other studies. As the solar altitude increases, more incoming light reaches the soil and water beneath the leaves and stalks, as indicated by the increase in observed light reflected from the ground surface and markedly different NDVI field values. On cloudy days, on the other hand, there is a greater amount of diffuse solar radiation, and subsequently no major change with time of day in the components of light reflected after reaching the ground; this probably accounts for the smaller changes in NDVI field . Figure 3 shows the observation-time series of NDVI field , as measured by the field sensor, from the early growth period to the heading stage. It can be seen that the 10 : 00 NDVI field rises smoothly. The 8 : 00, 12 : 00, and 16 : 00 NDVI field values also rise, but exhibit day-to-day fluctuations. In the correlation analysis between the LOWESS curve and the observed NDVI field shown in Table 1, the coefficient of determination was highest at the 10 : 00 NDVI field . This indicates that the time-series change of the 10 : 00 NDVI field was the smoothest. At the experimental field, the solar altitude at 10 : 00 was between 49.8 and 64.8 . Optical satellite observations are generally between 10 : 00 and 11 : 00. Hashimoto et al. 2019 , however, estimated that the best time for drone RS was when the solar altitude was between 25 and 45 based on a radiative transfer model simulation reasoning that the higher the solar altitude, the greater the observed light reflected from the ground surface. However, in June and July around the time of panicle formation stage an important period for observing growth and accordingly adjusting fertilization regimes solar altitude is between 25 and 45 at approximately 8 : 00. During this time window, rice plants are still coated with morning dew, which can cause errors in the spectral measurements Paul and Pinter, 1986;Maresma et al., Fig. 2. Changes in NDVI detected by a field sensor on clear and cloudy days. 2020 . Paul and Pinter 1986 reported that when morning dew was present on canopies, morning reflectance in wavelengths shorter than 700 nm was significantly different from that observed during the afternoon. In the present study, the morning dew evaporated between 7 : 40 and 9 : 40. Given that this morning dew had generally evaporated by 10 : 00, we reasoned that the appropriate observation time for low-altitude RS, such as drone RS, is 10 : 00, the same as satellite RS. Figure 4 shows the time series of NDVI P4M in the six plots situated within the experimental field, as observed on July 22 nd . Compared to 10 : 00, as the standard observation time, NDVI P4M variability fell by 1 at 12 : 00 in those plots with high NDVI P4M values, whereas there was a 6 decline in plots with low NDVI P4M values. This indicates that the measurements for these latter plots were more strongly affected by solar altitude.   These results are consistent with those obtained by Ishihara et al., 2015 , and indicate that higher NDVI P4M values are obtained at higher vegetation densities because a smaller fraction of incoming light is reflected after reaching the ground surface. We also found that when correcting NDVI P4M values according to sunlight conditions, it was necessary to adjust the amount of correction depending on the NDVI P4M values.

NDVI correction variables
Taking 10 : 00 as an appropriate RS observation time, the assumed standard solar altitude H o of Formula 8 was considered to be 60 , the average 10 : 00 solar altitude during the paddy rice cultivation months, May to September. Figure 5a shows the July 22 nd relationship between Cf the sunlight condition correction variable and the NDVI P4M variability for each observation time, relative to 10 : 00. In all of these relationships, p-values of 0.01 or lower were considered significant. Furthermore, the gradient slope of Cf against time was greater for plots with sparse growth, indicating that the vegetation images obtained for these plots were more notably influenced by the effects of sunlight conditions on NDVI P4M . This slope is the correction variable S in Formula 10 . Figure 5b shows the relationship between the correction variable S on July 22 nd and 10 : 00 NDVI P4M . The p-values were considered significant at 0.01 or lower. Using the regression equation shown in Fig. 5b

Corrections of NDVI observed with P4M
Then, the parameters were applied in Formula 11 , the NDVI P4M was corrected as follows: NDVI std = NDVI obs 1 + 0.014 Cf / 1 + 0.014 NDVI obs Cf , 14 Figure 6 shows the pre-correction time series of NDVI P4M on July 22 nd and August 5 th and post-correction by using Formula 14 . After these corrections, the NDVI P4M changes throughout the day were smaller in all surveyed plots Table 2 . On August 5 th , observations were made at 10 : 30, 11 : 00, 12 : 30, 13 : 00, 14 : 00, 15 : 00, and 16 : 00. There were no data for 12 : 00 the highest solar altitude ; therefore, in the time series of NDVI P4M , the NDVI decline at 12 : 00 was not clear. The NDVI P4M images indicate that the corrected NDVI P4M values account for the effects of sunlight conditions Fig. 7 .
It is generally assumed that the most accurate method for NDVI correction uses the radiative transfer model or the bidirectional reflectance distribution function model; however, these models require numerous complex parameters. For the NDVI P4M corrections of the present study given that solar altitude can be calculated from the location coordinates and the observation times the only necessary on-site observations were the NDVI P4M and global solar radiation. Corrections require only few parameters and, given the general reduction in the sunlight-condition effects on NDVI P4M , this may represent a simple and effective correction method. The results presented herein are based on an empirical formula specifically developed for paddy rice Koshihikari cultivar , and thus would not directly apply without modification to other vegetation types or cultivars. The amount of the correction corresponds to sunlight conditions and variations in growth; however, further study is required to verify the robustness. On the other hand, if the model parameters were to be adjusted for other vegetation types or cultivars, the procedure of NDVI correction used in this study could be applied to other cultivars.

Conclusion
In this study, we sought to develop an empirical correction formula that could be used to correct for sunlight conditions in NDVI images of paddy rice captured by drones.
Based on observations using a field sensor installed in a paddy field, the results of our analysis of the appropriate time for observing NDVI in paddy rice showed that the 10 : 00 NDVI field rose smoothly from the early growth period, whereas there were notable day-to-day fluctuations in the rise of the 8 : 00, 12 : 00, and 16 : 00 NDVI field values. Therefore, we conclude that 10 : 00 is the most appropriate time for RS observations of paddy rice. Although this time differs from the observation timing of a previously described radiative transfer model, we maintain that the effects of morning dew need to be considered; we established that this dew had mainly evaporated by 10 : 00.
Regarding the time series of drone-observed NDVI P4M , we found that the effects of differing sunlight were stronger in paddy plots with a lower vegetation density and a subsequently lower NDVI P4M value. We speculate that this is attributable to the fact that at higher vegetation densities, a smaller fraction of incoming light reaches the ground surface to be subsequently reflected back. We found also that when correcting NDVI P4M values to account for sunlight conditions, it is necessary to adjust the correction parameters depending on the NDVI P4M values.
Based on these considerations, we were able to apply corrections to reduce the variations in NDVI P4M associated with sunlight conditions. When using the correction method described herein, the only necessary on-site observations were NDVI P4M and global solar radiation, thereby making this a simple correction method with very few parameters. Despite its simplicity, it is effective in mitigating the influence of sunlight conditions on NDVI P4M , and thus has considerable practical utility. However, it should be emphasized that the results obtained in the present study are based on an empirical formula specifically developed for paddy rice; accordingly, caution should be exercised when applying the specified correction parameters to other types of vegetation or cultivars.

Acknowledgments
This work was supported by JSPS KAKENHI Grant Number 19J00437, 20H01388 , and was conducted as a joint research program of CEReS, Chiba University 2020 .