Towards an Improved High-Throughput Phenotyping Approach: Utilizing MLRA and Dimensionality Reduction Techniques for Transferring Hyperspectral Proximal-Based Model to Airborne Images

Abstract

At present, it is critical to accurately monitor wheat crops to help decision-making processes in precision agriculture. This research aims to retrieve various wheat crop traits from hyperspectral data using machine learning regression algorithms (MLRAs) and dimensionality reduction (DR) techniques. This experiment was conducted in an agricultural field in Arborea, Oristano-Sardinia, Italy, with different factors such as cultivars, N-treatments, and soil ploughing conditions. Hyperspectral data were acquired on the ground using a full-range Spectral Evolution spectrometer (350–2500 nm). Four DR techniques, including (i) variable influence on projection (VIP), (ii) principal component analysis (PCA), (iii) vegetation indices (VIs), and (iv) spectroscopic feature (SF) calculation, were undertaken to reduce the dimension of the hyperspectral data while maintaining the information content. We used five MLRA models, including (i) partial least squares regression (PLSR), (ii) random forest (RF), (iii) support vector regression (SVR), (iv) Gaussian process regression (GPR), and (v) neural network (NN), to retrieve wheat traits at either leaf and canopy levels. The studied traits were leaf area index (LAI), leaf and canopy water content (LWC and CWC), leaf and canopy chlorophyll content (LCC and CCC), and leaf and canopy nitrogen content (LNC and CNC). MLRA models were able to accurately retrieve wheat traits at the canopy level with PLSR and NN indicating the highest modelling performance. On the contrary, MLRA models indicated less accurate retrievals of the leaf-level traits. DR techniques were found to notably improve the retrieval accuracy of crop traits. Furthermore, the generated models were re-calibrated using soil spectra and then transferred to an airborne dataset collected using a CASI-SASI hyperspectral sensor, allowing the estimation of wheat traits across the entire field. The predicted crop trait maps illustrated consistent patterns while also preserving the real-field characteristics well. Lastly, a statistical paired t-test was undertaken to conduct a proof of concept of wheat phenotyping analysis considering the different agricultural variables across the study site. N-treatment caused significant differences in wheat crop traits in many instances, whereas the observed differences were less pronounced between the cultivars. No particular impact of soil ploughing conditions on wheat crop characteristics was found. Using such combinations of MLRA and DR techniques based on hyperspectral data can help to effectively monitor crop traits throughout the cropping seasons and can also be readily applied to other agricultural settings to help both precision farming applications and the implementation of high-throughput phenotyping solutions.

1. Introduction

Wheat is one of the leading cultivated crops across the globe [1]. Hence, accurate monitoring of wheat crop status over agricultural settings is essential to regulate the increasing trend of global food demand [2]. Traditional ground-based crop monitoring activities are often time-consuming, expensive, and destructive [3]. Whilst a new era of agricultural science is emerging owing to recent advances in remote and proximal sensing techniques. Hyperspectral sensors on both remote and proximal platforms can provide detailed spectral, spatial, and temporal information about crop status [4].

In the context of precision agriculture, prominent crop status monitoring relies on determining leaf biochemical properties and canopy structure [5]. In this regard, leaf area index (LAI), leaf and canopy water content (LWC and CWC), leaf and canopy chlorophyll content (LCC and CCC), and leaf and canopy nitrogen content (LNC and CNC) are key proxies for crop photosynthetic capacity, nutrient status, water use efficiency, and physiological state, e.g., [5,6,7]. One way to retrieve the aforementioned crop traits from remotely sensed data is the physically based radiative transfer modelling (RTM) [8]. RTMs simulate vegetation reflectance spectra based on leaf biophysical properties and canopy structure and then invert the simulated spectra, adopting different schemas to estimate crop traits [5,9]. In addition, crop traits can be retrieved using data-driven approaches through parametric or nonparametric regression methods [8]. Parametric regression simply links a derivate of spectra (such as vegetation indices) to the crop trait through a regression model [3]. Whilst, in non-parametric regression models, such as machine learning, the fitting model is not explicit and is defined using the method of empirical relationships [3,8].

In recent years, machine learning regression algorithms (MLRAs) have been successfully utilized to determine various crop properties based on remotely sensed data [6]. Until now, many previous studies have evaluated the potential of hyperspectral data in retrieving crop traits based on MLRA models, e.g., [8,10,11]. However, such high spectral dimensionality of hyperspectral data often leads to challenges associated with signal processing, spectral multicollinearity, computational capacity, and redundant information [12,13]. In this respect, dimensionality reduction (DR) is capable of transforming hyperspectral data into a lower-dimensional space while preserving the majority of the original spectral information [6].

To the best of our knowledge, no comprehensive study has yet been devoted to the performance of different MLRA models with respect to various DR techniques when retrieving wheat crop traits. Therefore, in this study, we have attempted to retrieve wheat traits at leaf and canopy levels, utilizing either full-range or DR spectral data as inputs to different MLRA models. The best-identified models for each crop trait were then transferred to the hyperspectral image acquired using a manned aircraft to perform wheat phenotyping analysis across the study field. Our methodology can help decision-makers with agricultural practices and, in principle, can also be applicable to other crop types.

2. Materials and Methods

2.1. Experimental Design

The study site (NW corner: 39°45′N, 8°35′E; SE corner: 39°45′N, 8°36′E) was an agricultural field in Arborea, Oristano—Sardinia, Itay, which was planted with durum wheat (Triticum aestivum L.) (Figure 1). To generate a wide variety of crop traits (allowing for retrieving a wide range of spectral information), various agricultural inputs were broadcasted over the field, as reported in Figure 2.

The field was composed of two soil ploughing characteristics: minimum tillage (L1) and traditional ploughing (L2) in the upper and lower portions of the field, respectively (Figure 2a). The experimental design constituted a block of three replicates: R1, R2, and R3, representing the western, middle, and eastern parts of the field. Each replicate included a fully randomized set of four different cultivars: Beltorax (V1), F. Camillo (V2), Giulio (V3), and M. Aurelio (V4) (Figure 2b), as well as four various N-treatments (i.e., ammonium nitrate) as Standard-30% (N1), Standard (N2), Standard + 30% (N3), and a variable rate according to the expert decision based on a prefertilization assessment of crop status with LAI and chlorophyll measurements (N4), totalling 96 experimental plots of 6 m × 3 m within the field.

2.2. Crop Trait Measurements

Non-destructive measurements of crop traits were performed in three campaigns on 28 April, 14 May, and 28 May 2022, as follows:

(i)

LAI: LAI was measured using an LI-COR LAI2200 plant analyser along transects according to an A-10 × B-A scheme, in which A and B represent above and below canopy measurements for each experimental plot. The average of 10 LAI2200 measurements made within each plot was used to describe the LAI in each experimental plot.

