Next Article in Journal
Comparison of Phenology Estimated from Reflectance-Based Indices and Solar-Induced Chlorophyll Fluorescence (SIF) Observations in a Temperate Forest Using GPP-Based Phenology as the Standard
Previous Article in Journal
Intercomparison and Validation of SAR-Based Ice Velocity Measurement Techniques within the Greenland Ice Sheet CCI Project
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Multi-Temporal and Spectral Analysis of High-Resolution Hyperspectral Airborne Imagery for Precision Agriculture: Assessment of Wheat Grain Yield and Grain Protein Content

1
International Maize and Wheat Improvement Center—CIMMYT, Texcoco 56237, Mexico
2
Food and Rural Development, School of Agriculture, Newcastle University, Newcastle NE1 7RU, UK
3
Earth and Life Institute, Université Catholique de Louvain, Croix du Sud L5.07.16, B-1348 Louvain-la-Neuve, Belgium
4
International Maize and Wheat Improvement Center—CIMMYT, Henan Agricultural University, 63 Nongye Road, Zhengzhou 450002, Henan, China
5
Instituto de Agricultura Sostenible (IAS), Consejo Superior de Investigaciones Científicas (CSIC), 14004 Cordoba, Spain
*
Author to whom correspondence should be addressed.
Remote Sens. 2018, 10(6), 930; https://doi.org/10.3390/rs10060930
Submission received: 14 May 2018 / Revised: 5 June 2018 / Accepted: 8 June 2018 / Published: 12 June 2018
(This article belongs to the Section Remote Sensing in Agriculture and Vegetation)

Abstract

:
This study evaluates the potential of high resolution hyperspectral airborne imagery to capture within-field variability of durum wheat grain yield (GY) and grain protein content (GPC) in two commercial fields in the Yaqui Valley (northwestern Mexico). Through a weekly/biweekly airborne flight campaign, we acquired 10 mosaics with a micro-hyperspectral Vis-NIR imaging sensor ranging from 400–850 nanometres (nm). Just before harvest, 114 georeferenced grain samples were obtained manually. Using spectral exploratory analysis, we calculated narrow-band physiological spectral indices—normalized difference spectral index (NDSI) and ratio spectral index (RSI)—from every single hyperspectral mosaic using complete two by two combinations of wavelengths. We applied two methods for the multi-temporal hyperspectral exploratory analysis: (a) Temporal Principal Component Analysis (tPCA) on wavelengths across all images and (b) the integration of vegetation indices over time based on area under the curve (AUC) calculations. For GY, the best R2 (0.32) were found using both the spectral (NDSI—Ri, 750 to 840 nm and Rj, ±720–736 nm) and the multi-temporal AUC exploratory analysis (EVI and OSAVI through AUC) methods. For GPC, all exploratory analysis methods tested revealed (a) a low to very low coefficient of determination (R2 ≤ 0.21), (b) a relatively low overall prediction error (RMSE: 0.45–0.49%), compared to results from other literature studies, and (c) that the spectral exploratory analysis approach is slightly better than the multi-temporal approaches, with early season NDSI of 700 with 574 nm and late season NDSI of 707 with 523 nm as the best indicators. Using residual maps from the regression analyses of NDSIs and GPC, we visualized GPC within-field variability and showed that up to 75% of the field area could be mapped with relatively good predictability (residual class: −0.25 to 0.25%), therefore showing the potential of remote sensing imagery to capture the within-field variation of GPC under conventional agricultural practices.

1. Introduction

Wheat (Triticum sp.) is one of the three most important cereals produced worldwide, along with maize (Zea mays) and rice (Oryza sp.). It is also one of the most important crops in Mexico, grown on more than 723,559 ha in 2016, with average yields of 5.3 Mg ha−1 [1]. Crop management is key to profitable and sustainable wheat production and must be able to respond to spatial and temporal variability in soil and climate with the aim of improving grain yield (GY) and quality.
Accurate and timely assessment of within-field GY and grain protein content (GPC) variations over the crop cycle generates different degrees of potential for N and water management/availability as indicated by Reference [2] and Reference [3], and also as a strategy for selective harvesting [4,5,6,7]. GY and GPC are two major factors for wheat production where the latter is an important determinant of the end-use value [8]. GPC is a function of the conversion of grain nitrogen (N) content into protein (further reading in References [9,10]), which is dependent on genotype and strongly influenced by environmental variables, such as timing and amount of nitrogen application, water access and temperature, especially during the grain filling period [11,12,13,14,15]. These factors influence the rate and duration of wheat grain development, protein accumulation and starch deposition [16,17,18]. The most influential environmental factor on wheat quality is the availability of soil N, which in turn is influenced by N fertilization. Therefore, proper management of N fertilizer is essential to ensure high quality wheat production.
The design of fertilizer application regimes should combine different factors (e.g., rate, placement, timing and splitting) with a view to optimize both wheat GY and GPC [19]. For example, the application of N fertilizer around heading increases GPC without reducing GY [20,21]. Hence, it is important to link N fertilization with N uptake by the plant.
The highest N uptake begins around jointing and reaches a plateau at heading. N uptake, N accumulation and the further partitioning of N within the plant are all important processes that determine GY and GPC [22,23]. N uptake, which has major influence on a crop’s green canopy and on the way N controls canopy growth and senescence, depends upon the stage of crop development—before, during and after stem extension. Canopy size determines the proportion of sunlight intercepted, and is directly related to dry matter and green biomass, with the latter shown to be strongly related to GY [24]. Differences in N uptake between the pre- and post-anthesis periods may affect N partitioning at the plant level [21] and N content in the grain as this originates from two different sources: partitioning from vegetative organs of N assimilated at pre-anthesis stages, and N uptake from the soil post-anthesis. It is understood that N accumulated before anthesis is the major source of grain N. In wheat, between 50–95% of the N content in the grain comes from the partitioning of N stored in shoots and roots before anthesis, with the leaves and stems being the most important sources of N for grain [25].
Crop nitrogen status has been widely reported to be linked to canopy reflectance because of its relation to chlorophyll content (Chlorophyll a and b; Ca+b) [3,26,27]. Remote sensing technologies have become less expensive in recent years, and this has improved their potential for operational purposes. Remote sensing-based approaches to vegetation monitoring use two main methods to estimate canopy biophysical and biochemical traits. The first is based on parameter retrievals through vegetation canopy reflectance modeling, while the second approach focuses on developing empirical relationships between remotely sensed canopy reflectance and biophysical variables.
The second approach, which is our focus, seeks to establish empirical relationships between plant biophysical variables and vegetation indices (VIs). VIs are essentially a way of combining multispectral observations in a single metric, aiming to minimize the effect of external factors on spectral data and to derive canopy characteristics [28], whereas these external factors influencing the reflectance values may come from soil brightness, soil color, atmospheric effects, sensor calibration and differences in the spectral responses of sensors to bidirectional effects [29].
A large number of VIs have been developed to minimize those external effects to better convey information about vegetation canopy [29,30]. Indices based on both ratios and normalization formulas are referred to as Ratio Spectral Indices (RSIs) and Normalized Difference Spectral Indices (NDSIs), of which many variants have been developed that combine wavelengths of different parts of the electromagnetic spectrum (EMS).
Chlorophyll related VIs have also been developed due its relation to crop N, as previously mentioned. Combinations of VIs, such as Modified Chlorophyll Absorption in Reflectance Index (MCARI) divided by Optimized Soil-Adjusted Vegetation Index (OSAVI), have shown good results for describing leaf chlorophyll concentration [31]. Also, more pigment specific VIs such as Pigment Specific Normalized Difference c (PSNDc), Pigment Specific Simple Ratio for chlorophyll a and b (PSSRa; PSSRb), and carotenoids (PSSRc) [32] have shown to be indicators of crop N status. Based on this, canopy reflectance at specific crop stages may be a proxy for GY and GPC and may be indirectly measured using remote sensing technologies.
There is spatial variability of GY and GPC within a crop field due to soil spatial variability, micro-climate conditions and crop management [4,33,34,35]. Crop yield assessment using remote sensing has been considered useful [36,37,38,39,40]; however, it is well documented that some normalized indices, such as the Normalized Difference Vegetation Index (NDVI) [41], saturate at high leaf area index (LAI) values and are also affected by other factors, such as soil background, canopy shadows, illumination, atmospheric conditions and variation in leaf chlorophyll concentration [28,42,43]. For GPC, there are only a few studies relating it with satellite imagery [2,3,44,45,46] or low-altitude aerial and terrestrial proximal sensing [2,27,47,48,49].
The majority of those studies are based on multispectral broad-band signals on experimental designs or at the regional scale. There is usually wider range variation of the response variable at the regional scale than at the field scale [50]. In several studies, this wide range was induced, for example, in the case of GPC with different N application levels [2,3,27,46,47,49]. To the best of our knowledge, in the literature there is no detailed description of in-situ GPC data estimated through remote sensing imagery. There are also only a few studies that have explored how multi-temporal information can enhance the assessment of GY and/or GPC at a regional scale using broad-band multispectral imagery [46,49,51,52].
New methods that use hyperspectral remote sensing promote further spectral exploration of the signal, allowing the calculation of several other narrow-band VIs, suggested as potentially useful for precision agriculture [42,53]. Besides the VIs currently presented in the literature, hyperspectral signals can be used to calculate complete combinations of all available wavelengths using generic formulas, e.g., NDSI and RSI [54,55,56,57,58,59]. Exploiting the hyperspectral signal, multi-temporal spectral exploratory analysis offers new potential for assessing within-field information, and identifying the most suitable combinations of spectral wavelengths and the optimal time of image acquisition for specific precision agriculture applications.
The present study analyses the respective contribution of the spectral, narrow-band physiological vegetative indices and temporal dimensions of the hyperspectral signal, to the assessment of the within-field spatial variability of GY and GPC in two commercial wheat fields. The feasibility of this study relies on a multi-temporal and spectral exploratory analysis.

2. Materials and Methods

2.1. Field Site and Data Collection

