Estimating the Protein Concentration in Rice Grain Using UAV Imagery Together with Agroclimatic Data

Global warming and climate change can potentially change not only rice production but also rice nutrient content. To adapt a rice-dependent country’s farming to the impacts of climate change, it is necessary to assess and monitor the potential risk that climate change poses to agriculture. The aim of this study was to clarify the relationship between rice grain protein content (GPC) and meteorological variables through unmanned aerial vehicle remote sensing and meteorological measurements. Furthermore, a method for GPC estimation that combines remote sensing data and meteorological variables was proposed. The conclusions of this study were as follows: (1) The accuracy and robustness of the GPC estimation model were improved by evaluating the nitrogen condition with the green normalized difference vegetation index at the heading stage (GNDVIheading) and evaluating photosynthesis with the average daily solar radiation during the grain-filling stage (SRgrain-filling). GPC estimation considering SRgrain-filling in addition to GNDVIheading was able to estimate the observed GPC under the different conditions. (2) Increased temperature from the transplantation date to the heading stage can affect increased GPC when extreme temperature does not cause the heat stress.


Introduction
The United Nation's 17 Sustainable Development Goals (SDGs) came into effect in January 2016, and they will continue until 2030. Many SDGs are related to agriculture (e.g., SDG2 Zero hunger: End hunger, achieve food security and improved nutrition and promote sustainable agriculture) that require urgent action from both developing and developed countries [1]. However, implications of climate change, such as global warming, have the potential to interrupt progress being made toward SDGs related to agriculture [2]. With respect to rice production in monsoon Asia, where rapid population growth is driving increased food demand, an increase in temperature and fluctuating precipitation were expected to reduce rice yield, whereas a rise in carbon dioxide (CO 2 ) concentration was expected to increase rice yield because of promoted photosynthesis [3][4][5][6]. To adapt agriculture to the impacts of climate change, it is necessary to assess and monitor the potential risk that climate change poses to agriculture. Remote sensing and imagery techniques are useful in monitoring and detecting the effects of climate change on crops [3,[7][8][9][10]; however, implementation of remote sensing and imagery techniques for monitoring and detecting the effects of climate change is still underdeveloped [11]. Furthermore, while much study has been done on changes in crop production, study on changes in crop nutrient content has been limited [12][13][14].      The study observation equipment included an electric-powered Multicopter (Zion QC630, enRoute) and a multispectral camera (Yubaflex, BIZWORKS). The Yubaflex is a modified version of the Canon PowerShot S100 camera, which takes images at the green, red, and near-infrared (NIR) bands. The bandwidth of each band was as follows: Green 500-600nm, Red 600-850nm, and NIR 700-1050nm (wavelengths showing the maximum spectral response of each band: Green 550 nm, Red 600 nm, and NIR 850 nm). The image was made up of 12 million pixels (4000 × 3000). The camera can also convert the observed digital number to radiance using the dedicated software Yubaflex 3.1 [33]. In this study, we used the Yubaflex for vegetation index monitoring and GPC estimation.
The UAV-based observations were acquired once a week between 10:00 and 10:30 a.m. local time, under both clear and cloudy sky conditions. The flight altitude was 50 m above ground level (ground resolution: 1.8 cm), and the overlap of each image was 70%. The settings of the Yubaflex were shutter speed priority mode, shutter speed set to 1/1000 sec, ISO sensitivity set to Auto, and interval of shooting set to the minimum value (approximately 4 sec).

