Enhancing the Nitrogen Signals of Rice Canopies across Critical Growth Stages through the Integration of Textural and Spectral Information from Unmanned Aerial Vehicle (UAV) Multispectral Imagery

: This paper evaluates the potential of integrating textural and spectral information from unmanned aerial vehicle (UAV)-based multispectral imagery for improving the quantiﬁcation of nitrogen (N) status in rice crops. Vegetation indices (VIs), normalized di ﬀ erence texture indices (NDTIs), and their combination were used to estimate four N nutrition parameters leaf nitrogen concentration (LNC), leaf nitrogen accumulation (LNA), plant nitrogen concentration (PNC), and plant nitrogen accumulation (PNA). Results demonstrated that the normalized di ﬀ erence red-edge index (NDRE) performed best in estimating the N nutrition parameters among all the VI candidates. The optimal texture indices had comparable performance in N nutrition parameters estimation as compared to NDRE. Signiﬁcant improvement for all N nutrition parameters could be obtained by integrating VIs with NDTIs using multiple linear regression. While tested across years and growth stages, the multivariate models also exhibited satisfactory estimation accuracy. For texture analysis, texture metrics calculated in the direction D3 (perpendicular to the row orientation) are recommended for monitoring row-planted crops. These ﬁndings indicate that the addition of textural information derived from UAV multispectral imagery could reduce the e ﬀ ects of background materials and saturation and enhance the N signals of rice canopies for the entire season.


Introduction
As one of the most important staple crops around the world, rice (Oryza sativa L.) crop feeds more than 50% of the world's population. Nitrogen is the primary element for crop growth due to the crucial impact on crop yield formulation and grain quality determination. Monitoring N status in rice leaves or plants could provide valuable information for growth diagnosis [1] and precise field N management [2], so as to improve the N use efficiency and reduce environmental pollution.

Ground Sampling and UAV Image Acquisition
The UAV campaigns and ground destructive samplings were taken at critical growth stages of rice (Table 1). After the UAV campaign, three hills of rice plants were randomly harvested in each plot and separated into different organs. All the samples were oven dried at 105 • C for 30 mins, followed by drying at 80 • C, until a constant weight was reached. Then, they were weighed, ground, and stored in plastic bags for laboratory chemical analysis. Aboveground biomass (AGB) was the total dry weight of all plant organs per unit ground area. Total N concentrations (%) were determined with the micro-Keldjahl analysis. The LNC and PNC were expressed on a dry weight basis (%), and LNA and PNA were expressed as N mass per unit ground area (g m −2 ). The LNA (g m −2 ) was obtained through multiplying LNC (%) by leaf dry biomass (LDB, t ha −1 ). The PNA (g m 2 ) was obtained through multiplying PNC (%) by AGB (t ha −1 ).  In this study, we used an eight-rotor aircraft (Mikrokopter OktoXL) which had a maximum payload capacity of 2.5 kg and a flight duration of 10−20 min, depending on the battery and actual payload. A multispectral camera (Mini-MCA6; Tetracam, Inc., Chatsworth, CA, USA) was mounted on the UAV to acquire images in this study ( Table 1). The specification of this multispectral camera can be found in Reference [17]. The UAV was flown at 100 m above ground level, resulting in a nominal resolution of 0.054 m in the multispectral images. The camera was set to continuous data capture at one frame per three seconds with a fixed aperture and exposure according to the light conditions. Aerial images were saved in the camera memory card in a 10 bit RAW format.
The UAV campaigns were conducted under a clear sky and calm wind conditions between 11:00 and 13:30 (local time). After each flight, only one image that covered all 36 plots in the nadir position was selected for subsequent analysis.