This study was carried out in a wheat field in the Yaqui Valley near Ciudad Obregón (Sonora), in northwestern Mexico (27°23′43.83″N and 109°55′0.90″W), during the 2014 wheat crop cycle. It consisted of two furrow irrigated 43-ha blocks (called ‘A’ and ‘B’) (Figure 1). The blocks were sown with the variety Cirno-C2008 in January 2014 and harvested at the end of May 2014. The climate in the Yaqui Valley is semi-arid (Köppen climate classification subtype “Bsh”) with variable precipitation rates averaging 280 mm per year and an average daily temperature of 24 °C. The soils in this region are coarse sandy-clay, mixed with montmorillonitic clay. The two major soil types in the Yaqui Valley are clayey and alluvial soils at a 3:2 ratio [60].
A weekly/biweekly flight campaign took place starting from stem elongation (GS31; February 14; [61]) until just prior to harvest (GS92; grain hard, 7 May), resulting in 10 airborne mosaics (Table 1). They were acquired with a push-broom micro-hyperspectral imaging sensor (Micro-Hyperspec VNIR model, Headwall Photonics, Bolton, MA, USA) (spectral region: 400–850 nm; 250 channels) flying at 1200 m above ground in a manned airplane, yielding a ground sampling distance of circa 1 m.
Ground control points (GCPs) were made using 3 × 3 m long × 1 m thick crosses of polyethylene sheet coated with aluminum foil fixed with duct tape and were placed around the study sites (Figure 1) and georeferenced with a global navigation satellite system (GNSS) receiver using the real-time kinematic (RTK) technique (Trimble R4 GNSS system, Trimble, Sunnyvale, CA, USA) to ensure accuracy. Georeference processing of each image was done ensuring root mean square errors (RMSE) lower than the image resolution.
The micro-hyperspectral instrument was radiometrically calibrated in the laboratory using derived coefficients with a calibrated uniform light source (integrating sphere, CSTM-USS-2000C Uniform Source System, LabSphere, North Sutton, NH, USA) at four levels of illumination and six integration times. Hyperspectral imagery was atmospherically corrected using the total incoming irradiance at 1 nm intervals simulated with the SMARTS model developed by the National Renewable Energy Laboratory, US Department of Energy [62,63]. Therefore, the aerosol optical depth was measured at 550 nm with a Micro-Tops II sun photometer (Solar LIGHT Co., Philadelphia, PA, USA) in the study area at the time of the flights. SMARTS computes clear sky spectral irradiance, including direct beam, circumsolar, hemispherical diffuse, and total irradiance on a tilted or horizontal plane for specified atmospheric conditions. The algorithms were developed to match the output from the MODTRAN complex band models to within 2%, using aerosol optical depth as an input. The spectral resolution was 0.5 nm for the 280–400 nm and 1 nm for the 400–1750 nm ranges of the electromagnetic spectrum. This radiative transfer model has been previously used in other studies [64,65,66,67,68] for the atmospheric correction of narrow-band multispectral imagery.
Manual grain sampling [69] took place just before harvest using a half regular and a half stratified grid of 50 sampling points in Block B (100 points in total). Reallocation of 10% of the regular grid points to shorter distances than the original grid was done to minimize the ratio of the smallest to largest sample distance [70,71]. This aimed to define the range of spatial dependence [72] of the response variables. In Block A, 14 sampling points were selected based on a visual inspection of the soil apparent electrical conductivity (ECa) map [73], to cover the full ECa range variation. The decision for different number of points in each block was due to budget constraints and to have at least a reasonable number of points (100) in one of the fields to capture its spatial variability and enable mapping using variography [74]. Each sampling point was based on a 2 m2 frame centered on the point geocoordinates, where all wheat plants were harvested, threshed and GY measured. A grain sub-sample was taken for laboratory quality analysis, where GPC (%) and moisture content (%) were determined by NIR spectroscopy (NIR Systems 6500, Foss, Hilleroed, Denmark) calibrated according to official AACC standard methods 39–10 and 46–11A [75]. The GY and GPC values were reported at 12.5% moisture basis. Descriptive statistics of the sample were computed using the corresponding variables.
Spectral binning was performed on each mosaic into 7.5 nm FWHM (Full Width at Half Maximum) to decrease noise effects, resulting in 62 wavelengths. From those, the 751, 759, 766, and 773 nm wavelengths were removed due to oxygen absorption by the sensor. Finally, 58 wavelengths were used for subsequent analyses. Spectral information was extracted from each image using the point sampling location from both blocks (n = 114). The extraction was done by taking the average of a 3 × 3-pixel window (9 m2) around each sampling point to account for minor shift of the images during orthomosaic processing. The resulting averaged spectral data set was used for the following statistical and remote sensing analyses regarding GY and GPC.

2.2. In Situ Data Description—GY and GPC Descriptive Statistics, Correlation Analysis and Hyperspectral Profiles

The GY and GPC data from Block B were interpolated onto a 3-m grid (pixels of 9 m2) using global block kriging, fitting the best global variogram according to the spatial variability of the data in VESPER [76]. The purpose of this mapping was to visualize the within-field spatial variability of the response variables. Map displays were done using the ArcGIS software suite (v10.1; ESRI, Redlands, CA, USA).
The relationship between GY and GPC was analyzed using Pearson’s correlation coefficient. As stated by many authors [4,21,77,78,79,80], GY is usually negatively correlated to GPC. However, [80] demonstrated that this correlation may vary from negative to positive within a farmer’s field due to soil spatial variability and crop management interactions. Based on this context, we decided to further explore the within-field relationship of GY and GPC by applying a moving window approach to check the spatial local correlation between both variables. The moving window approach was applied only to block B due to its denser data sampling (100 sampled points). It consisted of a 150 × 150 m window which moved across the block collecting the neighbors within this window, using as basis the same 100 sampled point locations. The window dimension was chosen to respect an average of 10 degrees of freedom for the correlation analysis. Correlation results were considered significant if p < 0.05.
The hyperspectral profiles of the highest and lowest GY and GPC sampling points only were extracted for each image and plotted on a single graph for each response variable. The choice to limit the graph to these two responses, and the decision not to plot the mean or median response, was made to visualize the possible reflectance differences for both levels (highest and lowest) of the response variables without generating an overcrowded and unintelligible graph.
Two exploratory analysis approaches were selected for data analysis of the acquired hyperspectral imagery time series: (a) a spectral approach; and (b) a multi-temporal approach. The spectral exploratory analysis approach is based on a complete two by two combination of all available wavelengths into NDSI and RSI spectral indices (SI). For the multi-temporal exploratory analysis approach, two methods were tested: temporal principal component analysis (tPCA) and the integration of VIs over time using an area under the curve (AUC) approach.

2.3. Spectral Exploratory Analysis Using Narrow-Band Physiological Spectral Indices

For the spectral exploratory analysis, we applied two types of generic formulas to generate SIs. Using complete two by two combinations of spectral wavelengths, the ratio spectral index (RSI) and the normalized difference spectral index (NDSI) were calculated. The RSI (Equation (1)) and NDSI (Equation (2)) are defined as:
RSI ( i , j ) = R i R j
NDSI ( i , j ) = R i R j R i + R j
where Ri and Rj are the reflectance for wavelengths i and j, respectively. Representations of RSI and NDSI are the PSSRa and NDVI, whereby both indices use NIR and R wavelengths. Spectral indices based on the complete two by two combinations of the hyperspectral signal were generated by applying the RSI and NDSI formulas, similar to studies by [54,55,56,57,58,59]. Subsequently, regression analyses were performed using RSI and NDSI indices as predictors for all GY and GPC data points. Contour maps of the coefficient of determination (R2) describing the relationships of these SIs with GY and GPC were generated. These maps provide inclusive information on the optimum pair of wavelengths to assess both response variables all along the crop cycle [54,55].

2.4. Multi-Temporal Spectral Exploratory Analysis

2.4.1. Temporal Principal Component Analysis (tPCA)

The relatively high temporal resolution flight campaign (here: 10 hyperspectral mosaics within one crop cycle—February to May) allowed for a multi-temporal approach to analyze the importance of the mosaics and wavelengths all along the crop cycle and at specific phenological stages (crop growth stages). For this purpose, a standardized principal component analysis (PCA) was applied to each individual wavelength across the 10 images, with 10 principal components (PCs) recorded for each wavelength. The first two PCs were considered in a Pearson’s correlation analysis with GY and GPC. The eigenvector and eigenvalue matrices were assessed to determine which image and wavelength across the crop cycle carried the highest proportion of variance within the dataset.

2.4.2. Integration of VIs Over Time

For further data exploration, four specific VIs for each response variable were chosen according to their potential to deliver structural and chlorophyll information (Table 2). Also, the relationship of VIs with GY and GPC was a selection criterion, as presented in preliminary results by [73,81]. Both canopy structure and its greenness are indicators of the N-status of the plants, where the leaves and stems are stated to be the most important sources of N for the grain [25].
The VIs’ temporal profiles based on the highest, median and lowest GY and GPC points within the field were plotted. Based on this information, the area under the curve (AUC) using the trapezoidal method (Riemann’s integrals) was calculated for each selected VI on each sampling point observation, aiming to integrate its temporal information across the crop cycle for different crop growth stages (Table 3).
The AUC for each VI and crop growth stage were calculated and the area regressed against GY and GPC using linear regression (software R: package “plyr” for AUC and package “stats” for regression analysis; [86,87]). Regression results were considered significant if p < 0.05.

3. Results and Discussion

3.1. In Situ Data Description—GY and GPC Descriptive Statistics, Correlation Analysis and Hyperspectral Profiles

Both GY and GPC showed a distribution with mean and median values close to each other. This reflects a broadly symmetrical (normal) distribution although with a high value of kurtosis for GPC (Table 4). The GPC data are concentrated between 12.03% (1st Quartile) and 12.56% (3rd Quartile) with the mean of 12.32%. Relatively low GPC ranges might occur for irrigated fields under natural conditions and conventional agricultural practices, as N and water are applied equally and at high quantities throughout the field and soil spatial variability may have a limited effect on canopy variability. However, higher GPC ranges have been reported in rain-fed crop management fields, where soil spatial variability strongly influences water and N availability [80].
Through the mapping exercise (Figure 2), it is possible to visualize the within-field spatial variability for GY and GPC, which ranged from around 4.6 to 8 Mg ha−1 and 11 to almost 15%, respectively. The threshold of 12.5% for GPC was chosen based on the farmer’s recommendation to differentiate between low and high wheat quality based on premium quality and better market prices. Although the maximum value of 14.96% for GPC may be considered an outlier, after a few analyses (data not shown), we decided to keep the data point. There was no further evidence this value was an error of measurement. GPCs of around 14–15% within farmers’ fields and on a regional scale are reported throughout the literature [4,45,80]. The highest and lowest GY and GPC points were identified among the sampled data points and maps, and afterwards used as references for extracting the reflectances from the images during the crop cycle.
GY and GPC did not show significant correlation (r = −0.09; p-value = 0.49) when we used all data points from the entire field. However, local spatial variation in the correlation between GY and GPC obtained through the moving window approach showed correlation coefficients varying from −0.76 to +0.61. Just two negative coefficients were significant at 5% probability (−0.76 and −0.69) and one positive (0.61) at 10% probability among the other 97 coefficients (Table 5).
These results agree with results obtained by Reference [80], where GY and GPC data (acquired from on-the-go sensors mounted on harvesters in 27 fields over 3 seasons) were evaluated. The majority of these fields showed bigger areas of negative local correlations for the first two seasons and positive for the third season; however, a considerable percentage of the area had non-significant coefficients for all seasons/fields.
Areas within the field where GY and GPC were negatively correlated could represent conditions where effective access to N has been relatively uniform, but with limited access to water due to soil spatial variability structure, which limited grain growth. Areas with positive coefficients could represent more access to water but limited N availability [80]. Soil apparent electrical conductivity (ECa), which may be a proxy for plant available soil water storage capacity (PAWC, [88]), was quite variable in this field [73]. This supports the within-field spatial variation of the correlation between GY and GPC and reinforces the need for exploring both variables separately.
The previously identified highest and lowest GY and GPC locations (Figure 2a,b) were used to extract the hyperspectral signal from the mosaics. The reflectance profile of the highest (dashed lines) and lowest (full lines) GY region is represented in Figure 3a. The dark red dashed and full lines (14 February) represent the GS31 stage, which is critical, as it is a stage where N is often applied and an early stage where nutrient stress diagnosis in wheat may be made. The dark green (dashed and full) lines represent anthesis (28 March), supposedly the peak of vegetative development. At GS31, a small reflectance difference in the NIR region is detectable. This difference increases up to 15 April, and decreases at physiological maturity (GS87, 25 April) and grain hard (GS92, 7 May).
The reflectance profile (Figure 3b) from the highest (dashed lines) and lowest (full lines) GPC region showed similar spectral behavior from 400 to 770 nm, except for the mosaic from 25 April (end of grain filling stage), which showed some difference in reflectance in the R region (around 680 nm). At 840 nm (NIR), it showed some differences in the reflectance for the highest and lowest GPC region.
Zarco-Tejada et al. [42] found similar reflectance behavior for low and high growth areas of cotton. In the case of wheat, the best time to diagnose for GY is at the beginning of stem elongation (GS31) and for GPC around heading/anthesis stage, as that is when GPC can be increased through nitrogen management [20]. Therefore, identifying differences in reflectance at an early stage (GS31) and at heading/anthesis makes it possible to diagnose nitrogen stress while there is still time to raise GY and GPC through crop management.