Image Processing and Analysis
The images taken with the Yubaflex, after being converted to radiance using the dedicated software Yubaflex3.1, were used to create orthomosaic images using Structure from Motion-Multi View Stereo (SfM-MVS) software (Agisoft PhotoScan Professional v1.4.1).
Yubaflex-based NDVI values are relatively lower than the other multispectral camera-based NDVI [34]. The bandwidth of Yubaflex red band is 600-850nm, and some of NIR wavelength are included. Therefore, Yubaflex red band continue to respond to high aboveground biomass, and NDVI values become relatively low. In this study, we used the green NDVI (GNDVI) [35] in order not to become vegetation index values low. GNDVI was calculated based on equation (1): where GNDVI is the Yubaflex-based GNDVI, and NIR Yubaflex and Green Yubaflex are the Yubaflex-based NIR band and Green band radiances, respectively. The average GNDVI was then calculated for each plot ( Figure 1) using ESRI ArcGIS 10.4. In general, the vegetation indices decreased with decreasing solar zenith angle. This response was affected significantly by the growth stage and diffuse/direct light conditions [36]. Ishihara et al. (2015) [36] reported that the decreasing response of the vegetation indices to the decreasing solar zenith angle was high during the middle growth stage and low at the heading stage. In addition, the response of vegetation indices to the solar zenith angle was evident under clear sky conditions at large solar zenith angles (less than 20 • ) but almost negligible under cloudy sky conditions [36]. In this study, UAV-based observations were acquired between 10:00 and 10:30 a.m. local time. The solar zenith angle at the observation was more than 20 • . The sunlight conditions appear to have a low effect on the consistency of GNDVI in this study.

Analysis of Collected Samples and Meteorological Data
At the Chiba site, GPC (converted to moisture content 15%) was observed from 2015 to 2017. Samples in which lodging occurred were excluded, and the analysis covered 24 samples from 2015, 23 samples from 2016, and 21 samples from 2017. Total nitrogen of grain adjusted with a grain thickness of 1.8 mm was observed using the NC analyzer (Sumigraph NC-900, Sumika Analysis Service [37]). Subsequently, the conversion factor (5.95) of nitrogen-protein was multiplied to covert total nitrogen to GPC.
In this study, the daily mean solar radiation and the daily mean temperature from agricultural weather data obtained from 1 km grid squares [38] were analyzed for GPC estimation. These agro-meteorological grid square data are based on data from the ground observation station, known as the Automated Meteorological Data Acquisition System network.

GPC Estimation
Two types models (simple linear regression model and multiple regression model) for the estimation of GPC were developed; simple linear regression (SLR) model was built by using GNDVI at the heading stage (GNDVI heading ) as explanatory variable, while multiple regression (MR) model was built by using GNDVI heading and average daily solar radiation during the grain-filling stage (SR grain-filling ) as explanatory variables. Figure 2 denotes the GPC estimation process. Service [37]). Subsequently, the conversion factor (5.95) of nitrogen-protein was multiplied to covert total nitrogen to GPC. In this study, the daily mean solar radiation and the daily mean temperature from agricultural weather data obtained from 1 km grid squares [38] were analyzed for GPC estimation. These agrometeorological grid square data are based on data from the ground observation station, known as the Automated Meteorological Data Acquisition System network.