(ii)

LCC and CCC: LCC was measured using an MC-100 chlorophyll content meter based on the ratio of light transmittance through the plant leaves at 653 and 931 nm spectral bands. LCC measurements were conducted on the last fully developed leaf of 10 randomly selected plants within each experimental plot. CCC was then computed by multiplying LCC and the corresponding LAI values [9].

(iii)

LWC and CWC: For LWC determination, five plants for each experimental plot of V3 and V4 varieties were considered, totalling 24 samples. A set of three leaf discs with an 8 mm diameter was sampled on the last fully developed leaf of the selected plants using a handheld punch. The samples were then oven-dried at 50 °C for about three days until reaching a constant weight to retrieve LWC. This destructive sampling of LWC and LMA took place on 29 April, 13 May, and 30 May. CWC was retrieved based on multiplying LWC by LAI values [9].

(iv)

LNC and CNC: The aforementioned leaf discs were used to determine leaf nitrogen concentration (Nmass) via dry combustion using the CN element analyser. LNC was then calculated by multiplying Nmass and leaf mass area. Subsequently, CNC was retrieved by multiplying LNC and the corresponding LAI values [9].

As stated before in Section 2.1, 96 experimental plots (i.e., 3 replicates × 2 tillages × 4 varieties × 4 N strategies) were designed throughout the study field. It should be noted that 6 plots along the outer edge of area L1 were excluded from ground sampling because of the excessive weed presence that effectively stifled crop emergence. Moreover, 2 plots in the final part of area L2 were omitted because of the low crop density due to heavy rainfall and consequent water stagnation in the pre-emergence phase. As such, a total of 88 plots were monitored three times during the 2022 cropping season. Destructive foliar measurements were limited to a set of plots of interest relative to L2 tillage and 2 varieties (Giulio (V3), recently introduced to the market, and Marco Aurelio (V4), characterized using excellent quanti-qualitative performance) for a total of 24 plots (3 replicate × 1 tillage × 2 variety × 4 N strategy).

2.3. Spectral Data Acquisition

2.3.1. Proximal Sensing—Spectrometer

Proximal spectral measurements were acquired using a handheld spectrometer (Spectral Evolution) over the three field campaigns. Spectral Evolution (SE) has a full spectral range of 350–2500 nm with a 1 nm spectral sampling interval (SSI). Four spectral measurements were conducted per each experimental plot at approximately 1 m above the canopy with a nadir viewing angle and under clear sky conditions. Care was taken to avoid the presence of weeds and shadows within the measuring footprint. A reference panel was sensed before and after performing the four measurements of each plot to derive reflectance. The average of the four reflectance measurements derived for each experimental plot was used to describe the spectra per plot (Figure 3). A spline smoothing filter was then applied to the entire spectra to reduce the high noise levels in the retrieved spectra. It is worth noting that during the first campaign, the spectral measurements of 50 plots were not recorded due to instrumental issues.

2.3.2. Manned Airborne Remote Sensing—Aircraft

Another set of hyperspectral images was collected close to the last campaign on 29 May by CzechGlobe using a Cessna 208B Grand Caravan manned aircraft carrying CASI-1500 and SASI-600 hyperspectral sensors simultaneously. CASI-1500 and SASI-600 are pushbroom hyperspectral sensors with a spectral range (resp. SSI) of 380–1050 nm (3.2 nm) and 1000–2450 nm (15 nm), respectively. The aircraft flight was performed at 3090 m AGL with an airspeed of 56 m/s, producing hyperspectral images as two datasets with a ground sampling distance (GSD) of 0.25 m and 0.63 m for VNIR and SWIR data, respectively. To maintain a high signal-to-noise ratio for the CASI instrument, spectral binning has been performed in our experiment. Due to this, CASI acquired data with 48 bands and a full-width half maximum (FWHM) corresponding to 14.3 nm. Aerial data were provided by CzechGlobe with radiometric and geometric correction, then in radiance at the georeferenced sensor. It is worth mentioning that comparison with ground radiance, as measured using a spectroradiometer, on standard tarps (i.e., white, grey, and black) highlighted the need to recalibrate SASI data adopting an empirical line approach (details are provided in Appendix A).

More precisely, prior to the atmospheric correction, each CASI and SASI line has been co-registered in a single hypercube using the geometric information of the CASI data (i.e., oversampling the SASI images to 0.25 m GSD with nearest neighbour interpolation in ENVI 5.6.3^® layer-stacking procedure) to process them simultaneously in ATCOR 4 and make use of the full spectrum information. From the layer stacks, the last NIR CASI bands overlapping with the first NIR SASI bands have been removed, yielding a VNIR-SWIR sensor with 141 spectral bands between 384.6 and 2442.50 nm. This synthetic sensor was then defined in ATCOR with the appropriate characteristics. The sensor’s spectral filter functions were defined following a Butterworth order 2 filter, with no spectral binning, and using the FWHM of each band as defined in the data headers.

Average ground spectral signatures from the black and white tarps, acquired with the SE spectroradiometer, were spectrally resampled to the airborne sensor using the ATCOR4 “RESPECT” routine, which uses the spectral response files of the airborne sensor to resample ground measurements. The resampled signatures were then used to perform an inflight calibration using the ATCOR4 Inflight Calibration module. Average spectra from the tarps on the flight line were regressed against resampled ground reflectance spectra to yield the coefficient of a linear regression that was afterward applied to the airborne radiance images.

After the generation of the calibration file with the Inflight Module, optimal visibility and atmospheric parameters were determined via an iterative procedure in the ATCOR4 “SPECTRA” module using ground reflectance signatures over the experimental fields as reference. The SPECTRA procedure yielded the following optimal parameters for the atmospheric correction: Atmosphere Aerosol Type: Rural, Atmosphere Water Vapor Content: 2.0 g·cm⁻², Adjacency Range: 0.1 Km, Adjacency Zones: 1, Visibility: 43.9 Km.

Finally, the atmospheric correction was run using no interpolation across the whole spectrum and with a water vapor algorithm with band regression in the 820 nm region (Left window: 797.9–798.0 nm, Right window: 840.7–840.8 nm, Absorption region: 825.9–826.9 nm). Given the flatness of the area, no digital elevation model was used for the atmospheric correction procedure.

2.3.3. Spectra Pre-Processing