3.2. Spectral Exploratory Analysis Using Narrow-Band Physiological Spectral Indices

Since the results of the regression analysis using RSI and NDSI were very similar, only NDSI results are shown here. The contour maps of the coefficients of determination (R2) are shown in Figure 4 and Figure 5, which are the results of the regression analyses of the NDSI indices versus GY and GPC, respectively.
Using the R2 contour maps, it is possible to infer which wavelength combination performed better when assessing GY and GPC (Figure 4 and Figure 5). The relationship between crop yield and VIs has been widely studied in the literature. Many authors have reported good correlations between the yield of major commodity crops such as wheat, maize, and cotton based on multispectral broadband data imagery and mostly using combinations of NIR, R and green (G) spectral regions [36,37,38,39,40,89,90]. In the current study, mosaics around booting and anthesis led to the highest relationship with GY (Figure 4). The R2 varied from 0.22 to 0.32 (booting to anthesis; GS41–GS65; 11–28 March). The best results across these images were with combinations between broad parts of the NIR region (Ri, 750 to 840 nm—except 751, 759, 766, and 773 nm which were previously removed) and narrow bands of the RE (Rj, ±720–736 nm) region. In these regions, reflectance is linked with biomass, which has been shown to be related to GY [91,92,93]. These results are followed by combinations of broad parts of NIR and the visible region of the spectrum. Here, all wavelengths between 400 and 680 nm (Rj) combined with NIR obtained slightly lower adjustments than NIR and RE combinations with peaks at NIR (Ri) and blue (B)/G regions (Rj). Similar results were highlighted in the studies by [30,94].
As well as revealing combinations between wavelengths in the NIR, RE and visible spectral regions, relationships between combinations solely within the visible spectrum were also detected. Although with a lower R2, G/R regions (Ri) combined with B (Rj) showed consistent relationships throughout the three images acquired from booting to anthesis. The correlations between GY and the combinations within the visible spectrum are supported by the fact that chlorophyll (Ca+b) has two peaks of light absorption in the B and R regions [95]. Also, carotenoids display peak absorption in the B region [96].
GPC is a function of the conversion of plant nitrogen content into protein, so it is expected that plant nitrogen concentration estimated through remote sensing techniques may be able to partially explain GPC variation [45]. A few studies have reported the potential use of individual wavelengths and/or different VIs such as PPR ([R550-R450]/[R550 + R450]), NDVI ([R810-R680]/[R810 + R680]), RVI (R810/R680), GNDVI ([R810-R560]/[R810 + R560]) and GRVI (R810/R560) to describe nitrogen status and GPC [3,27,45,46,47,49]. All of the former studies have used VIs derived from G, R and NIR wavelengths through normalized and simple ratio formulas. In our study, NDVI, for instance, was not one of the top 10 coefficients with GPC across the images and growth stages. However, most previous studies were carried out under controlled conditions with low and high levels of nitrogen application (plot or trial design). In contrast, this study was carried out in an irrigated farmer’s field without nitrogen level treatments (respectively under natural conditions), resulting in a low range of GPC variability (12.03–12.56% from the 1st to 3rd Quartile; Table 4). This may explain the rather low coefficients of determination (R2 ≤ 0.21) obtained in the results.
Through the R2-contour maps (Figure 5), it is also possible to visualize the wavelength region where R2 was predominantly high, such as from 574 to 700 nm on Ri and 400 to 574 nm on Rj at pre-booting (27 February), and 707 to 840 nm on Ri and 486 to 530 nm on Rj after physiological maturity (7 May). All dates show slight R2 differences within those regions of the spectrum. The highest R2 within the images from 27 February (pre-booting; GS41) (R2: 0.20) and 7 May (grain hard; GS92) (R2: 0.21) were obtained with narrow band combinations of the RE and G regions (27 February: 700, 574 nm; 7 May: 707, 523 nm). At booting, plants are increasing N uptake, and physiological maturity is the end of grain filling period, where N is redistributed from photosynthetic tissues to form GPC [97]. Although the complete two by two combinations of wavelengths provided improved SIs for the assessment of GPC against more traditional VIs, the results indicate that the visible near-infrared (VNIR) portion of the electromagnetic spectrum may not have enough potential to assess GPC at the field level. The shortwave infrared (SWIR) portion of the electromagnetic spectrum, where water content is the main determinant of leaf spectral behavior [98], may be of use for this specific trait. For technology development, the obtained information should support the wavelength decision for sensor development aimed at measuring GPC; however, different environmental factors, such as drought and heat, may also need to be taken into consideration.
In summary, the complete two by two combinations of spectral wavelengths (presented as R2-contour maps) proved to be a simple and robust method for exploring hyperspectral data. Here, combining NIR and RE spectral wavelengths into a normalized formula was shown to be more suitable for assessing GY and combinations of RE and G for assessing GPC. For the assessment of GY and GPC, better results were achieved with images acquired during anthesis (GY) and around booting and just after physiological maturity (GPC).

3.3. Multi-Temporal Spectral Analysis

3.3.1. Temporal Principal Component Analysis

The contribution to the total variance of each obtained PC for each wavelength is shown in Figure 6. Across all wavelengths, all PCs followed a similar course of fluctuation with (a) no data noise, (b) lower percentages at wavelengths in the violet (V) (~400 nm), G (~530 nm), and red-edge (RE) (~720 nm) regions, and (c) higher percentages at wavelengths in the B (~485 nm), R (~685 nm), and NIR (~830 nm) regions of the electromagnetic spectrum. Nevertheless, as the first PCs contribute more than the last PCs, generally most data noise is stored in the last PCs (here: PC6 to PC10) due to the PCA transformation process. PC1 (40–50% of the variation in each wavelength, except in some wavelengths in the V, G and RE regions) and PC2 (>25%, except in some wavelengths in the V, G and RE regions as well) explained the majority of the variation in each wavelength. PC1 made its highest contribution in wavelengths in the R and NIR region, and PC2 just in the R region. Together, >75% was obtained in wavelengths between 610 nm and 700 nm (R region). According to these results and taking into consideration Kaiser’s criterion [99] for PC selection (eigenvalue > 1), PC1 and PC2 of each wavelength were selected for subsequent statistical analyses.
To determine which mosaic across the crop cycle carries the highest proportion of variance within the dataset, the eigenvector and eigenvalue matrices were computed. For PC1 and PC2, Figure 7a,b depicts the percentage contribution to the total variance per wavelength of each image.
For PC1 (Figure 7a), all images (except 25 April and 7 May) followed the previously described course of fluctuation across all wavelengths with peaks in the B, R and NIR regions and depressions in the V, G and RE regions. Their contributions range almost equally between ~10 and 18%. Continuously over almost all wavelengths, the highest percentages (~18%) were obtained from images from 11 March (booting; GS41) and 17 March (heading; GS55), followed by images from 14 February (stem elongation; GS31) to 27 February (pre-booting; GS41), 28 March (anthesis; GS65) and 7 April (grain filling; GS71) (~10 to 15%). Very low percentages (~5%) were obtained from images taken on 15 April (late milk; GS77) to 7 May (grain hard, GS92), with the exception of the RE and NIR regions from the image of April 15 (~10%). The B and R spectral region is known for having absorbance peaks for Ca+b. NIR has high reflectance values from healthy vegetative organs with strong internal cellular geometry. At the booting stage (GS41), wheat plants are at approximately 40% of maximum growth and N uptake increases considerably. Meanwhile, the heading stage (GS55) comes just before the plant reaches its maximum green area index, which influences canopy reflectance and N uptake and reveals the importance of both stages in crop development [97].
In contrast to PC1, Figure 7b reveals that—based on PC2—late-season images between 15 April (late milk; GS77) and 25 April (physiological maturity; GS87) have a very high proportion of variance (~20 to 30%) across almost all wavelengths (except the NIR region of 15 April). Also, the formerly (based on PC1) lower contributing images from 14 February (stem elongation; GS31) to 27 February (pre-booting; GS41) show high percentages with peaks in the V, G and RE regions (~10 to 20%). Lower percentages (<10%) were found from images acquired between 11 March (booting; GS41) and 7 April (grain filling; GS71) across almost all wavelengths, except the RE (~720 nm) and NIR (~830 nm) regions. At the late milk stage (GS77), senescence and the rapid redistribution of soluble reserves have begun, which defines grain size and weight, while after physiological maturity (GS87), the grain will continue to lose moisture until it is ready for harvest.
Pearson’s correlation coefficient was calculated to understand the relationships of each wavelength from the selected PCs with GY and GPC, respectively. With GY, PC1 had weakly significant (−0.4 ≤ r ≤ −0.3) correlations across almost all wavelengths from the B to the R region, except in the G region. Moderately significant correlation coefficients (around +0.5) were found mainly for the NIR wavelengths (Figure 8a). For PC2, there were no significant correlations (−0.2 ≤ r ≤ +0.2) across almost all wavelengths, except very weak correlations (±0.2 < r ≤ ±0.3) with some RE and NIR wavelengths (Figure 8b). Overall, the strongest relationship was revealed with PC2 using the RE 722 nm wavelength (−0.6), where images from 17 March (GS55, heading) and 28 March (GS65, anthesis) revealed the highest loadings of 37% and 33%, respectively.
With GPC, PC1 as well as PC2 correlated weakly (±0.2 < r ≤ ±0.25) using some RE and all NIR wavelengths, with the 729 nm wavelength from PC1 having the highest coefficient (−0.25) and its surrounding wavelengths (722, 736 and 744 nm) slightly different coefficients. Images from 14 and 27 February and 11 March had the highest loadings (17%, 15% and 16%, respectively). Both PCs have in common that no correlations existed across all wavelengths from the V to the R region (Figure 8c,d).
In general, both PCs and their correlations with GY and GPC showed a smooth behavior across the available spectrum wavelengths, where peaks of correlation coefficients were shown at B, R and NIR spectral regions on PC1 for GY, with similar, although non-significant, inverted behavior for GPC. Although PC2 exhibited some noisy wavelengths within the G and RE regions for GPC, it still showed a continuous pattern across the available spectrum.