GPC Estimation
Two types models (simple linear regression model and multiple regression model) for the estimation of GPC were developed; simple linear regression (SLR) model was built by using GNDVI at the heading stage (GNDVIheading) as explanatory variable, while multiple regression (MR) model was built by using GNDVIheading and average daily solar radiation during the grain-filling stage (SRgrain-filling) as explanatory variables. Figure 2 denotes the GPC estimation process. The samples collected from 2015 to 2016 were used to build the empirical GPC estimation models for three cultivars (Koshihikari, Fusaotome, and Fusakogane). At first, a SLR analysis for three cultivars was conducted to explore the relationship between GPC and GNDVIheading. Following the SLR analysis, the GPC estimation models were obtained as: where GPC is the grain protein content (%), GNDVIheading is the Yubaflex-based GNDVI at the heading stage, m is the coefficient, and k is the intercept. Then, a MR analysis for three cultivars was conducted to explore the relationship between GPC, GNDVIheading, and SRgrain-filling. Following the multiple regression analysis, the GPC estimation models were obtained as: where GPC is the grain protein content (%), GNDVIheading is the Yubaflex-based GNDVI at the heading stage, SRgrain-filling is the average daily solar radiation during the grain-filling stage (MJ/m²), mi are the coefficients, and k is the intercept. GNDVIheading was used for evaluation of canopy nitrogen content, and SRgrain-filling was used for evaluation of photosynthesis. The results of both single and multiple regression analysis were evaluated by the coefficient of determination (R²), the P value. A p value lower than 0.05 was considered statistically significant. Finally, GPC estimation model robustness was validated by using the observed GPC collected in 2017. The validation results were evaluated in terms of the root mean square error (RMSE). A sensitivity analysis of the MR model was performed in terms of the rate of change in the output value resulting from a change of each input parameter while keeping all other parameters constant The samples collected from 2015 to 2016 were used to build the empirical GPC estimation models for three cultivars (Koshihikari, Fusaotome, and Fusakogane). At first, a SLR analysis for three cultivars was conducted to explore the relationship between GPC and GNDVI heading . Following the SLR analysis, the GPC estimation models were obtained as: where GPC is the grain protein content (%), GNDVI heading is the Yubaflex-based GNDVI at the heading stage, m is the coefficient, and k is the intercept. Then, a MR analysis for three cultivars was conducted to explore the relationship between GPC, GNDVI heading , and SR grain-filling . Following the multiple regression analysis, the GPC estimation models were obtained as: where GPC is the grain protein content (%), GNDVI heading is the Yubaflex-based GNDVI at the heading stage, SR grain-filling is the average daily solar radiation during the grain-filling stage (MJ/m 2 ), m i are the coefficients, and k is the intercept. GNDVI heading was used for evaluation of canopy nitrogen content, and SR grain-filling was used for evaluation of photosynthesis. The results of both single and multiple regression analysis were evaluated by the coefficient of determination (R 2 ), the p value. A p value lower than 0.05 was considered statistically significant. Finally, GPC estimation model robustness was validated by using the observed GPC collected in 2017. The validation results were evaluated in terms of the root mean square error (RMSE). A sensitivity analysis of the MR model was performed in terms of the rate of change in the output value resulting from a change of each input parameter while keeping all other parameters constant [32]. The maximum, minimum, and average values collected in the experiment from 2015 to 2017 were used for sensitivity analysis.  used temperature instead of SR grain-filling for GPC estimation [30]. In the case of Koshihikari, the average temperature from 0 to 20 days after the heading stage was used for the GPC estimation [30]. In the cases of Fusakogane and Fusaotome, the average temperatures from 0 to 30 days after heading stage were used for the GPC estimation. In this study, the determination of the SR grain-filling duration was based on   [30]. In the case of Koshihikari, the average daily solar radiation from 0 to 20 days after the heading stage was used as SR grain-filling . In the cases of Fusakogane and Fusaotome, the average daily solar radiation from 0 to 30 days after the heading stage was used as SR grain-filling . Table 2 denotes the results of the SLR analysis. The regression coefficient (m), intercept (k), R 2 , and p value are listed in Table 2. The results of SLR analysis demonstrate that the higher the GNDVI heading , the higher the GPC. The p values of GNDVI heading were lower than 0.05, and these were considered statistically significant.  Table 3 denotes the results of the MR analysis. The regression coefficients (m i ), intercept (k), R 2 , and p value are listed in Table 3. The results of MR analysis demonstrate that the higher the GNDVI heading , the higher the GPC, and the higher the SR grain-filling , the lower the GPC. Furthermore, with respect to the coefficients and intercepts, Fusaotome and Fusakogane were similar, but Koshihikari differed from the other two cultivars. The p values of GNDVI heading and SR grain-filling were lower than 0.05, and these were considered statistically significant.  Table 4 shows the results of sensitivity analysis of MR models. The GPC variation (%) in GNDVI heading of Koshihikari, Fusaotome, and Fusakogane were -7.7 to +11.2, -16.4 to +10.2, and -12.7 to +11.8, respectively. The GPC variation (%) in SR grain-filling of Koshihikari, Fusaotome, and Fusakogane were -6.3 to +6.6, -6.7 to +7.9, and -6.8 to +8.7, respectively. In all cultivars, sensitivity analysis revealed a higher GNDVI heading sensitivity to GPC. In other words, GNDVI heading had a greater effect on GPC.  The comparison between the estimated GPC and the observed GPC is shown in Figure 3. With respect to SLR-based estimation, the RMSE of Koshihikari was 0.61 (average observed GPC 7.08%), the RMSE of Fusaotome was 0.67 (average observed GPC 7.36%), and the RMSE of  The comparison between the estimated GPC and the observed GPC is shown in Figure 3. With respect to SLR-based estimation, the RMSE of Koshihikari was 0.61 (average observed GPC 7.08%), the RMSE of Fusaotome was 0.67 (average observed GPC 7.36%), and the RMSE of Fusakogane was 0.58 (average observed GPC 7.47%). With respect to MR-based estimation, the RMSE of Koshihikari, Fusaotome, and Fusakogane were 0.42, 0.34, and 0.33, respectively. The results of MR-based estimation outperformed the SLR-based estimation.   Figure 4 shows the GNDVI time series in Koshihikari, Fusaotome, and Fusakogane transplanted on May 13 2016. In all cultivars, GNDVI increased toward the heading stage and declined gradually after the heading stage. The peak of GNDVI in Koshihikari was 82 days after transplantation, whereas it was 68 days after transplantation for Fusaotome and Fusakogane. The time series patterns of the allied cultivars (Fusaotome and Fusakogane) were similar, and the peak of GNDVI was 14 days earlier than with Koshihikari. GNDVI time series clearly showed change speed with development depending on the cultivar. In addition, the peak of GNDVI was recorded at the heading stage in all cultivars.  with Koshihikari. GNDVI time series clearly showed change speed with development depending on the cultivar. In addition, the peak of GNDVI was recorded at the heading stage in all cultivars.       Koshihikari in 2017, and these were considered statistically significant. Although the cultivar was the same, the correlation between the average temperature from the transplantation date to the heading stage and the GNDVI heading was varied with respect to collected year.   Figure 7 shows the spatial and temporal variability in GNDVI. The value of GNDVI varied in plots with different cultivation conditions (transplantation date, cultivar, amount of fertilizer). Even if the transplantation date and cultivar were the same, the value of GNDVI increased as the amount of fertilizer increased.