Due to the different spectral characteristics of the acquired hyperspectral datasets (Section 2.3.1 and Section 2.3.2), we first resampled all the spectra data into the PRecursore IperSpettrale della Missione Applicativa (PRISMA) spectral configuration [14]. The PRISMA hyperspectral sensor (launched by the Italian Space Agency) is a new generation sensor that provides data useful for developing solutions that exploit space spectroscopy data. PRISMA has a spectral range of 400–2500 nm with an average SSI of less than 10 nm, yielding a total of 231 spectral bands. The spectral resampling was performed through the prospectr package [15] in R (R Core Team, Foundation for Statistical Computing, Vienna, Austria), using Gaussian spectral response functions (SRF) generated according to centre wavelengths and full-width-half-maximum values of PRISMA bands. Moreover, atmospheric regions between 1350–1500 nm and 1750–2000 nm and the last portion of the SWIR between 2400–2500 nm were excluded [16], resulting in a total of 170 so-called PRISMA-like spectral bands. This approach allows us to adequately transfer the generated MLRA models based on proximal data to the airborne dataset.

2.4. DR Analysis

Various DR analysis was performed in the present study, including (i) VIPs, (ii) PCA, (iii) VIs, and (iv) SF, as follows:

Variable influence on projection (VIP) indicates the importance of the input variables (regressors) for predicting the response variable [17]. In other words, VIP computes the importance of the input variables (in multivariate models) based on projections to latent structures to identify VIP scores. Then, a so-called VIP graph is defined according to the computed VIP scores. Based on the VIP graphs, user-defined thresholds (VIP scores > 1 as suggested by Akarachantachote et al. [18]) were considered to select the most important variables (i.e., highest VIP scores) for crop trait retrieval in this study.

Principal component analysis (PCA) transforms hyperspectral data into a new reduced-dimension dataset in which the first component contains the highest information of the original data [19]. More precisely, PCA computes the covariance matrix of the original dataset and, subsequently, the corresponding eigenvectors and eigenvalues. The eigenvectors related to the largest eigenvalues of the covariance matrix identify the principal components. Then, the original dataset is projected to a new orthogonal space defined using these eigenvectors, in which the new uncorrelated variables successively maximize the variance. Different numbers of components can be considered in PCA technique for hyperspectral studies, e.g., [9,16]. In this study, only PCA with 20 components (PCA-20), as proposed by [20], was computed. The mathematical description of VIP and PCA methods is depicted in Appendix B.

Computing vegetation indices (VIs) is a traditional feature extraction method that extracts a limited number of bands from the original spectral dataset to provoke spectral features that are sensitive to a specific crop trait [13]. Based on the literature, we chose to focus on a total of 53 VIs (Appendix C;

Table A1. List of hyperspectral vegetation indices (VIs) investigated in this study.

VI	Formula	Reference
CARI	R700 × abs(a × 670 + R670 + b)/R670 × (a2 + 1)0.5 a = (R700 − R550)/150 b = R550 − (a × 550)	[37]
CARTER1	R695/R420	[38]
CARTER2	R695/R760	[38]
CARTER3	R605/R760	[38]
CARTER4	R710/R760	[38]
CARTER5	R695/R670	[38]
CIgreen	(R780/R550) − 1	[39]
CIrededge	(R780/R710) − 1	[39]
DATT1	(R850 − R710)/(R850 − R680)	[40]
DATT2	R850/R710	[40]
DVI	R800 − R680	[41]
EVI	2.5 × ((R800 − R670)/(R800 − (6 × R670) − (7.5 × R475) + 1))	[42]
GI	R554/R677	[43]
GNDVI	(R800 − R550)/(R800 + R550)	[44]
MCARI	((R700 − R670) − 0:2 × (R700 − R550)) × (R700/R670)	[45]
MCARI_d_OSAVI	MCARI/OSAVI	[45]
mNDVI705	(R750 − R705)/(R750 + R705 − 2 × R445)	[46]
MSAVI	0.5 × (2 × R800 + 1 − ((2 × R800 + 1)2 − 8 × (R800 − R670))0.5)	[47]
MTCI	(R754 − R709)/(R709 − R681)	[48]
NDCI	(R762 − R527)/(R762 + R527)	[49]
NDNI	(log(1/R1510) − log(1/R1680))/(log(1/R1510) + log(1/R1680))	[50]
NRI1510	(R1510 − R660)/(R1510 + R660)	[51]
NDRE	(R800 − R739)/(R800 + R739)	[52]
NDVI	(R800 − R680)/(R800 + R680)	[53]
NDVI705	(R750 − R705)/(R750 + R705)	[54]
NVI	(R777 − R747)/R673	[55]
OSAVI	(1 + 0.16) × (R800 − R670)/(R800 + R670 + 0.16)	[56]
PRI	(R531 − R570)/(R531 + R570)	[57]
PSSR	R800/R635	[58]
REP	700 + (((40 × (R670 + R780)/2) − R700)/(R740 − R700))	[55]
SPVI	0.4 × 3.7 × (R800 − R670) − 1.2 × ((R530 − R670)2)0.5	[59]
SR	R800/R680	[60]
TCARI	3 × ((R700 − R670) − 0.2 × (R700 − R550) × (R700/R670))	[61]
TCARI2	3 × ((R750 − R705) − 0.2 × (R750 − R550) × (R750/R705))	[62]
TCARI_d_OSAVI	TCARI/OSAVI	[61]
TVI	0.5 × (120 × (R750 − R550) − 200 × (R670 − R550))	[63]
VOGI1	R740/R720	[64]
VOGI2	(R734 − R747)/(R715 + R726)	[64]
VOGI3	D715/D705	[64]
DWI	(R816 − R2218)/(R816 + R2218)	[40]
MSGR	(R753 − R708)/(R708 − R681)	[48]
NDII	(R819 − R1600)/(R819 + R1600)	[65]
NMDGI	(R860 − R1640 − R2130)/(R860 + R1640 − R2130)	[66]
NDWGI1	(R820 − R1650)/(R820 + R1650)	[65]
NDWGI2	(R860 − R1240)/(R860 + R1240)	[67]
RDGI	100 × (R1116 − (min(R1120, R1150)))/(R1116)	[68]
RGI1	R1600/R820	[69]
RGI2	R900/R970	[70]
RGI3	R860/R1240	[71]
TDGI	0.02 × (R670 − R550) + 0.01 × (R670 − R480)	[72]
WI	R900/R970	[73]
WLHGI	R676 − (0.5 × (R746 + R665))	[74]

) that were widely used in previous vegetation monitoring studies. Among them, 32 indices were related to the performance of the vegetation, and the other 13 indices were specifically pertinent to the vegetation water status.