3.3.2. Integration of VIs Over Time

Another approach for performing a temporal exploratory analysis of a hyperspectral image time series is the integration of VIs over time through the calculation of the area under the curve (AUC). The temporal profiles for the lowest, median and highest GY and GPC within-field values are shown in Figure 9.
For GY, all VI temporal profiles revealed that the curves of highest, median and lowest GY followed a very similar course over the crop season with (a) increasing VI values during early growth stages (steam elongation to pre-booting; GS31–GS41; 14–27 February), (b) peaks at booting (GS41; 11 March) and anthesis (GS65; 28 March), (c) a slight depression at heading (GS55; 17 March), and d) a strong descending response starting from late milk (GS77; 15 April) until grain hard (GS92; 7 May). Furthermore, the highest GY values always showed higher VI values and vice versa, which is expected due to their plant structural information and its correlation to biomass and, consequently, to GY [91,92,93].
In contrast to GY, the VI temporal profiles for GPC show a different pattern over time for each VI. With exception of the median GPC curve, the highest and lowest GPC curves followed a similar course, whereby—starting from pre-booting crop stage (GS41; 27 February)—the lowest GPC values are always related to higher VI values. The low range variation of GPC (refer to Section 3.1) may be the reason of the ‘noisy’ effect in the VI response of the median GPC. The peaks vary: For PSNDc and PSSRc at anthesis (GS65; 28 March), for MCARI/OSAVI at physiological maturity (GS87; 25 April), and for TCARI at booting (GS41; 11 March) and late milk (GS77; 15 April). The VIs MCARI/OSAVI and TCARI indicate a higher potential for GPC estimation at late season, PSNDc and PSSRc at mid-season. Pigment-related VIs have shown to be negatively related to chlorophyll at specific crop stages [27], which is enough to explain the high VI responses to the lowest GPC value.
To understand the prediction potential of multi-temporal hyperspectral image analysis for GY and GPC, the relationships between the AUC of selected VI temporal profiles and the values of both GY and GPC were statistically evaluated at different crop stages (AUC1 to AUC8; Table 3). Table 6 summarizes the R2 resulting from fitting values between AUCs from selected VI profiles with GY and GPC.
Overall, independently of the crop stages (AUC1 to AUC8), the relationships across all VI temporal profiles were weak for GY (R2 ≤ 0.32) and very weak for GPC (R2 ≤ 0.1).
However, when focusing on GY, the best R2 were found for AUC1 using EVI and OSAVI (R2: 0.32), followed by AUC4 (R2: 0.28; PSSRc), AUC7 (R2: 0.27; TCARI), and AUC3 (R2: 0.26; PSSRc). This means that the prediction potential improved slightly (by just 4–6%) using 10 hyperspectral mosaics over the whole crop cycle (AUC1) compared to using 2 mosaics between heading and anthesis (AUC4) or 3 mosaics between grain filling and physiological maturity (AUC7) or booting and anthesis (AUC3). As such, this may not justify making such an investment for highly temporal time series of images.
The best results for GPC were obtained with the temporal VI profiles of TCARI and MCARI/OSAVI for the following AUCs: AUC5 (R2: 0.1; TCARI, MCARI/OSAVI), AUC2 (R2: 0.09; TCARI) and AUC3 (R2: 0.08; MCARI/OSAVI). This indicates that the most interesting period is from early to midseason, respectively from steam elongation (GS31) to anthesis (GS65). This is also the period when most N uptake (about 60%) takes place and the period of major photosynthetic capacity [100], which influences canopy expansion and, consequently, the green area index [97]. This same AUC approach was used for all NDSI combinations (data not shown). The results did not show any improvement in comparison with the single mosaic NDSI approach (spectral exploratory analysis section).
Other authors have made use of remote sensing data acquired throughout the season to develop spectral growth profiles based on VIs. Using NDVI as an indicator, Dubey et al. [101] found that the area under the growth profile explained nearly 69% (R2 around 0.47) of GY variability in wheat production districts in Indian states. Similarly, Xue et al. [49] analyzed cumulative VIs for diverse crop growth periods (jointing to maturity, booting to maturity, heading to maturity) to estimate GPC. However, in our case, the use of accumulated VIs showed no significant improvement in the correlations with GPC compared to single VIs.
In summary, multi-temporal spectral exploratory analysis showed its potential for identifying the optimal image timing across the crop cycle and reflectance wavelengths. Temporal PCA revealed that B, R and NIR spectral information at booting and heading crop stages carried the highest variance of the dataset, followed by late season images across almost all wavelengths, except NIR. Moreover, this approach showed the smooth behavior of the spectrum across the crop cycle, where peaks of correlation coefficients occurred at B, R and NIR spectral regions on PC1 for GY, with similar inverted behavior for GPC. Independently of the crop stages, the AUC approach based on VI temporal profiles demonstrated for all profiles weak relationships with GY (R2 ≤ 0.32) and very weak with GPC (R2 ≤ 0.1), whereby the most suitable VIs and the best temporal window for GPC assessment were identified.

3.4. Grain Protein Content Estimation Maps

To address the main focus of this study—the potential for timely GPC assessment for crop management as well as selective harvesting—and because grain protein comes at a higher aggregate value to farmers than yield, only the estimation of GPC is targeted in this section. From all the previously presented mono/multi-temporal spectral exploratory analysis methods, the result with the best fitting for GPC estimation was chosen from each method (Figure 10).
For the spectral exploratory analysis approach, the highest correlations with GPC were achieved using NDSI of a 700 and 574 nm combination from the 27 February image (R2: 0.19; RMSE: 0.46, root mean square error) and a 707 and 523 nm combination from the 7 May image (R2: 0.21; RMSE: 0.45). In the following figure, both are presented because the 27 February image provides information around the beginning of booting (early season), which is a suitable temporal window for crop management intervention targeting GPC correction. The 7 May image, which comes after physiological maturity (grain hard, late season), offers potential for selective harvesting strategies.
Within the multi-temporal exploratory analysis approach, the tested methods were tPCA and AUC of selected VIs. For the tPCA method, the 729 nm wavelength within PC1 obtained the highest correlation with GPC. The linear regression adjustment, although statistically significant at p-value < 0.05, showed an R2 of 0.06 and an RMSE of 0.49. For the AUC method, the TCARI index and MCARI/OSAVI ratio index of AUC5, which represents the area under the curve for the indices response from steam elongation (GS31) to anthesis (GS65), revealed the highest relationships with GPC (R2: 0.10; Table 6). Due to its number of used wavelengths and simplicity, just AUC5-TCARI (R2: 0.10; RMSE: 0.49) was selected for further presentation.
The comparison of relationships between GPC and the best fit of each exploratory analysis methods revealed (a) a general low to very low coefficient of determination (R2), (b) a general relatively high overall prediction error (RMSE) considering the low range variability of GPC for this field, and (c) that the spectral exploratory analysis approach surpasses the multi-temporal approach.
As outlined in the introduction, most remote and proximal sensing studies for GPC estimation were conducted on (a) fields and/or regions with naturally higher GPC variation [3,45,46] and/or (b) fields with different level treatments inducing higher GPC variation [2,3,27,47,48,49,50]. Generally, a higher range of calibration parameters allows better model adjustment for more precise prediction results. As an example, when studying relatively high GPC prediction accuracy, Wang et al. [46] demonstrated at the regional scale the potential of multi-temporal spectral image analysis using fused multispectral, broad-band satellite data for the crop growth period from jointing to anthesis (R2: 0.64). Although the overall prediction error of our study (RMSE: 0.45–0.49%) was relatively high for the conditions in the study area, the selected methods seem to be relatively more precise when compared to other proximal and remote sensing studies ([47], RMSE: 0.4–0.79%; [48], RMSE: 0.89–0.96%; [46], RMSE: 1.28–2.86%).
Figure 11 depicts three GPC maps based on: (a) measured GPC data resulting from in-situ samples and laboratory analysis, (b) estimated GPC data resulting from NDSI 700/574 nm spectral combination of the 27 February image, (c) estimated GPC data resulting from NDSI 707/523 nm spectral combination of 7 May, and their respective residual maps (difference between the estimated map and the interpolated map using ground measurements).
For the detection of lower (<12.5%) and higher (≥12.5%) wheat grain quality areas within the field, both estimated maps show a very similar size and distribution of class zones, whereby higher GPC areas were estimated based on the late-season image. The late-season GPC estimation map reflects a more realistic situation of GPC at harvest time (as illustrated by the kriged GPC map based on in-situ measured data).
In general, the prediction quality parameters (R2, RMSE) indicated an overall low prediction potential, but without any specific information regarding the within-field variability. In contrast, the residual maps of Block B (Figure 11) give a spatially detailed picture of the prediction potential, making it possible to visualize the within-field variation of GPC. Using the NDSI 700/574 nm spectral combination of the 27 February image, 61.7% of the total field area could be mapped with relatively good prediction (residual class: −0.25–+0.25%) and an even larger field area (74.9% of total field area) using the NDSI 707/523 nm spectral combination from the 7 May image. Overestimations (residual class <−0.25%: 19.8% for February 27, 13.6% for 7 May) and underestimations (residual class >+0.25%: 18.5% for February 27, 11.4% for 7 May) were very similar for both NDSI combinations and their corresponding crop stage.

4. Conclusions