Image Processing
The UAV-based image preprocessing workflows were performed by following Reference [17] and processed in IDL/ENVI environment (Exelis Visual Information Solutions, Boulder, Colorado, USA). The pre-processing workflows included noise reduction, vignetting correction, and lens distortion correction. Subsequently, the six bands were co-registered with 25 ground control points (GCPs) evenly Remote Sens. 2020, 12, 957 4 of 17 distributed in the study area ( Figure 1). Then, we conducted radiometric correction with an empirical line correction method using a six flat calibration canvas at different reflectance intensities [17,18].
The grey level co-occurrence matrix (GLCM) was employed in this study to evaluate the potential of texture analysis on reflectance images for improving rice N status monitoring. Eight GLCM-based texture metrics including mean (MEA), variance (VAR), homogeneity (HOM), contrast (CON), dissimilarity (DIS), entropy (ENT), second moment (SEM), and correlation (COR) were computed using the ENVI software. A detailed description of the eight texture metrics can be found in Reference [19]. Since the smallest window size (3 × 3) was comparable to the row spacing at jointing stage and larger window sizes did not exhibit significant differences (data not shown), only the window size (3 × 3) was investigated. Next, inter-correlations between different texture metrics were analyzed to decrease the data dimensionality and improve the data processing efficiency. Finally, we chose four texture metrics (i.e., MEA, CON, DIS, and COR) at the minimal window size (3 × 3) for texture analysis with five bands (blue, green, red, RE, and NIR) in four directions ( Based on individual texture bands, the normalized difference texture index (NDTI) was also derived to improve the performance of texture analysis. The NDTI was calculated in the equation below: where T 1 and T 2 are a random texture feature from the five bands in four directions.

Statistical Analysis
The data collected from the two field experiments were pooled to examine the relationships of agronomic variables with VIs, NDTIs, and their combinations by simple regression (SR) and stepwise multiple linear regression (SMLR). Model validation was performed on the pooled data using a k-fold Remote Sens. 2020, 12, 957 5 of 17 (k = 10) cross-validation procedure and evaluated by the root mean square error (RMSE) and the relative RMSE (RRMSE). Table 2. List of vegetation indices (VIs) used in this study.

Index Formula Reference
Color indices Normalized green red difference index (NGRDI) NGRDI = (R 550 − R 680 )/(R 550 + R 680 ) [20] Visible atmospherically resistance index (VARI) Green leaf index (GLI) RE indices MERIS terrestrial chlorophyll index MTCI = (R 800 − R 720 )/(R 720 + R 680 ) [23] red edge chlorophyll index (CI RE ) CI RE = R 800 /R 720 − 1 [24] Normalized difference red edge index (NDRE) NDRE = (R 800 − R 720 )/(R 800 + R 720 ) [25] NIR indices Normalized difference vegetation index (NDVI) NDVI = (R 800 − R 680 )/(R 800 + R 680 ) [26] Renormalized difference vegetation index (RDVI) RDVI = (R 800 − R 680 )/ (R 800 + R 680 ) [27] Optimized soil adjusted vegetation index (OSAVI) Modified triangular vegetation index 2 (MTVI2) The bands used in the formula corresponded to the bands in the multispectral imagery. Table 3 shows the relationships among different agronomic parameters. The LNC and PNC and LNA and PNA were highly correlated to each other. In addition, biomass was highly correlated to both LNA and PNA. Table 4 shows the relationships between N nutrition parameters and VIs derived from UAV multispectral images. For foliar parameters, the relationships between VIs and LNC were remarkably low, and the highest R 2 was only 0.2 obtained by CI RE and NDRE ( Figure 1a). These two indices also exhibited the strongest relationships with LNA (R 2 = 0.77 for both) (Figure 1b). For plant parameters, PNC was most correlated to CI RE as seen for LNC but with a lower R 2 value (R 2 = 0.15) ( Figure 1c). The OSAVI exhibited the strongest relationships with PNA (R 2 = 0.73) which was also found for NDRE (R 2 = 0.73) (Figure 1d). In general, RE-based VIs exhibited better relationships with N nutrition parameters than other indices.   The numbers in bold correspond to the best-fit function for each N nutrition parameter. Significance level: "ns" represents no significance, * p < 0.01, ** p < 0.001, *** p < 0.0001. "L" and "E" indicate linear and exponential regression. Figure 2 shows the relationships between N nutrition parameters and single-band texture metrics. For LNC and PNC estimation, the texture metrics in the direction D3 were generally superior to those in other directions, while directions did not influence the performance of the optimal texture feature for LNA and PNA estimation. The DIS feature and the MEA feature exhibited a higher correlation than other texture metrics with N concentration and N accumulation indicators, respectively. The DIS feature in the RE band was superior to other bands in LNC estimation with an R 2 of only 0.26 ( Figure 2a). The MEA feature in the red band was most related to LNA (R 2 = 0.39) among all the texture metrics ( Figure 2b). Similar to LNC, the DIS feature in the RE band exhibited a higher correlation with PNC compared to their counterparts ( Figure 2c). The MEA feature in the NIR band obtained a higher correlation with PNA (R 2 = 0.55) than other texture metrics ( Figure 2d).

Relationships between N Nutrition Parameters and Texture Indices
Considerable improvement was obtained with texture indices for N nutrition parameters estimation comparing with individual texture metrics ( Table 5). The majority of optimal NDTIs were composed of texture metrics in D3 and D4. The optimal NDTIs for LNC and PNC estimation consisted of DIS, CON, and COR features calculated from RE and NIR bands. While the optimal NDTIs for LNA estimation consisted of MEA features calculated from RE and NIR bands, and those for PNA estimation were composed of MEA features calculated from green and NIR bands. The index (NDTI(DIS B5_D4 , CON B4_D3 )) was closely related to both LNC and PNC with R 2 of 0.29 and 0.41 (Figure 3a,c). The NDTI(MEA B5_D3 , MEA B4_D4 ) performed best in the LNA estimation (R 2 = 0.72) with exponential regression (Figure 3b). The NDTI(MEA B5_D3 , MEA B2_D3 ) performed best in the PNA estimation (R 2 = 0.75), and obvious saturation occurred at high PNA levels ( Figure 3d). feature for LNA and PNA estimation. The DIS feature and the MEA feature exhibited a higher correlation than other texture metrics with N concentration and N accumulation indicators, respectively. The DIS feature in the RE band was superior to other bands in LNC estimation with an R 2 of only 0.26 ( Figure 2a). The MEA feature in the red band was most related to LNA (R 2 = 0.39) among all the texture metrics ( Figure 2b). Similar to LNC, the DIS feature in the RE band exhibited a higher correlation with PNC compared to their counterparts ( Figure 2c). The MEA feature in the NIR band obtained a higher correlation with PNA (R 2 = 0.55) than other texture metrics ( Figure 2d).

Relationships between N Nutrition Parameters and Texture Indices
Considerable improvement was obtained with texture indices for N nutrition parameters estimation comparing with individual texture metrics ( Table 4). The majority of optimal NDTIs were composed of texture metrics in D3 and D4. The optimal NDTIs for LNC and PNC estimation consisted of DIS, CON, and COR features calculated from RE and NIR bands. While the optimal NDTIs for LNA estimation consisted of MEA features calculated from RE and NIR bands, and those for PNA estimation were composed of MEA features calculated from green and NIR bands. The index (NDTI(DISB5_D4, CONB4_D3)) was closely related to both LNC and PNC with R 2 of 0.29 and 0.41 ( Figure  3a,c). The NDTI(MEAB5_D3, MEAB4_D4) performed best in the LNA estimation (R 2 = 0.72) with exponential regression (Figure 3b). The NDTI(MEAB5_D3, MEAB2_D3) performed best in the PNA estimation (R 2 = 0.75), and obvious saturation occurred at high PNA levels ( Figure 3d). Table 4. Relationships (R 2 ) between N nutrition parameters and the top four best-performing texture indices.   Table 6 shows the performance of the combination of VIs and texture indices in the estimation of N nutrition parameters with stepwise multiple linear regression. The hybrid model with VIs and NDTIs could explain 70% of the variability in LNC which was significantly higher than using VIs or texture indices alone. Compared with the optimal VI (NDRE) for LNA estimation, the combination of VIs and NDTIs produced a slightly better relationship. Similar to the LNC estimation, the multivariate model improved the regression significantly from 0.41 to 0.68 with three texture indices. Furthermore, the hybrid model explained more than 80% of the variability in PNA (R 2 = 0.86) which was higher than using individual remote sensing variables alone. Interestingly, all the multivariate models for plant-level variables included NDTI(MEA B5_D3 , MEA B2_D4 ). Furthermore, all variables in the hybrid models were statistically significant and multicollinearity effects were not observed.   Table 5 shows the performance of the combination of VIs and texture indices in the estimation of N nutrition parameters with stepwise multiple linear regression. The hybrid model with VIs and NDTIs could explain 70% of the variability in LNC which was significantly higher than using VIs or texture indices alone. Compared with the optimal VI (NDRE) for LNA estimation, the combination of VIs and NDTIs produced a slightly better relationship. Similar to the LNC estimation, the multivariate model improved the regression significantly from 0.41 to 0.68 with three texture indices. Furthermore, the hybrid model explained more than 80% of the variability in PNA (R 2 = 0.86) which was higher than using individual remote sensing variables alone. Interestingly, all the multivariate models for plant-level variables included NDTI(MEAB5_D3, MEAB2_D4). Furthermore, all variables in the hybrid models were statistically significant and multicollinearity effects were not observed.    Table 7 shows the accuracy assessment for N nutrition parameters with the best-fit model in each category. The multivariate model yielded the significantly higher estimation accuracy (RMSE = 0.22 and RRMSE = 9.70%) for LNC estimation than using VIs and NDTIs alone (Figure 4a-c). Compared with the optimal VI, the combination of VI and NDTI exhibited a slightly higher estimation accuracy for LNA (RMSE = 1.16 g m −2 and RRMSE = 21.04%) estimation. The highest estimation accuracy (RMSE = 0.19 and RRMSE = 15.34%) for PNC estimation was obtained by the hybrid model, which was significantly higher than other techniques (Figure 4d-f). Similarly, the multivariate models yielded the highest estimation accuracy for PNA (RMSE = 2.38 g m −2 and RRMSE = 24.17%) estimation when compared with other techniques. Furthermore, the optimal multivariate models for LNC and PNC in Table 6 exhibited satisfactory performance when validated with different dataset groups (Table 8). each category. The multivariate model yielded the significantly higher estimation accuracy (RMSE = 0.22 and RRMSE = 9.70%) for LNC estimation than using VIs and NDTIs alone (Figure 4a-c). Compared with the optimal VI, the combination of VI and NDTI exhibited a slightly higher estimation accuracy for LNA (RMSE = 1.16 g m −2 and RRMSE = 21.04%) estimation. The highest estimation accuracy (RMSE = 0.19 and RRMSE = 15.34%) for PNC estimation was obtained by the hybrid model, which was significantly higher than other techniques (Figure 4d-f). Similarly, the multivariate models yielded the highest estimation accuracy for PNA (RMSE = 2.38 g m −2 and RRMSE = 24.17%) estimation when compared with other techniques. Furthermore, the optimal multivariate models for LNC and PNC in Table 5 exhibited satisfactory performance when validated with  different dataset groups (Table 7). When all the texture indices in the multivariate models used the texture metrics calculated in the direction D3 for simplicity, the validation accuracies for the majority of N nutrition parameters were marginally lower (data not shown). For example, the validation RMSE and RRMSE values of the LNC estimation decreased slightly to 0.24% and 10.52%, respectively.   When all the texture indices in the multivariate models used the texture metrics calculated in the direction D3 for simplicity, the validation accuracies for the majority of N nutrition parameters were marginally lower (data not shown). For example, the validation RMSE and RRMSE values of the LNC estimation decreased slightly to 0.24% and 10.52%, respectively.

Differences among N Nutrition Parameters: N Concentration versus N Accumulation
The estimation accuracies for LNC and PNC were low when using VIs and the best-performing VI could explain no more than 30% of the variability in N concentration. However, the estimation accuracies for LNA and PNA were significantly higher than those for LNC and PNC. The optimal VI (NDRE) could explain 77% and 73% of the variability in LNA and PNA, respectively. The LNC (or PNC) varied in a narrow range across the growing season of rice crops. It decreased from the beginning, leveled off in the middle and then decreased until harvesting [6]. This trend could be characterized by NDRE with a linear but weak relationship. However, LNA (or PNA) varied in a relatively wider range and kept increasing across the entire season. Particularly, NDRE tended to saturate when LNA (PNA) increased to a certain level. A linear regression would not be able to capture this decreasing rate of nitrogen accumulation. Therefore, LNC or PNC exhibited weak relationships with VIs in linear regressions and LNA or PNA in exponential regression [3,5]. Besides, the N absorption features are located in the shortwave infrared (SWIR) region rather than the visible and NIR region [30]. Li et al. [31] improved crop LNC estimation from the SWIR reflectance spectra of fresh leaves through enhancing the absorption features of nitrogen in the SWIR region. However, those VIs derived from UAV multispectral imagery were all based on visible and NIR bands, which could be easily affected by chlorophyll content and canopy structure [32]. Furthermore, the N dilution effect might be another reason for the low estimation accuracy of LNC and PNC with VIs due to the decrease of N concentration along with the increase of biomass [33]. The good relationships between NA and VIs might be attributed to the strong capabilities of VIs in retrieving biomass [34] and the close correlation between N accumulation and biomass (Table 3).
In texture analysis, the optimal NDTI for LNC and PNC estimation were all composed of NIR and RE bands with similar texture metrics (e.g., CON, DIS, COR). Although both LNC and PNC could be estimated with a common index NDTI (DIS B4_D3 , COR B5_D3 ), the relationships between texture indices and NC were not strong in both linear and exponential regressions. That might be because leaf color turned light to dark repeatedly in response to N fertilization at the vegetative stages and kept yellowing at reproductive stages [6]. This reduced the sensitivity of texture metrics to the heterogeneity of tonal variation caused by N status in crops across the whole growing season. Specifically, NDTI(MEA B5 , MEA B4 ) and NDTI(MEA B5 , MEA B2 ) were the optimal texture indices across the whole growing season in exponential regression for LNA and PNA estimation, respectively (Table 5). Because allometric variation of LNA or PNA in the season was largely attributed to biomass (Table 3), texture indices could characterize the biomass related variation with an exponential regression [7,8].
Furthermore, the variables in the multivariate model for LNA included one N-sensitive variable (NDTI(DIS B4_D3 , COR B5_D3 )) and two biomass-sensitive variables (NDRE and NGRDI). Therefore, LNA could reflect part of the variation in LNC [35,36], but the exact proportion needs to be investigated in future work. However, the optimal multivariate PNA estimation model consisted of entirely AGB-sensitive indices such as NDTI(MEA B5_D3 , MEA B2_D4 ), MTVI2 [29], and NGRDI [20]. This also suggests that the variation in PNA might be dominated by biomass and contain weak information on N concentration which corresponds well with the findings in Reference [5].

Directional Effect of Texture Analysis on Row-Planted Crops
The directional effect of texture analysis was rarely investigated in the existing literature, since the majority of previous studies executed texture analysis with the default direction (45 • ) [7,15]. Some studies calculated texture metrics with different directions but did not explicitly explain the reason [10,16]. That might be because most of them studied naturally grown forests, and the trees were distributed disorderly. In contrast, we found that texture metrics had a significant directional effect (Figure 2). That is because rice plants are grown in rows, and the local window sliding along the row orientation contains more homogeneous vegetation than from other directions. Furthermore, the rice plants in the same row grow more homogenously than those in different rows, resulting in lower contrast and higher correlation in the along-row direction than in other directions ( Figure 5). Texture indices calculated with texture metrics from different directions had different performances on the N nutrition parameter estimation (Table 6). However, the texture feature MEA represents the average values within the moving window and the textural information extracted from the UAV images is the average value from an ROI in the non-sampling area. As a result, MEA was not affected by directions, and the same R 2 value was obtained with the optimal NDTIs composed of MEA features calculated in different directions for LNA and PNA estimation ( Table 6).
sensitive indices such as NDTI(MEAB5_D3, MEAB2_D4), MTVI2 [29], and NGRDI [20]. This also suggests that the variation in PNA might be dominated by biomass and contain weak information on N concentration which corresponds well with the findings in Reference [5].

Directional Effect of Texture Analysis on Row-Planted Crops
The directional effect of texture analysis was rarely investigated in the existing literature, since the majority of previous studies executed texture analysis with the default direction (45°) [7,15]. Some studies calculated texture metrics with different directions but did not explicitly explain the reason [10,16]. That might be because most of them studied naturally grown forests, and the trees were distributed disorderly. In contrast, we found that texture metrics had a significant directional effect (Figure 3). That is because rice plants are grown in rows, and the local window sliding along the row orientation contains more homogeneous vegetation than from other directions. Furthermore, the rice plants in the same row grow more homogenously than those in different rows, resulting in lower contrast and higher correlation in the along-row direction than in other directions ( Figure 5). Texture indices calculated with texture metrics from different directions had different performances on the N nutrition parameter estimation (Table 5). However, the texture feature MEA represents the average values within the moving window and the textural information extracted from the UAV images is the average value from an ROI in the non-sampling area. As a result, MEA was not affected by directions, and the same R 2 value was obtained with the optimal NDTIs composed of MEA features calculated in different directions for LNA and PNA estimation ( Table 5). Texture metrics calculated in directions D3 and D4 constructed the optimal NDTIs for LNC and PNC estimation which might be explained by the strong capability of across-row texture metrics in differentiating the tonal variations caused by N status. Moreover, the estimation accuracy (RMSE = 0.24, RRMSE = 10.52%) of LNC decreased marginally if the texture indices in the multivariate model of LNC were derived from the direction D3 alone. Additionally, the same estimation accuracies were obtained for other N parameters (LNA, PNC, and PNA) through the multivariate models with all the texture indices calculated in the direction D3 when compared to the original multivariate models Texture metrics calculated in directions D3 and D4 constructed the optimal NDTIs for LNC and PNC estimation which might be explained by the strong capability of across-row texture metrics in differentiating the tonal variations caused by N status. Moreover, the estimation accuracy (RMSE = 0.24, RRMSE = 10.52%) of LNC decreased marginally if the texture indices in the multivariate model of LNC were derived from the direction D3 alone. Additionally, the same estimation accuracies were obtained for other N parameters (LNA, PNC, and PNA) through the multivariate models with all the texture indices calculated in the direction D3 when compared to the original multivariate models (data not shown). Therefore, all N nutrition parameters could be estimated at nearly the highest accuracies with the combination of VIs and texture indices calculated with texture metrics in the direction D3 which could simplify the use of texture analysis significantly.
Texture analysis also involved the optimal selection of texture calculation algorithm, spectral band, and window size [14]. Although window size has a considerable effect on the estimation accuracy of forest biomass due to the mismatch of spatial scale between remotely sensed pixel size and tree canopy [14,16], it showed no significant influence in the present study. This was because most window sizes contained a large proportion of crops due to the large canopy coverage since the jointing stage. Rice crops are often planted in the row distance of 24-30 cm, and multispectral images collected at 100 m usually possess a spatial resolution of 5-7 cm. Hence, a larger window size might be a good choice for texture analysis at early growth stages. Yue et al. [7] also found the optimal image resolution for using image textures to estimate AGB in winter wheat depended on the crop canopy size and row spacing. Therefore, the optimal window size has to be taken into consideration in terms of row space and image spatial resolution for other row-planted crops (e.g., soybean, corn).

The Benefits of Fused Information for Enhancing N Signals
Numerous studies were dedicated to improving N concentration estimation with different approaches due to the importance of N nutrition status in crop management. The most commonly used way was to propose new VIs for N concentration estimation [37,38]. Although those VIs yielded high estimation accuracies as shown in the literature, the unstable performance was reported in other studies [39]. In this study, VIs had poor performance in N concentration estimation, which is in line with References [3,5]. Texture indices showed comparable performance when compared to VIs which is in contrast to the performance of texture ratios in forest AGB estimation [10,16]. However, a significant improvement in LNC estimation was obtained when using the combination of texture indices and VIs as compared to using VIs or texture indices alone, with an increase of more than 35% in RMSE ( Table 7). This finding agrees well with the results of References [7] and [8] which reported that the combination of texture indices and VIs improved the estimation of crop AGB significantly. In contrast to relevant studies [6,40], a universal model could be used to estimate N concentration across the entire season, which was a substantial improvement in crop N status monitoring with remotely sensed data. Compared with the new spectral index proposed by Stroppiana et al. [38] with hyperspectral data (R 2 = 0.65), our multivariate model derived from multispectral data could even yield a similar R 2 value (R 2 = 0.68) for PNC. A universal model suitable for the entire season could not only be used to guide N fertilization applications at the early growth stages [1,41] but also to predict crop yield and quality before harvest [42,43].
The improvement in LNC and PNC estimation induced by the addition of texture indices stemmed from the enhancement of N signals for the early stages. VIs could not be used to estimate N concentration at the early growth stages (from tillering to booting) with a universal model, because N signals were hampered by the rapid increase of biomass and the large proportion of background materials [5,6]. However, leaf color changed remarkably due to the N fertilizer applications conducted at jointing stage [11,12], and texture metrics could characterize the spatial distribution of tonal variations caused by N status. The values of CON feature had a greater variation between jointing and booting stages than reflectance ( Figure 6). Furthermore, VIs had weak capability in N status monitoring at reproductive stages due to the high canopy coverage [2]. However, leaf tone had a visual color change from the vegetative period to the reproductive period, and the values of CON and DIS were significantly different at reproductive stages ( Figure 6). NDTI has a wider variation than NDRE at high N levels at late growth stages ( Figure 7). Therefore, the N signal could also be detected by those texture metrics at the late growth stages.
Moreover, texture metrics have the capability of smoothing the spatial heterogeneity between vegetation and background materials with a sliding window [9,16]. Therefore, the complementary information between texture and VIs might be useful for solving the saturation problem and reducing background interference. Although a similar conclusion has been drawn by Yue et al. [7] and Zheng et al. [8] for AGB estimation, this study represents the first solid evidence on improving the estimation of N nutrition parameters, especially leaf and plant N concentrations.
Remote Sens. 2020, 12, x FOR PEER REVIEW 13 of 17 et al. [8] for AGB estimation, this study represents the first solid evidence on improving the estimation of N nutrition parameters, especially leaf and plant N concentrations.  Multi-source data fusion was commonly used in previous studies to improve the estimation of crop biomass [44][45][46]. For example, Yue et al. [45] improved wheat biomass estimation by combining VIs and plant height derived from crop surface models (CSMs), but the construction of CSMs was tedious and time consuming. Although the fusion of data from different platforms could also acquire Remote Sens. 2020, 12, x FOR PEER REVIEW 13 of 17 et al. [8] for AGB estimation, this study represents the first solid evidence on improving the estimation of N nutrition parameters, especially leaf and plant N concentrations.  Multi-source data fusion was commonly used in previous studies to improve the estimation of crop biomass [44][45][46]. For example, Yue et al. [45] improved wheat biomass estimation by combining VIs and plant height derived from crop surface models (CSMs), but the construction of CSMs was tedious and time consuming. Although the fusion of data from different platforms could also acquire Multi-source data fusion was commonly used in previous studies to improve the estimation of crop biomass [44][45][46]. For example, Yue et al. [45] improved wheat biomass estimation by combining VIs and plant height derived from crop surface models (CSMs), but the construction of CSMs was tedious and time consuming. Although the fusion of data from different platforms could also acquire higher estimation accuracy [44], the data (e.g., LiDAR data) processing needs more professional skill.
However, UAV-based image texture and spectral data are more convenient to acquire over small fields. Furthermore, the fusion of texture and spectral data could improve the accuracy of N status estimation of multiple growth stages, which solves the problem that N status was often estimated at only one growth stage or a short time window [6,40]. Therefore, monitoring N status across the whole season could be realized through the integration of texture and spectral information into the crop growth monitoring systems.
There might be some concerns that it is challenging to build a universal model for predicting the N status of rice plants across the entire growing season [3]. With spectral information alone, VIs often exhibit weak capabilities of detecting N status at late growth stages due to the fact of crop senescence. Textural information could capture the variation of leaf color within and between plots, especially at the late stages when leaf color changes dramatically. However, our results demonstrated that the fusion of spectral and textural information could improve N status monitoring efficiently. Furthermore, the multivariate models for LNC and PNC were tested on different dataset groups with satisfactory validation performance ( Table 8). The applicability of those models still needs to be improved through further testing with more datasets from different geographic sites.

Potentials for Other Platforms
In this study, we found the combination of spectral and texture information in NIR and RE bands was superior to that in other bands in rice N status monitoring. As for monitoring crop N status in large areas, satellite imagery has been widely used but with low estimation accuracy. That is because satellite imagery has low spatial resolution and low spectral resolution normally with blue, green, red and NIR wavebands. RapidEye is the first launched satellite with a RE band and successfully applied to many aspects of precision agriculture [47,48]. Other satellites with RE bands (e.g., Sentinel-2, WorldView-2, and Gaofen-6) have been launched and have yielded significant improvement in the estimation accuracy on agronomic variables [49,50]. These findings from satellite observations were consistent with our results that the RE-based VIs were superior to other indices.
To date, texture analysis on satellite imagery was mostly used for forest AGB [14][15][16]. Our findings provide a strong reference on crop N status with texture analysis of satellite images equipped with the RE band. However, in terms of spatial resolution and crop row spacing, the spatial resolution of satellite imagery (highest spatial resolution is 0.3 m by WorldView 3) is insufficient for distinguishing the narrow row spacing in cereal crops (e.g., rice, wheat). Furthermore, crops could be planted in multiple directions over large areas and the direction effect of texture metrics may not be significant in satellite images with lower spatial resolutions. Given the suitability of UAV platforms for small-scale applications for the moment, the direction of texture analysis could still exist in UAV imagery. When using texture analysis on row-planted vegetation (e.g., staple crops, vegetable crops, and fruit trees) at the farm level with high resolution images, it is still beneficial to take consideration of the direction effect for improved nitrogen nutrition monitoring.

Conclusions
We investigated the potential of texture analysis from UAV-based multispectral imagery for N status monitoring in terms of texture metrics, texture directions, and texture indices. For N concentration estimation, low estimation accuracies were obtained using VIs, individual texture metrics or texture indices alone. However, considerable improvements were achieved with the combination of VIs and texture indices. The multivariate models with fused information yielded the highest accuracy for LNC (RMSE = 0.22 and RRMSE = 9.70%) and PNC (RMSE = 0.19 and RRMSE = 15.34%). Moreover, the multivariate models exhibited high estimation accuracy when tested on dataset in different years and growth stages. For N accumulation estimation, NDRE yielded high estimation accuracies on LNA (RMSE = 1.35 g m −2 and RRMSE = 24.56%) and PNA (RMSE = 3.16 g m −2 and RRMSE = 32.20%). Compared with the models based on VIs or texture indices alone, the multivariate model also yielded higher accuracies for LNA (RMSE = 1.16 g m −2 and RRMSE = 21.04%) and PNA (RMSE = 2.38 g m −2 and RRMSE = 24.17%) estimation. When calculating the texture metrics, the across-row direction (D3) should be used for the best performance over row-planted crops instead of the default diagonal direction (D2). These findings could serve as useful references for deriving appropriate texture indices from multispectral UAV imagery for N nutrition monitoring. This study provides new insights on enhancing the N signals from rice canopies by adding the textural information to conventional VIs. It has great potential in conducting rapid, accurate monitoring of N status monitoring over other crops and developing effective practices of improved crop management.