Several spectroscopic features (SFs) were retrieved based on continuum removal analysis. Continuum removal transformation highlights the absorption and reflection features of the spectrum by normalizing the reflectance spectra based on a common baseline [21]. For absorption features (AFs), the continuum-removed spectrum was calculated based on the segmented upper hull approach [7] to highlight the AFs that are related to vegetation properties. The same was performed using the bottom baseline to retrieve the Reflectance peak Features (RpFs) from the original spectrum. Moreover, the red-edge (RE) spectral region is a part of the spectrum where an abrupt change in vegetation reflectance from red to NIR exists, respectively, due to the high chlorophyll absorption in red and high vegetation reflectance in NIR [22]. The RE spectral information can provide important insights into vegetation characteristics such as chlorophyll and nitrogen contents, water status, and LAI [23]. As such, several properties of the identified features (AFs and RpFs), as well as the red-edge region (RE), were modelled to characterise each feature (

Table 1. List of feature properties evaluated in the present study.

Feature	Property	Description
Absorption Features (AFs) and Reflectance Peak Features (RpFs)	Area	sum of band depth values within the identified feature
	Width	wavelength difference between upper and lower FWHM
	RMSE_GC_L (Distance to left)	distances to the Gaussian curve from the maximum towards its left
	RMSE_GC_R (Distance to right)	distances to the Gaussian curve from the maximum towards its right
	Max value	observed maximum value within the feature
	Max value wavelength	wavelength of the Max value
Red-edge Region (RE)	RE_min_R	observed minimum reflectance value within the red spectral channel
	RE_min_wl	wavelength of the RE_min_R
	RE_infl_R	reflectance value at the red-edge inflection point
	RE_infl_wl	wavelength of the RE_infl_R
	RE_shld_R	reflectance value at the near-infrared shoulder
	RE_shld_wl	wavelength of the RE_shld_R

). The schematic representation of the spectroscopic features extraction is depicted in Appendix C; Figure A5. The analysis was carried out through the hsdar package [24] in R.

2.5. MLRA Model Generation

In the present study, we utilized five MLRA models that were extensively reported in the growing body of literature on crop traits estimation, including (i) partial least squares regression (PLSR), (ii) random forest (RF), (iii) support vector regression (SVR), (iv) Gaussian process regression (GPR), and (v) neural network (NN).

PLSR [25] reduces the input variables to a small set of uncorrelated components (such as the PCA approach) and then conducts least squares regression on these computed components to predict the response variable. We determined the optimal number of PLSR latent components by minimizing the root mean square error of prediction from cross-validation (RMSEPCV) for each of the crop traits [11].

RF [26,27] uses random samples during model training to develop decision trees. Each tree recursively partitions the data into tree nodes to ensure the predictability of the response variable and then reports the average value as the prediction of all the trees. We optimized the number of trees (Ntree) and the number of variables (Mtry) by minimizing the RMSE values of predictions [27].

SVR [3,28] uses a radial basis kernel function to transfer the data into a higher-dimensional space. We used the Gaussian kernel for SVR models in this research. As such, SVR separates the data into different parts (based on the radial distance between the observations) to minimize the complexity of the nonlinear trends, allowing the assessment of the relationship between training data samples (i.e., support vector) and the response variable. SVR models were optimized through a tune grid-search function identifying the best combinations of cost and gamma parameters [3].

GPR [6,29] with radial basis kernel function provides estimates alongside their uncertainties through a stochastic probability distribution-based Gaussian process, considering possible latent functions and their likelihood. A Gaussian kernel was used in the GPR models. Optimization of GPR models was performed by maximizing the log marginal likelihood according to the sigma values [9].

NN [20,30] is comprised of three main layers: input, hidden, and output. Many artificial neurons with specific rules and connections in the hidden layer exist that learn how to transfer the input data into outputs while minimizing errors and unwanted results. The hidden layer is comprised of two sets of neurons that are useful for considering different prior distributions for two groups of input variables. The Gauss–Newton algorithm was performed to optimize the number of neurons [30]

For each crop trait, the aforementioned MLRA models were generated using 7 different inputs consisting of (a) Full spectra, (b) VIPs of Full spectra, (c) DR PCA-20, (d) DR VIs, (e) VIPs of VIs, (f) DR Spectral Features, and (g) VIPs of SFs. We used the repeated k-fold cross-validation (with k = 10 folds and n = 5 repeats) approach, as proposed by Iqbal et al. [27], for traits with high cardinality (n > 100). Whilst the leave-one-out (LOO) cross-validation approach was performed [11] for low cardinality traits (i.e., CCC, LCC, CWC, LWC, CNC, and LNC).

To assess the quality of the generated MLRA models for retrieving crop traits, we computed the coefficient of determination (R²), the normalized root mean square error (nRMSE) normalized by the range [31] using Equation (1), and the ratio of the performance to deviation (RPD), which is the ratio of the standard deviation of the validation set to the standard error of performance [32] through Equation (2).

$$nRMSE \left(\%\right)=\frac{\sqrt{\frac{1}{n}\sum _{i=1}^n{(Yi-Xi)}^2}}{Xmax-Xmin}\times 100$$

(1)

and

$$RPD=\frac{SD}{\sqrt{\frac{1}{n}\sum _{i=1}^n{(Yi-Xi)}^2}},$$

(2)

where i identifies each sample, n represents the total number of samples, Xi and Yi are the observed and predicted crop trait values of sample i, and SD is the standard deviation of the predicted values. Based on the literature, an RPD > 2.0 indicates a model with excellent prediction power, 1.4 < RPD < 2.0 denotes an intermediate model for quantitative predictions, and RPD < 1.4 shows a nonreliable model, e.g., [33,34].

2.6. Model Transferability to Airborne Image

The robustness of the best-performing MLRA model generated for each crop trait was investigated for all 7 input scenarios while being transferred to the airborne dataset. In addition, the transferability of all the other excellent generated models (i.e., RPD > 2) was investigated. To evaluate the transferability of the generated models, the percent bias (PBIAS) was calculated [31] based on Equation (3). Besides the statistical metrics introduced in the model generation step (Section 2.5), this further indicator allowed us to compare the model’s persistent errors while being applied to airborne datasets.

$$PBIAS \left(\%\right)=\frac{\sum _{i=1}^n(Yi-Xi)}{\sum _{i=1}^n\left(Xi\right)}\times 100,$$

(3)

where i denotes each sample, n is the total number of samples, and Xi and Yi are the observed and predicted crop trait values for sample i.

Therefore, the best-identified combination of the MLRA model and input scenario for each crop trait was considered for mapping and hence, phenotyping analysis purposes. To ensure the models can deal with real-case remotely sensed images of agricultural fields that include either vegetation and soil pixels, 7 soil spectra (selected from the aircraft image) were added to the model generation step, similar to the approach proposed by Verrelst et al. [16]. As such, the aforementioned MLRA models were recalibrated using these 7 soil spectra while considering zero values of LAI, CWC, CCC, and CNC in the model training.