This research presents the evaluation of spectral and multi-temporal spectral exploratory analysis methods for the assessment of within-field variations of wheat GY and GPC using high-resolution hyperspectral airborne imagery.
For both parameters, all three tested exploratory analysis methods—(a) narrow-band physiological spectral indices (spectral approach), (b) temporal principal component analysis (multi-temporal approach), and (c) the integration of VIs over time (multi-temporal approach)—made it possible to identify the most valuable reflectance wavelengths matching wheat physiology as well as to detect the optimal time window across the crop cycle for image acquisition based on the most relevant crop stages for grain protein prediction. Moreover, this study revealed a generally low to very low coefficient of determination, whereby the spectral exploratory analysis approach surpasses the multi-temporal approaches. Considering the low range of GPC at the study site, the overall prediction error for GPC estimation is relatively high. However, compared to other proximal and remote sensing-based GPC assessment studies, the RMSE values from this study seem to be relatively low. Moreover, residual maps derived during the regression analysis of NDSIs and GPC proved to be good for within-field prediction assessment and that, in our study, an area of ~75% of the total field could be mapped with relatively good prediction.
For future studies, this study delivers valuable information on the use of remote sensing data for the assessment of within-field GY and GPC variability aiming N management and selective harvesting: (a) detecting the most important crop stage for image acquisition and (b) identifying the best reflectance wavelengths along the crop cycle. New studies targeting N management for grain protein should be conducted, specifically using long-term experiments with different varieties and locations and applying different amounts of N at different times to generate calibrations with RS data. Furthermore, investigating GPC variation over different scales using different proximal and remote sensing platforms and corresponding spectro-radiometers (e.g., hand-held and on-combine sensors, UAV, aircraft, satellite-based sensors) may enable transferring GPC prediction models from the field-scale to the regional scale.

Author Contributions

Conceptualization, F.A.R., Jr.; Data curation, F.A.R., Jr.; Formal analysis, F.A.R., Jr.; Investigation, F.A.R., Jr., P.D., J.I.O.-M., U.S., P.J.Z.-T., J.A.T. and B.G.; Methodology, F.A.R., Jr., G.B., P.D. and P.J.Z.-T.; Validation, F.A.R., Jr., P.D., J.I.O.-M., J.A.T. and B.G.; Visualization, G.B., P.D., J.I.O.-M., U.S., J.A.T. and B.G.; Writing—original draft, F.A.R., Jr. and G.B.; Writing—review & editing, F.A.R., Jr., G.B., P.D., J.I.O.-M., U.S., P.J.Z.-T., J.A.T. and B.G.

Funding

This work was funded partially by the CGIAR Research Program on Wheat (www.wheat.org) and by the Spurring Transformation in Agriculture Research (STARS) project—under number 1094229-2014 (www.starsproject.org). Blasch and Taylor were supported by United Kingdom Government funding through the Newton Fund for this work.

Acknowledgments

We are very thankful to Jose Alberto Mendoza, Jose Urrea, Lorena González Pérez and the whole Crop Nutrition team for their assistance with field measurements, flight campaign and image processing; to CIMMYT’s Biometrics and Statistics Unit (BSU), for their support on statistical programming; to Dimitris Stratoulias for his valuable advice on hyperspectral data analysis; to Gemma Molero for her valuable contribution to the wheat physiology discussion and to Carlos Guzman for providing the methodology applied in the GPC laboratory analysis.

Conflicts of Interest

The authors declare no conflict of interest.