Discussion
The regression coefficient both SLR and MR analysis showed that there was an increase in GPC with increasing GNDVIheading. The P values of GNDVIheading were significant at 5% or less in all cultivars. In the two types (SLR and MR) models for GPC estimation, GNDVIheading was used for evaluation of canopy nitrogen content. The relation between GNDVIheading and canopy nitrogen content is supported by the findings of Inoue et al. (2012) who examined simple and robust methods for remote sensing of canopy nitrogen content [39].
In the MR models for GPC estimation, SRgrain-filling was used for evaluation of photosynthesis. There was a decline in GPC with increasing SRgrain-filling. In addition, the P values of SRgrain-filling were significant at 5% or less in all cultivars. These findings were consistent with those of a previous study in which GPC was shown to decrease because of the increased carbohydrate production via

Discussion
The regression coefficient both SLR and MR analysis showed that there was an increase in GPC with increasing GNDVI heading . The p values of GNDVI heading were significant at 5% or less in all cultivars. In the two types (SLR and MR) models for GPC estimation, GNDVI heading was used for evaluation of canopy nitrogen content. The relation between GNDVI heading and canopy nitrogen content is supported by the findings of Inoue et al. (2012) who examined simple and robust methods for remote sensing of canopy nitrogen content [39].
In the MR models for GPC estimation, SR grain-filling was used for evaluation of photosynthesis. There was a decline in GPC with increasing SR grain-filling . In addition, the p values of SR grain-filling were significant at 5% or less in all cultivars. These findings were consistent with those of a previous study in which GPC was shown to decrease because of the increased carbohydrate production via photosynthesis in rice grains [13,14]. In addition, the R 2 and RMSE of the MR models including the SR grain-filling outperformed the previous study's method (SLR models) for GPC estimation, in which only the correlation between the GNDVI heading and the observed GPC was used (Figure 3). GPC estimation that considers only SR grain-filling in addition to GNDVI heading could improve estimation accuracy and improve the robustness of the estimation model. Furthermore, the observed GPC under the different conditions could be estimated with high accuracy without remaking the regression models for GPC estimation for each transplantation date. This was also an improvement on the previous studies of GPC estimation using remote sensing.
Sensitivity analysis revealed that GNDVI heading had a greater effect on GPC (Table 4). As previously stated, there is a correlation between GNDVI and canopy nitrogen content [39], these findings were consistent with those of a previous study in which nitrogen fertilizer strongly affected GPC [18,23].
It is worth noting that the parameters of the GPC estimation model are specific to the cultivar. The characteristics of each cultivar, such as growth speed and the nutrient translocation from source to sink, are different. According to previous studies, the effect of cultivar difference is as great as that of fertilization [19]. It is necessary to obtain parameters specific to the cultivar based on MR analysis using GNDVI heading and SR grain-filling ; therefore, it is necessary to perform tests with multiple transplantation dates, as done in this study, to optimize the parameters of the GPC estimation model. However, the parameters of the GPC estimation model might be similar in cases of allied cultivars, such as the Fusaotome and Fusakogane cultivars in this study.
MR models could estimate the observed GPC under the different conditions as described above; however, MR models would underestimate the GPC in case topdressing after the heading stage. The timing of topdressing affects protein, and topdressing after the heading stage greatly increases GPC [40][41][42]. MR models use GNDVI at the heading stage; therefore, the effect of topdressing after the heading stage cannot be considered.
As shown in Figure 4, the peak of GNDVI was recorded at the heading stage in all cultivars. This finding was consistent with those of a previous study in which monitoring the paddy field using remote sensing [43][44][45]. GNDVI time series (Figure 4) showed change speed with development depending on the cultivar.
The days after transplantation at the peak of GNDVI decreased when the transplantation date was later ( Figure 5). Changes in GNDVI time series due to differences in transplantation date are in agreement with the findings of the previous studies that the development speed of paddy rice becomes faster when the temperature increased during the growing season [46,47].
As shown in Figure 6, GNDVI heading increased when the average temperature from the transplantation date to the heading stage increased. In addition, the maximum value of GNDVI increased when the transplantation date was later ( Figure 5). Nitrogen absorbed by plants is decomposed into inorganic nitrogen in the soil. Previous study has demonstrated that the amount of inorganic nitrogen increases as the temperature rises [48]. Therefore, this study might show that as the average temperature from the transplantation date to the heading stage increased, the amount of inorganic nitrogen and nitrogen absorbed by the plants increased, even when the amount of fertilizer used was unchanged. Consequently, GNDVI heading increased as the transplantation date became later.
There are two different hypotheses concerning the effect of temperature on GPC: some studies suggest that increasing temperature can reduce GPC [15,16], whereas others suggest that increasing temperature can increase GPC [17,18]. As described above, the findings of the current study show that increasing temperature from the transplantation date to the heading stage can increase GPC. However, the findings of current study only follow the increase of temperature from spring to summer. Furthermore, the extreme temperature or temperature well above the average during crop development could cause the heat stress and crops can show rapid development without GPC increasing [32]. Therefore, the findings of current study suggested that increasing temperature from the transplantation date to the heading stage can increase GPC when extreme temperature does not cause the heat stress.