2.7. Phenotyping Analysis

For each crop trait, the best-identified MLRA-input model (recalibrated with soil spectra) was applied to the airborne image to predict a map of the intended crop trait across the experimental field. The average of the pixel values within the experimental plots was used to describe the predicted crop trait for each experimental plot. It is worth mentioning that an inner buffer of 0.5 m from the borders of the plots was considered to avoid the gradient effects at the edge of the plots and potential mixed pixels.

As mentioned in Section 2.1, cultivars, N-treatments, and soil conditions were different agricultural variables in this study. To study the sole impact of each variable, statistical paired t-tests [35] were undertaken to assess the significance of the differences between different groups of each variable. The independence of the considered pairs while comparing different groups of each variable is the key advantage of the paired t-test. More precisely, paired t-test allows comparing different cultivars (i.e., the target genotypes for which performances according to different managements want to be analysed) while considering similar N-treatments and soil conditions in each of the pairs (the same method for evaluating N-treatments as well as soil conditions). As such, paired t-tests were applied to the predicted crop traits from the airborne image. Four levels of significance were considered in the examined differences with asterisks *, **, ***, and ****, denoting the significance levels of p-values less than 0.05, 0.01, 0.001, and 0.0001, respectively.

3. Results

3.1. MLRA Model Generation for Crop Trait Estimation

Figure 4 represents the performance assessment of the five generated MLRA models with respect to the seven different input scenarios for all the crop traits. A comprehensive overview of all the statistical metrics is presented in Appendix D,