References

  1. FAO. FAOSTAT. 2017. Available online: http://www.fao.org/faostat/en/#data/QC (accessed on January 2018).
  2. Wright, D.; Rasmussen, V.; Ramsey, R.; Baker, D.; Ellsworth, J. Canopy reflectance estimation of wheat nitrogen content for grain protein management. GISci. Remote Sens. 2004, 41, 287–300. [Google Scholar] [CrossRef]
  3. Zhao, C.; Liu, L.; Wang, J.; Huang, W.; Song, X.; Li, C. Predicting grain protein content of winter wheat using remote sensing data based on nitrogen status and water stress. Int. J. Appl. Earth Obs. Geoinf. 2005, 7, 1–9. [Google Scholar] [CrossRef]
  4. Stewart, C.M.; Mcbratney, A.B.; Skerritt, J.H. Site-Specific durum wheat quality and its relationship to soil properties in a single field in Northern New South Wales. Precis. Agric. 2002, 3, 155–168. [Google Scholar] [CrossRef]
  5. Long, D.S.; Engel, R.E.; Carpenter, F.M. On-Combine Sensing and Mapping of Wheat Protein Concentration. Crop Manag. 2005, 4. [Google Scholar] [CrossRef]
  6. Long, D.S.; Engel, R.E.; Siemens, M.C. Measuring grain protein concentration with in-line near infrared reflectance spectroscopy. Agron. J. 2008, 100, 247–252. [Google Scholar] [CrossRef]
  7. Bramley, R.; Mowat, D.; Gobbett, D.; Branson, M.; Wakefield, A.; Wilksch, R. Mixing grapes and grain-Scoping the opportunity for selective harvesting in cereals. In Capturing Opportunities and Overcoming Obstacles in Australian Agronomy. Proceedings of 16th Australian Agronomy Conference 2012, Armidale, NSW. Australian, 14–18 October 2012; Australian Society of Agronomy Inc.: Armidale, Australia, 2012; Volume 16. [Google Scholar]
  8. Shewry, P.R. Improving the protein content and composition of cereal grain. J. Cereal Sci. 2007, 46, 239–250. [Google Scholar] [CrossRef]
  9. McMullan, P.M.; McVetty, P.B.E.; Urquhart, A.A. Dry matter and nitrogen accumulation and redistribution and their relationship to grain yield and grain protein in wheat. Can. J. Plant Sci. 1988, 68, 311–322. [Google Scholar] [CrossRef]
  10. Mariotti, F.; Tomé, D.; Mirand, P.P. Converting nitrogen into protein—Beyond 6.25 and Jones’ factors. Crit. Rev. Food Sci. Nutr. 2008, 48, 177–184. [Google Scholar] [CrossRef] [PubMed]
  11. Daniel, C.; Triboi, E. Effects of temperature and nitrogen nutrition on the grain composition of winter wheat: Effects on gliadin content and composition. J. Cereal Sci. 2000, 32, 45–56. [Google Scholar] [CrossRef]
  12. Ottman, M.J.; Doerge, T.A.; Martin, E.C. Durum grain quality as affected by nitrogen fertilization near anthesis and irrigation during grain fill. Agron. J. 2000, 92, 1035–1041. [Google Scholar] [CrossRef]
  13. Rharrabti, Y.; Villegas, D.; Garcia Del Moral, L.F.; Aparicio, N.; Elhani, S.; Royo, C. Environmental and genetic determination of protein content and grain yield in durum wheat under Mediterranean conditions. Plant Breed. 2001, 120, 381–388. [Google Scholar] [CrossRef]
  14. Altenbach, S.B.; DuPont, F.M.; Kothari, K.M.; Chan, R.; Johnson, E.L.; Lieu, D. Temperature, water and fertilizer influence the timing of key events during grain development in a US spring wheat. J. Cereal Sci. 2003, 37, 9–20. [Google Scholar] [CrossRef]
  15. Triboï, E.; Martre, P.; Triboï Blondel, A.M. Environmentally-induced changes in protein composition in developing grains of wheat are related to changes in total protein content. J. Exp. Bot. 2003, 54, 1731–1742. [Google Scholar] [CrossRef] [PubMed] [Green Version]
  16. Jamieson, P.D.; Semenov, M.A. Modelling nitrogen uptake and redistribution in wheat. Field Crop. Res. 2000, 68, 21–29. [Google Scholar] [CrossRef]
  17. Martre, P.; Porter, J.R.; Jamieson, P.D..; Triboï, E. Modeling grain nitrogen accumulation and protein composition to understand the sink/source regulations of nitrogen remobilization for wheat. Plant Phisiol. 2003, 133, 1959–1967. [Google Scholar] [CrossRef] [PubMed]
  18. Dupont, F.; Altenbach, S. Molecular and biochemical impacts of environmental factors on wheat grain development and protein synthesis. J. Cereal Sci. 2003, 38, 133–146. [Google Scholar] [CrossRef]
  19. Blankenau, K.; Olfs, H.W.; Kuhlmann, H. Strategies to improve the use efficiency of mineral fertilizer nitrogen applied to winter wheat. J. Agron. Crop Sci. 2002, 188, 146–154. [Google Scholar] [CrossRef]
  20. Fischer, R.A.; Howe, G.N.; Ibrahim, Z. Irrigated spring wheat and timing and amount of nitrogen fertilizer. I. Grain yield and protein content. Field Crop. Res. 1993, 33, 37–56. [Google Scholar] [CrossRef]
  21. Bogard, M.; Allard, V.; Brancourt-Hulmel, M.; Heumez, E.; MacHet, J.M.; Jeuffroy, M.H.; Gate, P.; Martre, P.; Le Gouis, J. Deviation from the grain protein concentration-grain yield negative relationship is highly correlated to post-anthesis N uptake in winter wheat. J. Exp. Bot. 2010, 61, 4303–4312. [Google Scholar] [CrossRef] [PubMed]
  22. Hirel, B.; Le Gouis, J.; Ney, B.; Gallais, A. The challenge of improving nitrogen use efficiency in crop plants: Towards a more central role for genetic variability and quantitative genetics within integrated approaches. J. Exp. Bot. 2007, 58, 2369–2387. [Google Scholar] [CrossRef] [PubMed]
  23. Gaju, O.; Allard, V.; Martre, P.; Snape, J.W.; Heumez, E.; LeGouis, J.; Moreau, D.; Bogard, M.; Griffiths, S.; Orford, S.; et al. Identification of traits to improve the nitrogen-use efficiency of wheat genotypes. Field Crop. Res. 2011, 123, 139–152. [Google Scholar] [CrossRef]
  24. Marti, J.; Bort, J.; Slafer, G.A.; Araus, J.L. Can wheat yield be assessed by early measurements of Normalized Difference Vegetation Index? Ann. Appl. Biol. 2007, 150, 253–257. [Google Scholar] [CrossRef]
  25. Gaju, O.; Allard, V.; Martre, P.; Le Gouis, J.; Moreau, D.; Bogard, M.; Hubbart, S.; Foulkes, M.J. Nitrogen partitioning and remobilization in relation to leaf senescence, grain yield and grain nitrogen concentration in wheat cultivars. Field Crop. Res. 2014, 155, 213–223. [Google Scholar] [CrossRef]
  26. Reyniers, M.; Vrindts, E. Measuring wheat nitrogen status from space and ground-based platform. Int. J. Remote Sens. 2006, 27, 549–567. [Google Scholar] [CrossRef]
  27. Wang, Z.J.; Wang, J.H.; Liu, L.Y.; Huang, W.J.; Zhao, C.J.; Wang, C.Z. Prediction of grain protein content in winter wheat (Triticum aestivum L.) using plant pigment ratio (PPR). Field Crop. Res. 2004, 90, 311–321. [Google Scholar] [CrossRef]
  28. Baret, F.; Guyot, G. Potentials and limits of vegetation indices for LAI and APAR assessment. Remote Sens. Environ. 1991, 35, 161–173. [Google Scholar] [CrossRef]
  29. Bannari, A.; Morin, D.; Bonn, F.; Huete, A.R. A review of vegetation indices. Remote Sens. Rev. 1995, 13, 95–120. [Google Scholar] [CrossRef]
  30. Merriman, J. Remote Sensing and Hyperspectral Data to Estimate Wheat and Maize Crop Characteristics in the Yaqui Valley, Mexico. Master’s Thesis, Université Catholique de Louvain, Louvain-la-Neuve, Belgium, 2017. [Google Scholar]
  31. Daughtry, C.S.T.; Walthall, C.L.; Kim, M.S.; Brown de Colstoun, E.; McMurtrey, J.E., III. Estimating corn leaf chlorophyll concentration from leaf and canopy reflectance. Remote Sens. Environ. 2000, 74, 229–239. [Google Scholar] [CrossRef]
  32. Blackburn, G.A. Spectral indices for estimating photosynthetic pigment concentrations: A test using senescent tree leaves. Int. J. Remote Sens. 1998, 19, 657–675. [Google Scholar] [CrossRef]
  33. Reyns, P.; Spaepen, P.; De Baerdemaeker, J. Site-specific relationship between grain quality and yield. Precis. Agric. 2000, 2, 231–246. [Google Scholar] [CrossRef]
  34. Delin, S. Within-field variations in grain protein content—Relationships to yield and soil nitrogen and consistency in maps between years. Precis. Agric. 2004, 5, 565–577. [Google Scholar] [CrossRef]
  35. Diacono, M.; Castrignanò, A.; Troccoli, A.; De Benedetto, D.; Basso, B.; Rubino, P. Spatial and temporal variability of wheat grain yield and quality in a Mediterranean environment: A multivariate geostatistical approach. Field Crop. Res. 2012, 131, 49–62. [Google Scholar] [CrossRef]
  36. Shanahan, J.F.; Schepers, J.S.; Francis, D.D.; Varvel, G.E.; Wilhelm, W.W.; Tringe, J.M.; Schlemmer, M.R.; Major, D.J. Use of remote-sensing imagery to estimate corn grain yield. Agron. J. 2001, 93, 583–589. [Google Scholar] [CrossRef]
  37. Leon, C.T.; Shaw, D.R.; Cox, M.S.; Abshire, M.J.; Ward, B.; Wardlaw, M.C.; Watson, C. Utility of Remote Sensing in Predicting Crop and Soil Characteristics. Precis. Agric. 2003, 4, 359–384. [Google Scholar] [CrossRef]
  38. Inman, D.; Khosla, R.; Reich, R.; Westfall, D.G. Normalized difference vegetation index and soil color-based management zones in irrigated maize. Agron. J. 2008, 100, 60–66. [Google Scholar] [CrossRef]
  39. Blaes, X.; Chomé, G.; Lambert, M.J.; Traoré, P.S.; Schut, A.G.T.; Defourny, P. Quantifying fertilizer application response variability with VHR satellite NDVI time series in a rainfed smallholder cropping system of Mali. Remote Sens. 2016, 8, 531. [Google Scholar] [CrossRef]
  40. Blaes, X.; Lambert, M.J.; Chomé, G.; Traore, P.S.; De By, R.A.; Defourny, P. Yield mapping for different crops in Sudano-Sahelian smallholder farming systems: Results based on metric Worldview and decametric SPOT-5 Take5 time series. In ESA Living Planet: Proceedings of ESA Living Planet Symposium 2016, Prague, Czech Republic, 9–13 May 2016; ESA: Paris, France, 2016; Volume SP-740. [Google Scholar]
  41. Rouse, J.W.; Hass, R.H.; Schell, J.A.; Deering, D.W.; Harlan, J.C. Monitoring the Vernal Advancement of Retrogradation of Natural Vegetation; Type III, Final Report; NASA’s Goddard Space Flight Center: Greenbelt, MD, USA, 1974; p. 371.
  42. Zarco-Tejada, P.J.; Ustin, S.L.; Whiting, M.L. Temporal and spatial relationships between within-field yield variability in cotton and high-spatial hyperspectral remote sensing imagery. Agron. J. 2005, 97, 641–653. [Google Scholar] [CrossRef]
  43. Eitel, J.U.H.; Long, D.S.; Gessler, P.E.; Hunt, E.R. Combined spectral index to improve ground-based estimates of nitrogen status in dryland wheat. Agron. J. 2008, 100, 1694–1702. [Google Scholar] [CrossRef]
  44. Basnet, B.B.; Apan, A.A.; Kelly, R.M.; Jensen, T.A.; Strong, W.M.; Butler, D.G. Relating satellite imagery with grain protein content. In Proceedings of the Spatial Sciences 2003 Conference, Canberra, Australia, 22–26 September 2003; pp. 1–11. [Google Scholar]
  45. Feng, M.C.; Xiao, L.J.; Zhang, M.J.; Ding, G.W. Integrating remote sensing and GIS for prediction of winter wheat (Triticum aestivum) protein contents in Linfen (Shanxi), China. PLoS ONE 2014, 9, e80989. [Google Scholar] [CrossRef] [PubMed]
  46. Wang, L.; Tian, Y.; Yao, X.; Zhu, Y.; Cao, W. Predicting grain yield and protein content in wheat by fusing multi-sensor and multi-temporal remote-sensing images. Field Crop. Res. 2014, 164, 178–188. [Google Scholar] [CrossRef]
  47. Hansen, P.M.; Jørgensen, J.R.; Thomsen, A. Predicting grain yield and protein content in winter wheat and spring barley using repeated canopy reflectance measurements and partial least squares regression. J. Agric. Sci. 2002, 139, 307–318. [Google Scholar] [CrossRef]
  48. Jensen, T.; Apan, A.; Young, F.; Zeller, L. Detecting the attributes of a wheat crop using digital imagery acquired from a low-altitude platform. Comput. Electron. Agric. 2007, 59, 66–77. [Google Scholar] [CrossRef] [Green Version]
  49. Xue, L.H.; Cao, W.X.; Yang, L.Z. Predicting grain yield and protein content in winter wheat at different N supply levels using canopy reflectance spectra. Pedosphere 2007, 17, 646–653. [Google Scholar] [CrossRef]
  50. Perry, E.M.; Fitzgerald, G.J.; Nuttall, J.G.; O’Leary, G.J.; Schulthess, U.; Whitlock, A. Rapid estimation of canopy nitrogen of cereal crops at paddock scale using a Canopy Chlorophyll Content Index. Field Crop. Res. 2012, 134, 158–164. [Google Scholar] [CrossRef]
  51. Becker-Reshef, I.; Vermote, E.; Lindeman, M.; Justice, C. A generalized regression-based model for forecasting winter wheat yields in Kansas and Ukraine using MODIS data. Remote Sens. Environ. 2010, 114, 1312–1323. [Google Scholar] [CrossRef]
  52. Dempewolf, J.; Adusei, B.; Becker-Reshef, I.; Hansen, M.; Potapov, P.; Khan, A.; Barker, B. Wheat yield forecasting for Punjab Province from vegetation index time series and historic crop statistics. Remote Sens. 2014, 6, 9653–9675. [Google Scholar] [CrossRef]
  53. Willis, P.R.; Carter, P.G.; Johannsen, C.J. Assessing yield parameters by remote sensing techniques. In Precision Agriculture; Robert, P.C., Rust, R.H., Larson, W.E., Eds.; American Society of Agronomy, Crop Science Society of America, Soil Science Society of America: Madison, WI, USA, 1999; pp. 1465–1473. [Google Scholar]
  54. Inoue, Y.; Miah, G.; Sakaiya, E.; Nakano, K.; Kawamura, K. NDSI Map and IPLS Using Hyperspectral Data for Assessment of Plant and Ecosystem Variables—With a Case Study on Remote Sensing of Grain Protein Content, Chlorophyll Content and Biomass in Rice. J. Remote Sens. Soc. Jpn. 2008, 28, 317–330. [Google Scholar]
  55. Inoue, Y.; Peñuelas, J.; Miyata, A.; Mano, M. Normalized difference spectral indices for estimating photosynthetic efficiency and capacity at a canopy scale derived from hyperspectral and CO2 flux measurements in rice. Remote Sens. Environ. 2008, 112, 156–172. [Google Scholar] [CrossRef]
  56. Stagakis, S.; Markos, N.; Sykioti, O.; Kyparissis, A. Monitoring canopy biophysical and biochemical parameters in ecosystem scale using satellite hyperspectral imagery: An application on a Phlomis fruticosa Mediterranean ecosystem using multiangular CHRIS/PROBA observations. Remote Sens. Environ. 2010, 114, 977–994. [Google Scholar] [CrossRef]
  57. Inoue, Y.; Sakaiya, E.; Zhu, Y.; Takahashi, W. Diagnostic mapping of canopy nitrogen content in rice based on hyperspectral measurements. Remote Sens. Environ. 2012, 126, 210–221. [Google Scholar] [CrossRef]
  58. Stratoulias, D.; Balzter, H.; Zlinszky, A.; Tóth, V.R. Assessment of ecophysiology of lake shore reed vegetation based on chlorophyll fluorescence, field spectroscopy and hyperspectral airborne imagery. Remote Sens. Environ. 2015, 157, 72–84. [Google Scholar] [CrossRef] [Green Version]
  59. Inoue, Y.; Guérif, M.; Baret, F.; Skidmore, A.; Gitelson, A.; Schlerf, M.; Darvishzadeh, R.; Olioso, A. Simple and robust methods for remote sensing of canopy chlorophyll content: A comparative analysis of hyperspectral data for different types of vegetation. Plant. Cell Environ. 2016, 39, 2609–2623. [Google Scholar] [CrossRef] [PubMed]
  60. Meisner, C.A.; Acevedo, E.; Flores, D.; Sayre, K.; Ortiz-Monasterio, J.I.; Byerlee, D. Wheat Production and Grower Practices in the Yaqui Valley, Sonora, Mexico; CIMMYT: Texcoco, Mexico, 1992; Volume 6, ISBN 9686127992. [Google Scholar]
  61. Zadoks, J.C.; Chang, T.T.; Konzak, C.F. A decimal code for the growth stages of cereals. Weed Res. 1974, 14, 415–421. [Google Scholar] [CrossRef]
  62. Gueymard, C. SMARTS2: A Simple Model of the Atmospheric Radiative Transfer of Sunshine: Algorithms and Performance Assessment; Florida Solar Energy Center: Cocoa, FL, USA, 1995. [Google Scholar]
  63. Gueymard, C.A. Interdisciplinary applications of a versatile spectral solar irradiance model: A review. Energy 2005, 30, 1551–1576. [Google Scholar] [CrossRef]
  64. Berni, J.A.J.; Zarco-Tejada, P.J.; Suárez, L.; Fereres, E.; Suarez, L.; Fereres, E. Thermal and narrowband multispectral remote sensing for vegetation monitoring from an unmanned aerial vehicle. IEEE Trans. Geosci. Remote Sens. 2009, 47, 722–738. [Google Scholar] [CrossRef] [Green Version]
  65. Zarco-Tejada, P.J.; González-Dugo, V.; Berni, J.A.J. Fluorescence, temperature and narrow-band indices acquired from a UAV platform for water stress detection using a micro-hyperspectral imager and a thermal camera. Remote Sens. Environ. 2012, 117, 322–337. [Google Scholar] [CrossRef] [Green Version]
  66. Zarco-Tejada, P.J.; González-Dugo, M.V.; Fereres, E. Seasonal stability of chlorophyll fluorescence quanti fi ed from airborne hyperspectral imagery as an indicator of net photosynthesis in the context of precision agriculture. Remote Sens. Environ. 2016, 179, 89–103. [Google Scholar] [CrossRef]
  67. Calderón, R.; Navas-Cortés, J.A.; Lucena, C.; Zarco-Tejada, P.J. High-resolution airborne hyperspectral and thermal imagery for early detection of Verticillium wilt of olive using fluorescence, temperature and narrow-band spectral indices. Remote Sens. Environ. 2013, 139, 231–245. [Google Scholar] [CrossRef] [Green Version]
  68. Calderón, R.; Navas-Cortés, J.A.; Zarco-Tejada, P.J. Early detection and quantification of Verticillium Wilt in olive using hyperspectral and thermal imagery over large areas. Remote Sens. 2015, 7, 5584–5610. [Google Scholar] [CrossRef]
  69. Pask, A.; Pietragalla, J.; Mullan, D. Physiological Breeding II: A Field Guide to Wheat Phenotyping; CIMMYT: Texcoco, Mexico, 2012; ISBN 9789706481825. [Google Scholar]
  70. Bramley, R.G.V.; White, R.E. An analysis of variability in the activity of nitrifiers in a soil under pasture. II. Some problems in the geostatistical analysis of biological soil properties. Aust. J. Soil Res. 1991, 29, 109–122. [Google Scholar] [CrossRef]
  71. Bramley, R.G.V. Understanding variability in winegrape production systems 2. Within vineyard variation in quality over several vintages. Aust. J. Grape Wine Res. 2005, 11, 33–42. [Google Scholar] [CrossRef]
  72. Isaaks, E.H.; Srivastava, R.M. Applied Geostatistics; Oxford University Press: New York, NY, USA, 1989. [Google Scholar]
  73. Rodrigues, F.A., Jr.; Ortiz-Monasterio, I.; Zarco-Tejada, P.J.; Schulthess, U.; Gérard, B. High resolution remote and proximal sensing to assess low and high yield areas in a wheat field. In Precision Agriculture 2015—Papers Presented at the 10th European Conference on Precision Agriculture, ECPA 2015; Wageningen Academic Publishers: Wageningen, The Netherlands, 2015. [Google Scholar]
  74. Webster, R.; Oliver, M.A. Geostatistics for Environmental Scientists; John Wiley & Sons: Hoboken, NJ, USA, 2007; ISBN 0470517263. [Google Scholar]
  75. AACC. Approved Methods of the AACC. In American Association of Cereal Chemists International; AACC: St. Paul, MN, USA, 2010; ISBN 978-1-891127-68-2. [Google Scholar]
  76. Minasny, B.; McBratney, A.B.; Whelan, B.M. VESPER Version 1.62; Australian Centre for Precision Agriculture: St. Louis, MO, USA, 2005. [Google Scholar]
  77. Feil, B. The inverse yield-protein relationship in cereals: Possibilities and limitations for genetically improving the grain protein yield. Trends Agron. 1997, 1, 103–119. [Google Scholar]
  78. Simmonds, N.W. The relation between yield and protein in cereal grain. J. Sci. Food Agric. 1995, 67, 309–315. [Google Scholar] [CrossRef]
  79. Slafer, G.A.; Andrade, F.H.; Feingold, S.E. Genetic improvement of bread wheat (Triticum aestivum L.) in Argentina: Relationships between nitrogen and dry matter. Euphytica 1990, 50, 63–71. [Google Scholar] [CrossRef]
  80. Whelan, B.M.; Taylor, J.A.; Hassall, J.A. Site-specific variation in wheat grain protein concentration and wheat grain yield measured on an Australian farm using harvester-mounted on-the-go sensors. Crop Pasture Sci. 2009, 60, 808–817. [Google Scholar] [CrossRef]
  81. Rodrigues, F.A.; Ortiz-Monasterio, I.; Zarco-Tejada, P.J.; Toledo, F.H.R.B. High resolution hyperspectral imagery to assess wheat grain protein in a farmer’s field. In Precision Agriculture 2016—Papers Presented at the 13th International Conference on Precision Agriculture, ECPA 2016; International Society of Precision Agriculture: Monticello, IL, USA, 2016; pp. 1–13. [Google Scholar]
  82. Huete, A.R.; Liu, H.Q.; Batchily, K.; Van Leeuwen, W. A comparison of vegetation indices over a global set of TM images for EOS-MODIS. Remote Sens. Environ. 1997, 59, 440–451. [Google Scholar] [CrossRef]
  83. Haboudane, D.; Miller, J.R.; Pattey, E.; Zarco-Tejada, P.J.; Strachan, I.B. Hyperspectral vegetation indices and novel algorithms for predicting green LAI of crop canopies: Modeling and validation in the context of precision agriculture. Remote Sens. Environ. 2004, 90, 337–352. [Google Scholar] [CrossRef]
  84. Rondeaux, G.; Steven, M.; Baret, F. Optimization of soil-adjusted vegetation indices. Remote Sens. Environ. 1996, 55, 95–107. [Google Scholar] [CrossRef]
  85. Haboudane, D.; Miller, J.R.; Tremblay, N.; Zarco-Tejada, P.J.; Dextraze, L. Integrated narrow-band vegetation indices for prediction of crop chlorophyll content for application to precision agriculture. Remote Sens. Environ. 2002, 81, 416–426. [Google Scholar] [CrossRef]
  86. Team, R.C. R: A Language and Environment for Statistical Computing; R Foundation for Statistical Computing: Vienna, Austria, 2016. [Google Scholar]
  87. Wickham, H. The split-apply-combine strategy for data analysis. J. Stat. Softw. 2011, 40, 1–29. [Google Scholar] [CrossRef]
  88. Wong, M.T.F.; Asseng, S.; Zhang, H. A flexible approach to managing variability in grain yield and nitrate leaching at within-field to farm scales. Precis. Agric. 2006, 7, 405–417. [Google Scholar] [CrossRef]
  89. Wiegand, C.L.; Richardson, A.J.; Escobar, D.E.; Gerbermann, A.H. Vegetation indices in crop assessments. Remote Sens. Environ. 1991, 35, 105–119. [Google Scholar] [CrossRef]
  90. Moriondo, M.; Maselli, F.; Bindi, M. A simple model of regional wheat yield based on NDVI data. Eur. J. Agron. 2007, 26, 266–274. [Google Scholar] [CrossRef]
  91. Zheng, T.C.; Zhang, X.K.; Yin, G.H.; Wang, L.N.; Han, Y.L.; Chen, L.; Huang, F.; Tang, J.W.; Xia, X.C.; He, Z.H. Genetic gains in grain yield, net photosynthesis and stomatal conductance achieved in Henan Province of China between 1981 and 2008. Field Crop. Res. 2011, 122, 225–233. [Google Scholar] [CrossRef]
  92. Xiao, Y.G.; Qian, Z.G.; Wu, K.; Liu, J.J.; Xia, X.C.; Ji, W.Q.; He, Z.H. Genetic gains in grain yield and physiological traits of winter wheat in Shandong province, China, from 1969 to 2006. Crop Sci. 2012, 52, 44–56. [Google Scholar] [CrossRef]
  93. Aisawi, K.A.B.; Reynolds, M.P.; Singh, R.P.; Foulkes, M.J. The physiological basis of the genetic progress in yield potential of CIMMYT spring wheat cultivars from 1966 to 2009. Crop Sci. 2015, 55, 1749–1764. [Google Scholar] [CrossRef]
  94. Thenkabail, P.S.; Smith, R.B.; De Pauw, E. Hyperspectral vegetation indices and their relationships with agricultural crop characteristics. Remote Sens. Environ. 2000, 71, 158–182. [Google Scholar] [CrossRef]
  95. Richardson, A.D.; Duigan, S.P.; Berlyn, G.P. An evaluation of noninvasive methods to estimate foliar chlorophyll content. New Phytol. 2002, 153, 185–194. [Google Scholar] [CrossRef] [Green Version]
  96. Sims, D.A.; Gamon, J.A. Relationships between leaf pigment content and spectral reflectance across a wide range of species, leaf structures and developmental stages. Remote Sens. Environ. 2002, 81, 337–354. [Google Scholar] [CrossRef]
  97. Sylvester-Bradley, R.; Berry, P.; Blake, J.; Kindred, D.; Spink, J.; Bingham, I.; McVittie, J.; Foulkes, J.; Edwards, C.; Dodgson, G. The Wheat Growth Guide; H.G.C.A: London, UK, 2008. [Google Scholar]
  98. Hurcom, S.J.; Harrison, A.R.; Taberner, M. Assessment of biophysical vegetation properties through spectral decomposition techniques. Remote Sens. Environ. 1996, 56, 203–214. [Google Scholar] [CrossRef]
  99. Kaiser, H.F. The application of electronic computers to factor analysis. Educ. Psychol. Meas. 1960, 20, 141–151. [Google Scholar] [CrossRef]
  100. Slafer, G. Wheat development: Its role in phenotyping and improving crop adaptation. In Physiological Breeding I: Interdisciplinary Approaches to Improve Crop Adaptation; Reynolds, M.P., Pask, A.J.D., Mullan, D.M., Eds.; CIMMYT: Texcoco, Mexico, 2012; ISBN 9706481818. [Google Scholar]
  101. Dubey, R.P.; Ajwani, N.D.; Navalgund, R.R. Relation of wheat yield with parameters derived from a spectral growth profile. J. Indian Soc. Remote Sens. 1991, 19, 27–44. [Google Scholar] [CrossRef]