Conclusions
We used UAV-RS data and meteorological measurements to clarify the relationship between GPC and meteorological variables. Furthermore, a method for GPC estimation that combines remote sensing data and meteorological variables has been proposed.
We proposed the simple method for GPC estimation by using GNDVI heading together with SR grain-filling . GNDVI heading was used for evaluation of canopy nitrogen content, and SR grain-filling was used for evaluation of photosynthesis. MR analysis and the GPC estimation models for three cultivars showed that the higher the GNDVI heading , the higher the GPC, and the higher the SR grain-filling , the lower the GPC. Additionally, the p value of GNDVI heading and SR grain-filling were significant at 5% or less in all cultivars. The validation of GPC estimation showed that the RMSE of Koshihikari was 0.42 (average observed GPC 7.08%), the RMSE of Fusaotome was 0.34 (average observed GPC 7.36%), and the RMSE of Fusakogane was 0.33 (average observed GPC 7.47%). Although the parameters of the GPC estimation model are specific to the cultivar, it was possible to improve GPC estimation accuracy and model robustness. Estimation of GPC that only considers SR grain-filling in addition to GNDVI heading could estimate the observed GPC under the different conditions without remaking the regression models for each transplanting date.
The GNDVI time series and the correlation between the average temperature from the transplantation date to the heading stage and the GNDVI heading showed that GNDVI heading increased when the temperature from the transplantation date to the heading stage increased although the amount of fertilizer was the same. In addition, MR models for GPC estimation showed that GPC increased when GNDVI heading increased. Therefore, this study has shown that increasing temperature from the transplantation date to the heading stage can affect increased GPC when extreme temperature does not cause the heat stress.