Table A2. Statistics comparing the performance of different MLRA models for different input scenarios for all the studied crop traits. The coefficient of determination (R2), the normalized root mean square error (nRMSE), and the performance-to-deviation ratio (RPD) are presented. The best-identified models for each of the crop traits are colour-coded with dark green, moderate green, and red indicating excellent (RPD > 2.0), intermediate (1.4 < RPD < 2.0), and 312 nonreliable (RPD < 1.4) models, respectively.

		LAI (#232)			CWC (#61)			CCC (#64)			CNC (#44)			LWC (#61)			LCC (#64)			LNC (#44)
		R2	nRMSE	RPD	R2	nRMSE	RPD	R2	nRMSE	RPD	R2	nRMSE	RPD	R2	nRMSE	RPD	R2	nRMSE	RPD	R2	nRMSE	RPD
Full spectra	PLSR	0.72	11.97	1.89	0.77	10.89	2.12	0.7	14.61	1.83	0.74	14.38	1.97	0.58	17.55	1.91	0.34	16.69	1.41	0.2	17.61	1.59
	RF	0.59	14.59	1.55	0.42	17.42	1.32	0.5	18.86	1.42	0.42	21.41	1.32	0.38	22.23	1.27	0.14	19.31	1.06	0.06	17.61	1.06
	SVR	0.69	12.75	1.77	0.46	16.69	1.38	0.52	18.59	1.43	0.37	22.46	1.25	0.36	22.23	1.25	0.23	18.11	1.14	0.01	17.61	1.06
	GPR	0.68	13.46	1.68	0.45	17.42	1.32	0.53	18.78	1.42	0.49	20.38	1.38	0.37	23.23	1.24	0.21	18.2	1.13	0.01	17.61	1.06
	NN	0.82	9.66	2.34	0.63	14.52	1.58	0.66	16	1.67	0.69	15.86	1.78	0.39	22.23	1.25	0.34	17.37	1.18	0.09	17.61	1.07
Full Spectra—VIPs	PLSR	0.72	11.9	1.9	0.76	10.89	2.27	0.72	13.9	2.01	0.74	14.42	1.96	0.62	23.4	2.45	0.35	16.48	1.43	0.21	17.6	1.6
	RF	0.61	14.24	1.58	0.47	16.69	1.38	0.58	17.3	1.54	0.44	20.96	1.35	0.4	23.4	1.22	0.22	18.09	1.14	0.09	17.61	1.06
	SVR	0.69	12.7	1.78	0.45	16.69	1.38	0.62	16.61	1.6	0.47	20.55	1.37	0.37	22.7	1.26	0.34	17.03	1.21	0.07	17.61	1.06
	GPR	0.68	13.16	1.71	0.51	15.97	1.44	0.61	17.12	1.56	0.51	19.78	1.43	0.43	22	1.3	0.29	17.32	1.19	0.02	17.61	1.06
	NN	0.79	10.33	2.18	0.72	14.51	1.59	0.68	15.12	1.76	0.8	12.57	2.25	0.62	17.55	1.63	0.38	16.28	1.26	0.15	17.61	1.06
Dimensionally-reduced spectra—PCA20	PLSR	0.81	9.9	2.28	0.8	7.26	3.18	0.76	12.92	2.29	0.74	14.38	1.96	0.58	18.25	1.91	0.33	16.68	1.41	0.2	17.61	1.59
	RF	0.76	11.39	1.98	0.7	14.51	1.59	0.73	14.23	1.87	0.62	17.36	1.62	0.4	21.06	1.29	0.28	17.4	1.18	0.06	17.62	1.06
	SVR	0.76	11.18	2.02	0.63	14.51	0.59	0.65	16.5	1.61	0.47	21.91	1.29	0.39	22.23	1.28	0.27	17.46	1.18	0.22	17.62	1.06
	GPR	0.75	12.07	1.87	0.7	14.51	1.59	0.67	17.22	1.55	0.59	22.35	1.26	0.46	21.06	1.31	0.27	17.67	1.16	0.06	17.62	1.06
	NN	0.82	9.48	2.38	0.79	7.26	3	0.72	14.12	1.89	0.67	16.28	1.73	0.54	18.72	1.48	0.27	17.74	1.16	0.09	17.62	1.06
Dimensionally-reduced spectra—VIs	PLSR	0.74	11.5	1.96	0.77	7.26	3.18	0.76	11.77	2.26	0.69	15.63	1.81	0.51	18.72	1.51	0.39	14.14	1.46	0.08	17.62	1.06
	RF	0.75	11.39	1.98	0.76	14.51	1.59	0.68	14.99	1.78	0.61	17.56	1.61	0.52	23.4	1.22	0.22	18.45	1.12	0.2	17.61	1.06
	SVR	0.74	11.51	1.96	0.53	14.51	1.59	0.58	17.56	1.52	0.51	19.66	1.43	0.53	18.72	1.46	0.25	17.89	1.15	0.13	17.61	1.06
	GPR	0.72	12.23	1.85	0.62	14.51	1.59	0.57	17.56	1.52	0.48	20.71	1.36	0.45	21.06	1.34	0.22	18.08	1.14	0.16	17.61	1.06
	NN	0.8	10.03	2.25	0.77	14.51	1.59	0.72	13.99	1.9	0.67	16.17	1.74	0.64	16.38	1.66	0.31	17.19	1.2	0.41	15.47	1.2
Dimensionally-reduced spectra—VIs—VIPs	PLSR	0.74	11.45	1.97	0.76	7.26	3.18	0.75	13.37	2.15	0.73	14.68	1.92	0.62	16.36	1.75	0.29	17.25	1.34	0.35	17.61	1.59
	RF	0.75	11.21	2.01	0.8	10.16	2.27	0.73	13.91	1.92	0.65	16.5	1.71	0.55	23.4	1.22	0.25	18	1.14	0.23	17.61	1.06
	SVR	0.77	10.79	2.09	0.67	13.06	1.76	0.68	15.03	1.77	0.57	18.51	1.52	0.57	23.4	1.22	0.26	17.19	1.61	0.24	17.61	1.06
	GPR	0.77	11.18	2.02	0.66	13.79	1.67	0.7	14.91	1.79	0.59	18.27	1.54	0.56	23.4	1.22	0.22	18.19	1.13	0.24	17.61	1.06
	NN	0.8	10.18	2.22	0.78	20.89	2.12	0.73	13.79	1.93	0.71	15.02	1.88	0.64	11.7	2.45	0.33	16.83	1.22	0.28	17.61	1.06
Dimensionally-reduced spectra—Features	PLSR	0.77	10.75	2.1	0.62	7.26	3.18	0.53	13.48	1.98	0.63	17.04	1.66	0.49	15.21	1.82	0.22	14.2	1.45	0.09	17.61	0.97
	RF	0.77	10.76	2.1	0.72	14.51	1.59	0.69	14.69	1.81	0.59	17.84	1.58	0.58	23.4	1.22	0.29	17.25	1.19	0.21	17.61	1.06
	SVR	0.74	11.53	1.96	0.59	14.51	1.56	0.58	17.27	1.54	0.55	18.81	1.5	0.52	18.72	1.44	0.25	17.8	1.6	0.07	17.81	1.04
	GPR	0.73	12	1.88	0.66	7.26	1.67	0.61	17.23	1.55	0.58	18.67	1.51	0.59	18.72	1.5	0.21	18.23	1.13	0.08	27.37	0.68
	NN	0.78	10.76	2.1	0.65	7.26	1.7	0.66	15.57	1.71	0.65	16.79	1.68	0.58	18.72	1.53	0.19	18.97	1.08	0.12	17.61	1.06
Dimensionally-reduced spectra—Features—VIPs	PLSR	0.77	10.75	2.1	0.71	14.51	1.59	0.72	13.98	2.14	0.68	15.91	1.77	0.67	16.38	1.88	0.3	17.16	1.34	0.17	17.61	1.06
	RF	0.79	10.39	2.17	0.77	10.89	2.12	0.73	13.79	1.93	0.63	17.1	1.65	0.66	11.7	2.45	0.37	16.27	1.26	0.32	17.61	1.06
	SVR	0.77	11.01	2.05	0.72	12.34	1.87	0.68	14.96	1.78	0.61	17.46	1.62	0.6	23.4	1.22	0.42	15.58	1.32	0.18	17.61	1.06
	GPR	0.76	11.2	2.02	0.72	12.34	1.87	0.66	16.95	1.57	0.6	17.93	1.57	0.66	11.7	2.45	0.39	16.13	1.28	0.18	17.61	1.06
	NN	0.82	9.64	2.34	0.78	10.89	2.12	0.73	13.71	1.94	0.64	16.78	1.68	0.69	11.7	2.45	0.36	16.59	1.24	0.11	17.61	1.06

. Overall, the accuracy of the models was remarkably higher for canopy-level traits (i.e., LAI, CWC, CCC, and CNC) compared to leaf-level traits (i.e., LWC, LCC, and LNC). What was noteworthy was the impact of DR on the accuracy of crop trait retrieval. More specifically, DR MLRA models were found to attain a higher accuracy compared to the full spectra models for retrieving wheat traits. PCA-20 and Features-VIPs seemed to be the best-examined DR techniques, with DR PCA-20 indicating the best retrieval models for LAI, CWC, and CCC, and DR Features-VIs revealing the best LWC and LCC retrieval models (Figure 5).

In general, PLSR and NN models worked better for the examined traits compared to RF, SVR, or GPR models. PLSR and NN yielded 44 and 34 good (considering also excellent) models (i.e., RPD > 1.4), followed by 26, 26, and 25 in SVR, GPR, and RF, respectively.

All LAI retrieval models indicated either good or excellent prediction capability, with RPD values greater than 1.4 (Figure 4). It is worth noting that NN on all the input scenarios yielded excellent models (RPD > 2.0) for LAI retrievals. The best-performing LAI model was attained using DR PCA-20 in NN with an RPD of 2.38 (R² = 0.82, nRMSE = 9.48%; Figure 5a).

The majority of CWC retrieval models were categorized as good, though a few unreliable instances were also found (Figure 4). PLSR showed excellent models for CWC retrieval for all the input scenarios except on DR Features-VIPs when a good model was achieved (Figure 4). CWC was most accurately retrieved using DR PCA-20 in PLSR with an RPD of 3.18 (R² = 0.80, nRMSE = 7.26%; Figure 5b).

Similar to LAI, all the CWC retrieval models showed good or excellent prediction capability (Figure 4). PLSR generated excellent CCC models for most of the input scenarios except for using Full spectra and DR Features where the models categorised as good (Figure 4). As shown in Figure 5c, DR PCA-20 in PLSR indicated the best CCC retrieval model (RPD = 2.29, R² = 0.76, nRMSE = 12.92%).

Most of the CNC retrieval models were classified as good, with a few unreliable models (Figure 4). NN and PLSR indicated relatively better models for CNC retrieval (Figure 4). It is worth stating that the best-performing CNC model was achieved using Full spectra-VIPs in NN as an excellent model with an RPD of 2.25 (R² = 0.80, nRMSE = 12.57%; Figure 5d).

3.2. MLRAs Model Transferability

The transferability results of the best-performing MLRA models (RPD > 2 during the model generation phase), applied to the airborne dataset, are presented in

Table A3. Transferability results of the best-performing MLRA models (RPD > 2 during the model generation phase) for the canopy-level traits. The coefficient of determination (R2), the normalized root mean square error (nRMSE), the performance-to-deviation ratio (RPD), and the percent bias (PBIAS) are presented. The best-identified models for each of the crop traits are colour-coded with dark green, moderate green, and red indicating excellent (RPD > 2.0), intermediate (1.4 < RPD < 2.0), and 312 nonreliable (RPD < 1.4) models, respectively.

			LAI (#89)				CWC (#22)				CCC (#24)				CNC (#22)
		R2	nRMSE	RPD	PBIAS	R2	nRMSE	RPD	PBIAS	R2	nRMSE	RPD	PBIAS	R2	nRMSE	RPD	PBIAS
Full spectra	PLSR	-	-	-	-	0.19	45.75	0.58	−30.50	0.43	43.88	0.71	−26.00	0.51	21.27	1.23	10.00
	RF	-	-	-	-	-	-	-	-	-	-	-	-	-	-	-	-
	SVR	-	-	-	-	-	-	-	-	-	-	-	-	-	-	-	-
	GPR	-	-	-	-	-	-	-	-	-	-	-	-	-	-	-	-
	NN	0.62	19.46	1.15	−12.50	-	-	-	-	-	-	-	-	-	-	-	-
Full Spectra—VIPs	PLSR	-	-	-	-	0.63	30.50	0.87	−16.10	0.00	85.02	0.36	19.20	0.64	18.34	1.43	9.00
	RF	-	-	-	-	-	-	-	-	-	-	-	-	-	-	-	-
	SVR	-	-	-	-	-	-	-	-	-	-	-	-	-	-	-	-
	GPR	-	-	-	-	-	-	-	-	-	-	-	-	-	-	-	-
	NN	0.69	14.95	1.50	7.30	-	-	-	-	-	-	-	-	0.49	52.85	0.50	−42.10
Dimensionally-reduced spectra—PCA20	PLSR	0.52	20.18	1.11	−10.00	0.21	45.70	0.60	−30.10	0.45	37.26	0.83	−19.70	0.53	21.23	1.25	−9.30
	RF	-	-	-	-	-	-	-	-	-	-	-	-	-	-	-	-
	SVR	0.52	22.34	1.00	−13.50	-	-	-	-	-	-	-	-	-	-	-	-
	GPR	-	-	-	-	-	-	-	-	-	-	-	-	-	-	-	-
	NN	0.55	18.56	1.21	9.80	0.26	30.50	0.87	9.80	-	-	-	-	-	-	-	-
Dimensionally-reduced spectra—VIs	PLSR	-	-	-	-	0.39	15.25	1.74	1.20	0.00	85.02	0.36	19.20	0.29	49.07	0.54	40.00
	RF	-	-	-	-	-	-	-	-	-	-	-	-	-	-	-	-
	SVR	-	-	-	-	-	-	-	-	-	-	-	-	-	-	-	-
	GPR	-	-	-	-	-	-	-	-	-	-	-	-	-	-	-	-
	NN	0.00	25.59	0.88	0.40	-	-	-	-	-	-	-	-	-	-	-	-
Dimensionally-reduced spectra—VIs—VIPs	PLSR	-	-	-	-	0.55	30.50	0.87	16.20	0.02	44.61	0.69	21.30	0.04	35.45	0.74	−20.00
	RF	0.32	39.64	0.57	−31.90	0.22	61.00	0.43	−40.60	-	-	-	-	-	-	-	-
	SVR	0.03	37.12	0.60	−27.40	-	-	-	-	-	-	-	-	-	-	-	-
	GPR	0.12	27.57	0.81	−15.90	-	-	-	-	-	-	-	-	-	-	-	-
	NN	0.00	26.13	0.86	−7.80	0.00	76.25	0.35	−47.60	-	-	-	-	-	-	-	-
Dimensionally-reduced spectra—Features	PLSR	0.55	34.41	0.65	−27.80	0.23	61.00	0.43	−35.10	0.00	85.02	0.36	19.20	0.25	53.13	0.9	−44.00
	RF	0.61	33.33	0.67	−27.60	-	-	-	-	-	-	-	-	-	-	-	-
	SVR	-	-	-	-	-	-	-	-	-	-	-	-	-	-	-	-
	GPR	-	-	-	-	-	-	-	-	-	-	-	-	-	-	-	-
	NN	0.31	47.21	0.48	−35.80	-	-	-	-	-	-	-	-	-	-	-	-
Dimensionally-reduced spectra—Features—VIPs	PLSR	0.18	72.61	0.31	−63.00	0.67	30.50	0.87	5.10	0.13	85.22	0.36	−57.80	0.55	23.82	1.10	−10.00
	RF	0.66	27.21	0.83	−21.40	0.46	61.00	0.43	−36.60	-	-	-	-	-	-	-	-
	SVR	0.31	35.5	0.63	−26.10	-	-	-	-	-	-	-	-	-	-	-	-
	GPR	0.31	27.93	0.80	−16.70	-	-	-	-	-	-	-	-	-	-	-	-
	NN	0.07	63.06	0.36	−47.10	0.34	45.75	0.58	−26.30	-	-	-	-	-	-	-	-

of Appendix E. It is worth stating that validation metrics were obtained with independent data with respect to the ones used for model training. Spectra were derived using airplane images, and ground reference data were collected in the field. As described in Section 2.6, the best combination of the MLRA model and input scenario was selected for mapping the studied crop traits. Figure 6 illustrates the transferability assessment of the best MLRA-input cases for the examined canopy-level traits, with and without adding soil spectra during the model generation phase. LAI revealed the highest transferability accuracy with an RPD value of 1.5 (R2 = 0.69, nRMSE = 14.95%, PBIAS = 7.30%) when utilizing the Full spectra-VIPs in NN (Figure 6a). CNC also displayed good transferability accuracy with an RPD of 1.43 (R2 = 0.64, nRMSE = 18.34%, PBIAS = 9%) while using the Full spectra-VIPs in PLSR (Figure 6g). Transferability accuracy was, however, less significant for CWC (Figure 6c) and CCC (Figure 6e) with RPD values below 1.4. Moreover, adding the soil spectra in the model training phase slightly diminished the model accuracy (RPD values) when being transferred to the aircraft dataset (Figure 6).

As such, the generated crop trait maps (also using the soil spectra) are depicted in Figure 7. Overall, LAI, CWC, CCC, and CNC maps revealed rather consistent patterns that also well represented the real-field situation. The LAI map indicated better wheat development in the middle of the field and, more specifically, within the bottom-middle part (Figure 7). This was also corroborated by higher values of CWC, CCC, and CNC in the bottom–middle zone of the field. In addition, the left side of the field represented a relatively worse crop status for all the studied crop traits.

3.3. Wheat Phenotyping

No statistical difference was found between the LAI values of the different wheat cultivars (Figure 8), except between V2 and V4, in which V2 exhibited better crop development. However, significant differences were found in LAI values in accordance with the various N-treatment rates. The highest LAI values were observed for N4, followed by N1 and N2 treatments.

Furthermore, no impact of soil tillage type (minimum vs. traditional) was found on wheat LAI values. However, it is worth mentioning that in the L1 (minimum tillage) south–western part, an excess of weeds was found.

CWC was not statistically different amongst V1, V2, and V3 cultivars, whilst V4 indicated significantly fewer CWC values when compared to the other cultivars (Figure 8). CWC, however, was not different between N1, N2, and N4 treatments, while only N3 exhibited significantly fewer CWC values than the other classes. Similar to LAI, CWC did not undergo any particular impact on soil preparation.

CCC was only significantly different between V1 and V2 cultivars (Figure 8). N3 denoted significantly fewer CCC values compared to the other N-treatments in agreement with other traits. Moreover, L2 soil preparation (traditional deep ploughing) was found to significantly increase CCC compared to L1 (minimum tillage).

Eventually, V2 showed significantly higher CNC values than V1 and V4 cultivars (Figure 8). In accordance with the observed trends in LAI, CWC, and CCC, N3 denoted significantly fewer CNC values compared to the other N-treatments. In addition, soil preparation (L1 and L2) did not have any particular effect on the observed CNC values.

4. Discussion

4.1. Canopy vs. Leaf Traits Accuracy Retrieval

Experimental results demonstrate the high accuracy that can be obtained from hyperspectral data and appropriate MLRA coupled with DR technique for the estimation of canopy level crop traits. The models for retrieving leaf-level traits were less accurate (Figure 4). LWC, LCC, and LNC were most accurately retrieved using DR Features-VIPs in NN, DR Features-VIPs in SVR, and DR VIs in NN, respectively, with the corresponding RPD (resp. R² and nRMSE) values of 2.45 (0.69, 11.70%), 1.32 (0.42, 15.58%), and 1.20 (0.41, 15.47%), as shown in Figure 5e–g. This corroborates the previous finding of Candiani et al. [9] where retrieving the leaf-level traits from canopy spectra from remote was more critical than canopy-level trait estimation.

4.2. Tranferability from Ground to Aerial Data

Results demonstrate that MLRA derived from ground data are robust and can be transferred to other hyperspectral imaging sensors for mapping application. Independent validation reveal reduced performance, with respect to cross validation results, in our experimental case. Despite an expected certain level of uncertainty when transferring models from one calibration data set to other independent data, results show satisfying performance for mapping demonstration useful for phenotyping activity.

It is worth mentioning that a reduction in transferability performance observed could also be related to the sub-optimal spectral resampling of the PRISMA-like CASI-SASI dataset in the SWIR region useful for nitrogen and, in particular, water detection. In our opinion, DR techniques, such as PCA that exploits spectral information from all available bands, can minimize such problems. In addition, the modelling accuracy decreased while adding the soil spectra in the training (Figure 6). The latter is in line with the recent finding of Verrelst et al. [16], such decision however increases model robustness in case of sparse canopy with soil patches. This is the case of the eastern part of the field that indicated the least modelled crop trait values (Figure 7), in which wheat plants were suffering at emergence from high water accumulations that also caused crop failure, and hence, a sparser vegetation cover (last four southern plots).

4.3. Phenotyping Considerations

The performed phenotyping analysis represents a proof of concept of the potentiality of the hyperspectral imaging system coupled with MLRA and DR techniques. The experiment allowed to investigate the effect of different agro-practices on different genotypes development and growth. Concerning the fertilisation strategy, it is worth stating that N3 caused remarkably lower LAI and CCC values compared to the other N-treatments (Figure 8); this was the strategy with a 30% excess of fertilization. Results demonstrate that conservative approaches with high fertilisation can produce negative effect on crop performance beside increasing environmental and economic impact. On the other hand, the N4 strategy was based on a variable rate approach for which the second fertilization dose, overall less than the standard one, was modulated using expert knowledge according to LAI and LCC values measured in the field (see Figure 2d). This approach overall produces comparable or even better results on crop traits with respect to the standard treatments showing the importance of site specific management based on remote sensing data for precision farming application.

Regarding genetic differences, significantly less CWC values in V4 (Figure 8) can be related to a faster phenological development of V4 that determines not-insignificant number of plants already in grain-filling conditions.

Although we demonstrated a strong capability of MLRAs and DR techniques to estimate wheat crop traits, future work should still test the proposed methodology over other study areas and datasets to draw more generalized findings. In this research, we chose only to focus on estimating crop traits based on MLRAs and DR techniques that are more appropriate for a limited training dataset. Future research may also consider using deep learning (DL) methods for similar purposes. To test DL solutions, which are usually data-hungry [16,36], appropriate methods able to exploit limited data should also be investigated.

This paper provided a proof of concept of utilizing MLRA and DR techniques based on hyperspectral data to help crop traits monitoring throughout the cropping seasons for precision farming applications and the implementation of high-throughput phenotyping solutions.

5. Conclusions

In this paper, we focused on modelling different wheat crop traits based on MLRA and DR techniques using hyperspectral data. Canopy-level traits were accurately retrieved using MLRA models, particularly through PLSR and NN. Although we denoted that MLRA models were less accurate while retrieving leaf-level traits, more work is yet needed to draw generalized conclusions. It is also important to highlight the promising role of DR techniques when conducting phenotyping studies using hyperspectral data. Our results indicated that MLRA models based on DR data can better retrieve crop traits than when using the full hyperspectral data. PCA-20 and Features-VIPs were found to be the best DR techniques amongst the tested methods.

In the present study, we examined the performance of the generated models when being transferred to an airborne dataset. As such, the studied crop traits were modelled across the entire wheat field (96 plots vs. 24 ground measurements) based on the airborne image as an overarching objective for field phenotyping. The predicted crop trait maps matched well with the real-field characteristics, and amongst the studied traits, LAI had the highest transferability accuracy. This approach allowed us to perform accurate wheat phenotyping analysis with respect to different agricultural variables (i.e., cultivars, N-treatments, and soil conditions). N-treatment was found to have a significant impact on wheat crop traits, whilst less significant differences were reported between the cultivars. It is worth mentioning that the variable rate strategy (N4) that used ground LAI and LCC measurements to support expert-based decisions in defining the appropriate dose resulted in efficiency for crop growth and, at the same time, reduced nitrogen input. In this framework, crop parameter maps obtained through the MLRA strategy can represent a fundamental input to creating prescription maps to implement precision farming solutions. Moreover, no particular effect of soil ploughing condition on wheat crop traits was observed. In principle, the tested methodology can be applied to other crops and agricultural treatments. Future research should be aimed at repeating this approach over other hyperspectral datasets and across other agricultural fields to draw stronger conclusions.