Figure 1. Location of the study area in Mexico, showing a hyperspectral early-season colour infrared (CIR) image (February 14) of the study fields (Block A and B) and used ground control points (GPCs).
Figure 1. Location of the study area in Mexico, showing a hyperspectral early-season colour infrared (CIR) image (February 14) of the study fields (Block A and B) and used ground control points (GPCs).
Remotesensing 10 00930 g001
Figure 2. (a) GY and (b) GPC maps based on manually-collected grain samples at the two study sites (Block A: only point values, n = 14; Block B: spatially interpolated point values using block kriging, n = 100) (GPC-threshold = 12.5% according to the farmer’s recommendation).
Figure 2. (a) GY and (b) GPC maps based on manually-collected grain samples at the two study sites (Block A: only point values, n = 14; Block B: spatially interpolated point values using block kriging, n = 100) (GPC-threshold = 12.5% according to the farmer’s recommendation).
Remotesensing 10 00930 g002
Figure 3. Reflectance profile from the highest and lowest GY (a) and GPC (b) measured sampling points across blocks.
Figure 3. Reflectance profile from the highest and lowest GY (a) and GPC (b) measured sampling points across blocks.
Remotesensing 10 00930 g003
Figure 4. Contour maps of the coefficient of determination (R2) between NDSI (Ri, Rj) and GY using the complete combinations of two wavelengths at i and j nm.
Figure 4. Contour maps of the coefficient of determination (R2) between NDSI (Ri, Rj) and GY using the complete combinations of two wavelengths at i and j nm.
Remotesensing 10 00930 g004
Figure 5. Contour maps of the coefficient of determination (R2) between NDSI (Ri, Rj) and GPC using the complete combinations of two wavelengths at i and j nm.
Figure 5. Contour maps of the coefficient of determination (R2) between NDSI (Ri, Rj) and GPC using the complete combinations of two wavelengths at i and j nm.
Remotesensing 10 00930 g005
Figure 6. Eigenvalues/percentage of variance of each PC within each wavelength.
Figure 6. Eigenvalues/percentage of variance of each PC within each wavelength.
Remotesensing 10 00930 g006
Figure 7. Percent contribution of each image to the total variance, per wavelength, for PC1 (a) and PC2 (b).
Figure 7. Percent contribution of each image to the total variance, per wavelength, for PC1 (a) and PC2 (b).
Remotesensing 10 00930 g007
Figure 8. Correlation coefficients between each wavelength from PC1 with GY and GPC (a,c) and wavelengths from PC2 with GY and GPC (b,d). The dashed line represents a 5% significance level.
Figure 8. Correlation coefficients between each wavelength from PC1 with GY and GPC (a,c) and wavelengths from PC2 with GY and GPC (b,d). The dashed line represents a 5% significance level.
Remotesensing 10 00930 g008
Figure 9. Temporal profiles of selected vegetation indices (VIs) to the within-field highest, median and lowest measured GY and GPC.
Figure 9. Temporal profiles of selected vegetation indices (VIs) to the within-field highest, median and lowest measured GY and GPC.
Remotesensing 10 00930 g009
Figure 10. Regression plots between GPC and the best fittings of each exploratory approach, where p-value: p-values of the regression analysis; R2: coefficient of determination of the regression analysis; RMSE: root mean square error of the regression analysis.
Figure 10. Regression plots between GPC and the best fittings of each exploratory approach, where p-value: p-values of the regression analysis; R2: coefficient of determination of the regression analysis; RMSE: root mean square error of the regression analysis.
Remotesensing 10 00930 g010
Figure 11. Measured interpolated GPC map (as shown in Figure 1) (a); estimated early-season GPC map and derived residual map (b); estimated late-season GPC map and derived residual map (c). White spots in the estimated maps are weed spots treated by the farmer.
Figure 11. Measured interpolated GPC map (as shown in Figure 1) (a); estimated early-season GPC map and derived residual map (b); estimated late-season GPC map and derived residual map (c). White spots in the estimated maps are weed spots treated by the farmer.
Remotesensing 10 00930 g011
Table 1. Image dates and respective crop growth stages.
Table 1. Image dates and respective crop growth stages.
Image DateCrop Growth Stage
14 FebruaryInitiation of stem elongation (GS31)
19 FebruaryStem elongation period
27 FebruaryStem elongation period
11 MarchBooting (GS41)
17 MarchHeading (GS55)
28 MarchAnthesis (GS65)
7 AprilGrain filling (GS71)
15 AprilLate milk (GS77)
25 AprilPhysiological maturity (GS87)
7 MayGrain hard (GS92)
Table 2. Selected VIs for GY and GPC.
Table 2. Selected VIs for GY and GPC.
IndexFormulaReference
GY
Enhanced vegetation index (EVI) 2.5 ( R 800 R 770 R 800 + 6 R 670 7.5 R 400 + 1 ) [82]
Modified triangular vegetation index 2# (MTVI2) 1.2 [ 1.2 ( R 800 R 550 ) 2.5 ( R 670 R 550 ) ] ( 2 R 800 + 1 ) 2 ( 6 R 800 5 R 670 ) 0.5 [83]
Normalized difference vegetation index (NDVI) R 800 R 680 R 800 + R 680 [41]
Optimized soil-adjusted vegetation index (OSAVI) ( 1 + 0.16 ) ( R 800 R 670 ) ( R 800 + R 670 + 0.16 ) [84]
GPC
Ratio: Modified chlorophyll absorption ratio index/Optimized soil-adjusted vegetation index (MCARI/OSAVI) [ ( R 700 R 670 ) 0.2 ( R 700 R 550 ) ] R 700 R 670 ( 1 + 0.16 ) ( R 800 R 670 ) / ( R 800 + R 670 + 0.16 ) [31]
Pigment specific normalized difference c (PSNDc) R 800 R 470 R 800 + R 470 [32]
Pigment specific simple ratio for carotenoids (PSSRc) R 800 R 470 [32]
Transformed chlorophyll absorption in reflectance index (TCARI) 3 [ ( R 700 R 670 ) 0.2 ( R 700 R 550 ) ] R 700 R 670 [85]
Table 3. Range of crop growth stages used to calculate the area under the curve (AUC) for different vegetation indices.
Table 3. Range of crop growth stages used to calculate the area under the curve (AUC) for different vegetation indices.
CodeDescription of Crop Growth StageDates of Used ImagesNumber of Mosaics
AUC1Whole time series14 February to 7 May10
AUC2Steam elongation (GS31) to booting (GS41)14 February to 11 March4
AUC3Booting (GS41) to anthesis (GS65)11 March to 28 March3
AUC4Heading (GS55) to anthesis (GS65)17 March to 28 March2
AUC5Steam elongation (GS31) to anthesis (GS65)14 February to 28 March6
AUC6Grain filling (GS71) to late milk (GS77)7 April to 15 April2
AUC7Grain filling (GS71) to physiological maturity (GS87)7 April to 25 April3
AUC8Grain filling (GS71) to grain hard (GS92)7 April to 7 May4
Table 4. Descriptive statistics for GY and GPC from both blocks (n = 114).
Table 4. Descriptive statistics for GY and GPC from both blocks (n = 114).
GYGPC
Maximum8.0214.96
3º Quartile6.8412.56
Median6.3612.25
1º Quartile5.9812.03
Minimum4.6610.87
Mean6.4112.32
Skewness0.021.25
Kurtosis−0.145.67
CV11.514.15
GY—Mg ha−1; GPC—%.
Table 5. Descriptive statistics of Pearson’s correlation coefficients between GY and GPC obtained through the moving window process.
Table 5. Descriptive statistics of Pearson’s correlation coefficients between GY and GPC obtained through the moving window process.
Remotesensing 10 00930 i001 r CoefficientRespective p-ValueRespective Degrees of Freedom
Maximum0.610.077
3º Quartile0.080.809
Median−0.120.719
1º Quartile−0.260.3015
Minimum−0.760.045
Mean−0.100.359
Skewness0.10--
Kurtosis−0.04--
n100--
Table 6. R-squared of fitting between the areas under the curve (AUC) of the selected vegetation index profiles for different crop stages and GY and GPC.
Table 6. R-squared of fitting between the areas under the curve (AUC) of the selected vegetation index profiles for different crop stages and GY and GPC.
AUC1; #10AUC2; #4AUC3; #3AUC4; #2AUC5; #6AUC6; #2AUC7; #3AUC8; #4
GYGPCGYGPCGYGPCGYGPCGYGPCGYGPCGYGPCGYGPC
NDVI0.280.010.120.060.250.010.260.010.200.040.130.010.130.010.100.02
EVI0.320.020.150.080.250.020.250.010.220.060.22-0.220.010.180.01
OSAVI0.320.010.140.070.250.020.260.010.220.050.20-0.180.010.140.02
MTVI20.310.020.150.080.250.020.250.010.220.070.21-0.210.010.160.01
PSNDc0.280.010.110.050.250.010.260.010.190.040.14-0.150.010.120.02
MCARI/OSAVI0.080.060.060.060.010.08-0.070.050.10-0.070.090.070.03-
PSSRc0.240.020.100.060.260.020.280.010.190.050.13-0.140.010.130.01
TCARI0.240.070.100.090.110.070.090.050.130.100.120.040.270.010.16-
AUC1—Whole time series of images; AUC2—images from steam elongation (GS31) to booting (GS41); AUC3—images from booting (GS41) to anthesis (GS65); AUC4—images from heading (GS55) to anthesis (GS65); AUC5—images from stem elongation (GS31) to anthesis (GS65); AUC6—images from grain filling (GS71) to late milk (GS77); AUC7—images from grain filling (GS71) to physiological maturity (GS87); AUC8—images from grain filling (GS71) to grain hard (GS92). # indicates the number of mosaics used in each AUC. Bold r-squared are significant at p < 0.05.

Share and Cite

MDPI and ACS Style

Rodrigues, F.A., Jr.; Blasch, G.; Defourny, P.; Ortiz-Monasterio, J.I.; Schulthess, U.; Zarco-Tejada, P.J.; Taylor, J.A.; Gérard, B. Multi-Temporal and Spectral Analysis of High-Resolution Hyperspectral Airborne Imagery for Precision Agriculture: Assessment of Wheat Grain Yield and Grain Protein Content. Remote Sens. 2018, 10, 930. https://doi.org/10.3390/rs10060930

AMA Style

Rodrigues FA Jr., Blasch G, Defourny P, Ortiz-Monasterio JI, Schulthess U, Zarco-Tejada PJ, Taylor JA, Gérard B. Multi-Temporal and Spectral Analysis of High-Resolution Hyperspectral Airborne Imagery for Precision Agriculture: Assessment of Wheat Grain Yield and Grain Protein Content. Remote Sensing. 2018; 10(6):930. https://doi.org/10.3390/rs10060930

Chicago/Turabian Style

Rodrigues, Francelino A., Jr., Gerald Blasch, Pierre Defourny, J. Ivan Ortiz-Monasterio, Urs Schulthess, Pablo J. Zarco-Tejada, James A. Taylor, and Bruno Gérard. 2018. "Multi-Temporal and Spectral Analysis of High-Resolution Hyperspectral Airborne Imagery for Precision Agriculture: Assessment of Wheat Grain Yield and Grain Protein Content" Remote Sensing 10, no. 6: 930. https://doi.org/10.3390/rs10060930

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop