Retrieval of Maize Leaf Area Index Using Hyperspectral and Multispectral Data

Field spectra acquired from a handheld spectroradiometer and Sentinel-2 images spectra were used to investigate the applicability of hyperspectral and multispectral data in retrieving the maize leaf area index in low-input crop systems, with high spatial and intra-annual variability, and low yield, in southern Mozambique, during three years. Seventeen vegetation indices, comprising two and three band indices, and nine machine learning regression algorithms (MLRA) were tested for the statistical approach while five cost functions were tested in the look-up-table (LUT) inversion approach. The three band vegetation indices were selected, specifically the modified difference index (mDId: 725; 715; 565) for the hyperspectral dataset and the modified simple ratio (mSRc: 740; 705; 865) for the multispectral dataset of field spectra and the three band spectral index (TBSIb: 665; 865; 783) for the Sentinel-2 dataset. The relevant vector machine was the selected MLRA for the two datasets of field spectra (multispectral and hyperspectral) while the support vector machine was selected for the Sentinel-2 data. When using the LUT inversion technique, the minimum contrast estimation and the Bhattacharyya divergence cost functions were the best performing. The vegetation indices outperformed the other two approaches, with the TBSIb as the most accurate index (RMSE = 0.35). At the field scale, spectral data from Sentinel-2 can accurately retrieve the maize leaf area index in the study area.

The combination of POSPECT and SAIL (PROSAIL) interconnects the variation of the spectral reflectance of the canopy, which is mainly related to the biochemical attributes of the leaf, with its directional variation, which is mainly related to the canopy architecture and the contrast between soil/vegetation. This linkage is essential for the simultaneous estimation of the biophysical and structural variables of the vegetation [21]. However, the parameterization of PROSAIL seems quite challenging in such a way that most of the studies use nominal values or ranges applied in other studies carried out with the same crops, but in different places. In Table A1 (Appendix A), we present the PROSPECT and SAIL models, including parameters used in several crops and vegetation types.
Overall, although several mentioned studies have achieved good results estimating and mapping crop biophysical properties through different remote sensing data and methods, according to [22], there are still gaps to be explored in terms of the use and operationalization of these data and methods. The diversity of bands selected for estimating the same variable, for the same crop, in different studies, and the variability of the intervals of parameters used in the calibration of the RTMs for the same crop, emphasize the impact of the radiometric characteristics of the sensors, the agricultural practices (which determine the canopy structure), the environmental conditions, and sampling conditions on the models' results. This assumption is even more important when remote sensing data and methods are to be applied in low-input crop systems and highly heterogeneous farming systems, as is the case of the present study area. In this type of farming system, the limitation of the production factors, such as fertilizers and irrigation, leads to crops with very low LAI and high infield heterogeneity, which constitutes an additional challenge in the application of remote sensing techniques.
In this context, the main goal of this work is testing the applicability of different types of remote sensing data (hyperspectral and multispectral) and methods (statistical and physical) for developing models to predict LAI in maize crop systems characterized by low input factors, high spatial and intra-annual variability, and low yield in the south of Mozambique. The specific objectives include (i) to compare the performance of multi and hyperspectral data in estimating the LAI of maize in the study region; and (ii) to test the performance of the RTMs and the statistically-based models in the estimation of LAI of maize.

Study Area Description
The current study was conducted in maize crop fields located in the district of Vilankulo, within the province of Inhambane in southern Mozambique (lat.: 21 • 58 S, long.: 035 • 09 E, 31.83 m a.s.l) ( Figure 1).
The district of Vilankulo is characterized by a semi-arid to arid climate, with sandy soils of low fertility. Figure 2 presents the average precipitation and temperature in the study area for a reference period and for the three years of data collection. For the period studied (2015-2018), the average of total annual rainfall was 403 mm, which is 50% below the 804.4 mm registered in the period . The average annual temperature of 26.9 • C during the study period was warmer than that recorded in the period of 1991-2015 (24.5 • C). The summer and rainy season (January, February, March, October, November, and December) accounted for the highest amount of rainfall for both the reference period (79.8%) and study period (81.2%). The driest conditions were recorded in 2015 with 335.3 mm of total annual rainfall and an average temperature of 29.0 • C.     [23].

Characteristics of the Sampled Fields
The data collection took place in four maize fields in 2015 (field 1), 2016 (field 2 and field 3) and 2018 (field4). In 2015 and 2016, the crop cycle coincided with the cold and dry season (April-September) while in 2018, the crop was grown during hot and rainy season. In each field, two 400 m 2    [23].

Characteristics of the Sampled Fields
The

Characteristics of the Sampled Fields
The data collection took place in four maize fields in 2015 (field 1), 2016 (field 2 and field 3) and 2018 (field 4). In 2015 and 2016, the crop cycle coincided with the cold and dry season (April-September) while in 2018, the crop was grown during hot and rainy season. In each field, two 400 m 2 parcels were established for monitoring and data collection. Table 1 summarizes the cropping practices in each sampled field. The field LAI data was obtained from an allometric model, which estimates the plant total leaf area based on biometric variables (length and width) and the number of leaves. The model variables were collected during the years, 2015, 2016, and 2018, in seven plants randomly selected and marked within each parcel. The maize measurements (allometrics descriptors) were recorded in different crop phenologic stages ( Table 2) described in accordance with the leaf collar method [24]. This allometric model was calibrated and validated during the 2015 field campaign using the maize variety, 'PAN 53', and explained 90% of the maize total leaf area variability at different stages of crop development (R 2 = 0.90; n = 30; p < 0.000) [1]. Additionally, an independent validation of the model was performed using data from the crop variety, 'PAN 67', grown in 2018 (R 2 = 0.914, n = 16, p < 0.000; data not published).

Field Equivalent Water Thickness and Dry Mater Content
For the estimation of the equivalent water thickness (Cw,) and dry mater content (Cm, g/cm 2 ), the first full developed leaf, defined from the top to the bottom of each sampled plant, was cut and enclosed in a sealed plastic bag and brought to the laboratory inside a cooler. In the laboratory, a 2 cm diameter disc was cut in each leaf and the wet weight taken before the discs were dried at 65 • C in an oven until constant weight (dry weight). The Cw (cm) was then calculated according to the following equation: where F w (g) and D w (g) are, respectively, the fresh and dry weights of the leaf discs; A (cm 2 ) is the area of the disc; and dw (1 g cm −3 ) is the water density. For the estimation of Cm, we used the following equation:

Field Spectral Data Collection
The field spectral data were acquired using an Apogee hyperspectral spectroradiometer with sensitivity in the range of 236 nm-1100 nm and a spectral resolution of 1 nm. The instrument is attached to an optic cable with a field of view (FOV) of 30 • . The measurements were taken at the nadir position, keeping a constant distance of about 5 cm to the spot, and were carried out around solar noon (between 11 a.m. and 1 p.m.), when the changes in solar zenith are minimal. Standard procedures of the equipment calibration were followed, namely the dark object correction and the acquisition of white reference reflectance of the Spectralon. The reference reflectance was updated every 10 min during the measurements. The spectra readings were done in the same dates (Table 2) and points where the data for LAI determination were collected. In the phenological stages, V3, V6, Remote Sens. 2018, 10, 1942 7 of 30 and V8, while the fraction of canopy cover was still very low, two additional readings were done, one over the interface between two plants and another over bare soil. In each point, 10 replications were made, with each an average of 20 automatic readings.

Spectral Data from Sentinel-2
The multispectral sensors imaging (MSI) on board of the Sentinel-2A and Sentinel-2B, launched, respectively, in June 2015 and March 2017, have a temporal resolution of 10 days and present 4 bands (2, 3, 4, and 8) with 10 m of spatial resolution and 6 bands (5, 6, 7, 8a, 11, and 12) with 20 m of spatial resolution. The high spatial and temporal resolutions and the spectral intervals covered by these bands make Sentinel-2 potentially more interesting than other multispectral sensors for studies related with LAI [25][26][27].
Atmospherically corrected Sentinel-2A and Sentinel-2B images of the study area, including their value added product bands, were acquired from the data service platform developed by [28]. For all the pixels covering each sampling parcel, the spectral data of all the original Sentinel-2 bands (SP_S2) and of the Sentinel LAI products (LAI_S2) were extracted and then averaged, using a 10 m spatial filter applying the zonal statistics tool of QGIS software, version 2.14.14. The Sentinel-2 images were selected considering the closeness to field spectral data collection dates and low cloud cover. The images from the following dates were then selected and downloaded: 30/06/2016; 20/07/2016; 30/07/2016/; 28/09/2016; 21/01/2018; 05/02/2018, and 02/03/2018. The dates of satellite acquisition and field data collection did not always match due to frequent high cloud cover conditions in the study area at the time of satellite overpass. However, nearby satellite image dates, corresponding to the same crop phenological stage, were always considered.

Field Spectral Data Pre-Processing
The reflectance of each sampling point was achieved by averaging the 10 replications. For the stages, V3, V6, and V8, the average included the reflectance of the plant, bare soil, and interface between two plants. A weighted average of the three ground cover components was considered to obtain a reflectance value representative of maize growth conditions at each sampling point and allow the comparability with satellite data. Due to noise in the initial portion of the spectra, data analysis was performed considering data in the interval of 400 nm-1100 nm of the spectra. Afterwards, these data were spectrally aggregated to simulate hyperspectral and multispectral datasets. For the hyperspectral dataset, a 10 nm spectral resolution aggregation (FSP_10) was considered while for the multispectral dataset, the Sentinel-2 band setting was used for aggregation (FSP_S2). Because the field spectra does not match the whole Sentinel-2 spectra, the aggregation was limited to 8 bands (2 A-8 A) acknowledged as being more convenient in LAI estimations [2,29]. Table 3 presents the spectral vegetation indices (VI) assessed as predictors for LAI statistical modelling. The VIs included band combinations in the form of the simple ratio (SR), modified simple ratio (mSR), normalized difference (ND), modified normalized difference (mND), difference between bands (DI), and modified difference between bands (mDI). For each VI, all possible band combinations were tested considering the three types of data: Hyperspectral dataset (FSP_10), the multispectral dataset (FSP_S2), and the Sentinel-2 data (SP_S2). The FSP_10 and FSP_S2 were correlated only to the field leaf area index (Field_LAI) while the SP_S2 was correlated to both the average Field_LAI and Sentinel-2 leaf area index product (LAI_S2). The linear, exponential, and second degree polynomial functions were evaluated for modelling the LAI based on VIs. The analysis of the best band combination per VI and the modelling approach were carried out in the Spectral Index Toolbox of the software Automated Radiative Transfer Models Operator (ARTMO) [15,30].  [45] ρ-Reflectance at given wavelength; SR-simple ratio; ND-normalized difference; mSR-modified simple ratio; mND-modified normalized difference; TBSI-three band spectral index; DI-difference index; EVI-enhanced vegetation index; SAVI-soil adjusted vegetation index; NDVI-normalized difference vegetation index; VARI-visible atmospherically resistant index; CARI-chlorophyll absorption ratio index; mCARI-modified chlorophyll absorption ratio index; CI Green -green chlorophyll index; ARI-anthocyanin reflectance index; mARI-modified anthocyanin reflectance index; SIPI-structure insensitive pigment index; PSRI-plant senescence reflectance index; RVSI-reflectance vegetation stress index. Figure 3 presents the flowchart of the two statistical modelling approaches.

Statistical Modelling of LAI Based on Machine Learning Regression Algorithms
The machine learning regression algorithms (MLRA) are nonparametric models that are adjusted to predict a variable of interest using a training dataset of input-output data pairs, derived from synchronised measurements of the parameter and the corresponding reflectance observations. Numerous nonparametric regression algorithms are available in the statistics and Remote Sens. 2018, 10,1942 9 of 30 machine learning literature, and they are being increasingly applied for vegetation biophysical parameters' retrieval [14,16,17,46,47].
Similarly to the VIs, the MLRA algorithms were evaluated using the three types of data: Hyperspectral dataset (FSP_10), multispectral dataset (FSP_S2), and Sentinel-2 data (SP_S2). The FSP_10 and FSP_S2 were correlated only to the field leaf area index (Field_LAI) and the SP_S2 was correlated to both the Field_LAI and Sentinel-2 LAI product (LAI_S2). Table 4 presents the MLRA evaluated in the present study and the modelling process is presented in Figure 3. Prior to the application of the MLRA algorithms, an analysis of the best number and combination of bands for estimating the LAI using the different datasets was performed based on the band analysis toolbox. The algorithms were also implemented on the MLRA toolbox [17] of the software Automated Radiative Transfer Models Operator (ARTMO). Table 4. Machine learning regression algorithms tested for LAI modelling.

Retrieval of LAI Based on Radiative Transfer Models
To perform an LAI retrieval through RTMs, one need to first build up a look-up-table (LUT), which is a database of simulated canopy reflectance and their corresponding set of input parameters. Afterwards, an inversion technique is applied to search and identify the parameters' combination that yields the best fit between the measured and LUT reflectance [18,56].
Simulation of the Look-Up- Table (LUT) The PROSAIL model, a widely used RTM model for biophysical parameters' retrieval [14,57], was used to simulate reflectance data that were stored in an LUT. Table 5 presents the input parameters' settings used for the LUT simulations. The input parameters relative to LAI, water equivalent thickness (Cw), dry matter content (Cm), hotspot (hspot), and the reflectance of dry and wet soil were set according to field data. The parameter of leaf structure (N) was set to the 1-1.4 interval as proposed by [2,58,59] for various crops, including maize. The interval of values considered for the chlorophyll a and b content (Ca+b) was defined slightly below the ranges presented in the current literature considering the local conditions, which are prone to stress occurrence due to irrigation deficit conditions and damages by insects. As mentioned by several authors, when plants are exposed to stress (water deficit, insect damage, nutrient deficit, etc.), the pigments' concentration is reduced and/or degraded [60][61][62], lowering the energy absorption and thus increasing the reflectance in the visible domain. The observer zenith angle (tto), the azimuth angle (psi), and the sun zenital angle (tts) were set according to the geometry of field spectra readings.
A total number of 100,000 random simulations was used for generating the LUT spectral reflectance data. This simulations' number is considered appropriate for good computational efficiency and high accuracy in the parameter estimation [63][64][65]. The simulated data were aggregated into the eight Sentinel-2 bands (2-8 A) in order to match the FSP_S2 dataset. The PROSAIL simulations were carried out in the software, ARTMO [30].

LUT Inversion for LAI Retrieval
A look-up table (LUT) based inversion technique was applied for the LAI retrieval. This technique is considered easy-to-use and potentially outperforms the limitations of iterative techniques [63,66]. The LUT inversion consists of a direct comparison of simulated spectral data with the observed spectral data (field or satellite data), through one or various cost functions, aiming to identify the best parameter combinations that yield minimal differences between the observed and simulated data [14,67]. The LUT inversion toolbox [68] of the ARTMO software was used to perform this inversion task. The toolbox enables testing simultaneously a wide diversity of cost functions making it possible to optimize the models for different assumptions on the nature and proprieties of the residuals [68,69]. Several classes of cost functions were considered, including minimum contrast estimation, information measures, and M-estimates.
The minimum contrast estimation considers the bidirectional reflectance as a spectral density of the stochastic process and aims at minimizing the distance between a parametric model and a non-parametric spectral density [69]. Within this class, three cost functions were tested, which are detailed in [69,70]: where, D[E, L] is the distance between two functions, E = {e(λi), . . . , e(λn)} is the FSP_S2 or the SP_S2 datasets, L = {l(λi), . . . , l(λn)} is the PROSAIL simulated reflectance (LUT), and, λi, . . . , λn is the spectral bands. The information measures represent measures of distance between two probability distributions. The reflectance is interpreted as a function of probability distribution [68]. In this class, the Bhattacharyya divergence [69] was tested: The M-estimate measures rely on minimizing a dispersion function. Even though the M-estimates are the most commonly used measures, they can be misleading if their basic assumptions of residual normality and absence of extreme values are not met [69]. In this class, we tested the root mean square error (RMSE), which is calculated as: Figure 4 presents the scheme of the LUT generation and the inversion process for LAI retrieval.
The information measures represent measures of distance between two probability distributions. The reflectance is interpreted as a function of probability distribution [68]. In this class, the Bhattacharyya divergence [69] was tested: The M-estimate measures rely on minimizing a dispersion function. Even though the Mestimates are the most commonly used measures, they can be misleading if their basic assumptions of residual normality and absence of extreme values are not met [69]. In this class, we tested the root mean square error (RMSE), which is calculated as: Figure 4 presents the scheme of the LUT generation and the inversion process for LAI retrieval.

LAI Model Calibration, Cross Validation, and External Validation
Data from 2016 were used for the calibration of statistical models (VI and MLRA) because it covered more crop phenological stages and a larger number of observations than the datasets of the other study years (2015 and 2018) ( Table 2), resulting in high data variability.
In the VI approach, the calibration consists in a systematic assessment of all possible band combinations, vegetation indices formulations (Table 3), and curve fitting options. This process yields a table of goodness-of-fit measures for the values estimated by the calibrated model against the measured calibration measures [15]. A validation process is simultaneously performed. In this study, we used a k-fold cross validation method, which involves randomly dividing the dataset into k equalsized sub-datasets. From these k sub-datasets, k − 1 sub-datasets are used as a training dataset and a single k sub-dataset is used as a validation dataset for model testing. The cross-validation process is then repeated k times, with each of the k sub-datasets used as a validation dataset [14]. The results from each of these iterative validation steps are then combined to derive a single estimation value. At the end, all data are used for both training and validation, and each single observation is used for

LAI Model Calibration, Cross Validation, and External Validation
Data from 2016 were used for the calibration of statistical models (VI and MLRA) because it covered more crop phenological stages and a larger number of observations than the datasets of the other study years (2015 and 2018) ( Table 2), resulting in high data variability.
In the VI approach, the calibration consists in a systematic assessment of all possible band combinations, vegetation indices formulations (Table 3), and curve fitting options. This process yields a table of goodness-of-fit measures for the values estimated by the calibrated model against the measured calibration measures [15]. A validation process is simultaneously performed. In this study, we used a k-fold cross validation method, which involves randomly dividing the dataset into k equal-sized sub-datasets. From these k sub-datasets, k − 1 sub-datasets are used as a training dataset and a single k sub-dataset is used as a validation dataset for model testing. The cross-validation process is then repeated k times, with each of the k sub-datasets used as a validation dataset [14]. The results from each of these iterative validation steps are then combined to derive a single estimation value. At the end, all data are used for both training and validation, and each single observation is used for validation exactly once [14,71]. In this study, we used k = 2 and run the calibration/validation process several times until no improvement of the goodness-of-fit statistics were detected in both the calibration and validation processes.
In the MLRA approach, we first used the Gaussian process regression band analysis tool (GPR-BAT) to identify the most sensitive bands to LAI and the minimum number of bands for an acceptable estimate accuracy [17]. Subsequently, we tested nine algorithms using the identified bands for each dataset. We equally used a k = 2 k-fold cross validation process and ran it several times until no improvement of the results were obtained.
For both VI and MLRA, an external validation was performed using a dataset comprising data collected in 2015 and 2018 to assess the model transferability. However, the external validation was not run for the SP_S2 dataset due to the limited number of observations (n = 22).

Model Assumptions, Accuracy Assessment, and Selection
The assumption of residual normal distribution was assessed by the Jarque-Bera (JB) test [72], and the Breusch-Pegan (BP) test [73] was used to investigate the presence of homoscedasticity, testing the dependence of the residuals' variance on the independent variables. For JB tests, the null hypothesis is that the variance of the residuals is normally distributed while for the BP, the null hypothesis is that the residuals are homogeneous. Thus, for p-values less than 5%, we reject the null hypothesis. The models' performance was assessed through three widely used statistical measures: The root mean square error (RMSE), the relative root mean square, and the coefficient of determination (R 2 ). These measures were compared within the different retrieval approaches (VI, MLRA, and LUT inversion) and between them to decide which approach performed better in estimating LAI in the study area. several times until no improvement of the goodness-of-fit statistics were detected in both the calibration and validation processes.

LAI Ground Measurements and Derived from Sentinel-2
In the MLRA approach, we first used the Gaussian process regression band analysis tool (GPR-BAT) to identify the most sensitive bands to LAI and the minimum number of bands for an acceptable estimate accuracy [17]. Subsequently, we tested nine algorithms using the identified bands for each dataset. We equally used a k = 2 k-fold cross validation process and ran it several times until no improvement of the results were obtained.
For both VI and MLRA, an external validation was performed using a dataset comprising data collected in 2015 and 2018 to assess the model transferability. However, the external validation was not run for the SP_S2 dataset due to the limited number of observations (n = 22).

Model Assumptions, Accuracy Assessment, and Selection
The assumption of residual normal distribution was assessed by the Jarque-Bera (JB) test [72], and the Breusch-Pegan (BP) test [73] was used to investigate the presence of homoscedasticity, testing the dependence of the residuals' variance on the independent variables. For JB tests, the null hypothesis is that the variance of the residuals is normally distributed while for the BP, the null hypothesis is that the residuals are homogeneous. Thus, for p-values less than 5%, we reject the null hypothesis. The models' performance was assessed through three widely used statistical measures: The root mean square error (RMSE), the relative root mean square, and the coefficient of determination (R 2 ). These measures were compared within the different retrieval approaches (VI, MLRA, and LUT inversion) and between them to decide which approach performed better in estimating LAI in the study area.  In 2016 the data collection covered the entire crop growing cycle ( Table 2) and thus has an extended range and high average value of LAI, when compared with the data collected in the other In 2016 the data collection covered the entire crop growing cycle ( Table 2) and thus has an extended range and high average value of LAI, when compared with the data collected in the other two years. Also, the diversity between fields (sowing and irrigation conditions and crop treatments; Table 2) as well as the inter-annual weather conditions (Figure 2) likely contributed for the inter-annual LAI variability. The LAI_S2 that comprises data extracted from Sentinel-2 images acquired in 2016 and in 2018 ranged from 0.21 to 1.90 (averaging 0.61) and is very low if compared to the Field_LAI. A correlation between the Field_LAI and LAI_S2 yielded r = 0.69 (n = 30; p-value = 0.0001). This low correlation between measured and Sentinel derived LAI could be explained by several factors, including the scale mismatch and the fact that the LAI_S2 product was not validated in the context of our study area, which has very peculiar characteristics in terms of cropping systems (large in-field heterogeneity) with low LAI and soil background effects (soils with white colour). Also, the temporal mismatch between some LAI_S2 and Field_LAI data, due to cloud cover issues, may have also impacted the correlation between the two LAI products. Nevertheless, the S2 images' acquisition dates are in agreement with the same crop growing stages as those of field data collection dates, minimizing this time effect.

Vegetation Index Based LAI Models and Validation
3.2.1. Field Spectral Data Resampled to 10 nm (FSP_10) The model (FSP_10 vs. Field_LAI) goodness-of-fit statistics of the three best performing VIs in the calibration, cross validation, and external validation dataset are shown in Table 6. The mSR b and mDI d fitted to a second degree polynomial function outperformed the other VI, showing consistent results for the LAI estimations when applied to the calibration, cross validation, and external validation datasets. For the external validation dataset, the mSR b presents R 2 = 0.80, comparatively higher than the mDI d (R 2 = 0.62), while the RMSE = 0.80 of mSR b shows lower accuracy than the mDI d with RMSE = 0.58. The three bands of VI, mSR b, and mDI d include bands within the same spectral regions of the green (565 nm) and red edge (715 nm, 725 nm, and 735 nm), suggesting reliability of the band optimization process for LAI estimation. Figure 6 presents the agreement between the LAI measured and estimated with FSP_10 using the mSR b and the mDI d for the cross validation and external validation datasets. For the cross-validation dataset ( Figure 6 left panel), the two indices exhibit good matching in the whole range of LAI values, but there is an underestimation for values greater than three. For the independent validation dataset ( Figure 6 right panel), the two VIs estimate well the lower values of LAI, but for values above one, the estimation is poorer with an underestimation for the mDI d and an overestimation for the mSR b .  The goodness-of-fit statistics of the three best performing VIs for LAI estimation in calibration, cross validation, and external validation with the FSP_S2 is also presented in Table 6 (FSP_S2 vs. Field_LAI). The mSRc, based on the polynomial model with bands located at the red edge (705 nm, 740 nm) and near infrared (865 nm) region of the spectra, presents better LAI estimates with R 2 ≥ 0.82 and RMSE of 0.40, 0.43, and 0.62, respectively, for the calibration, cross validation, and external  CV-cross validation; EV-external validation; FSP_10-field spectral data resampled to 10 nm; FSP_S2-field spectral data resampled to Sentinel-2 bands; SP_S2-Sentinel-2 spectral data (SP_S2); Field_LAI-field leaf area index; LAI_S2-Sentinel-2 leaf area index product data; mSR-modified simple ratio; mDI-modified difference index; SR-simple ratio; TBSI-three band spectral index; a 0 , a 1 , a 2 -coefficients of the second degree polynomial function (y = a 2 X 2 + a 1 X + a 0 ); m, b-coefficients of the linear function (y = mX + b); k, n-coefficients of the exponential equation (y = ne kx ); m, b-coefficients of the linear function (y = mX + b).

Field Spectral Data Resampled to Sentinel-2 (FSP_S2)
The goodness-of-fit statistics of the three best performing VIs for LAI estimation in calibration, cross validation, and external validation with the FSP_S2 is also presented in Table 6 (FSP_S2 vs. Field_LAI). The mSR c , based on the polynomial model with bands located at the red edge (705 nm, 740 nm) and near infrared (865 nm) region of the spectra, presents better LAI estimates with R 2 ≥ 0.82 and RMSE of 0.40, 0.43, and 0.62, respectively, for the calibration, cross validation, and external validation datasets. The mDI c presents very similar statistics for the calibration and cross validation datasets, but relatively lower R 2 (0.71) and higher RMSE (0.77) for the independent dataset. It is interesting to note that all the validated VIs derived from FSP_S2 were based on bands within the red edge and near infrared region of the spectra, which may suggest consistency of the band optimization and selection process. Figure 7 depicts the agreement between measured and estimated LAI for the cross validation and external validation datasets. In the cross validation, there is a good agreement between measured and estimated LAI for the two indices, especially for LAI values below 2. As with the FSP_10 dataset, there is a miss match trend for LAI values higher than 2, which is more pronounced when the external dataset is used (Figure 7). there is a miss match trend for LAI values higher than 2, which is more pronounced when the external dataset is used (Figure 7).

Sentinel-2 Spectral Data (SP_S2)
The goodness-of-fit statistics of the three best performing VIs for LAI estimation in calibration and cross validation for both SP_S2 vs. Field_LAI and SP_S2 vs. LAI_S2 are also presented in Table  6. For the SP_S2 vs. Field_LAI combination, the selected VIs (TBSIb, mDIc, TBSIc) present very similar statistics and all were constructed with the same spectral bands, centred at red (665 nm), the red edge (783 nm), and near infrared (865 nm). However, the TBSIb fitted to a second degree polynomial function outperforms the other VIs with an RMSE = 0.32, R 2 = 0.79 and RMSE = 0.35, R 2 = 0.74 for calibration and cross validation, respectively.
The SP_S2 vs. LAI_S2 combination presents excellent statistics for the three selected VIs (TBSIb, SR, TBSIa), with the TBSIb the outperforming one with an RMSE = 0.14, R 2 = 0.83 and RMSE = 0.18, R 2 = 0.76, respectively, for calibration and cross validation. This was expected because the SP_S2 was

Sentinel-2 Spectral Data (SP_S2)
The goodness-of-fit statistics of the three best performing VIs for LAI estimation in calibration and cross validation for both SP_S2 vs. Field_LAI and SP_S2 vs. LAI_S2 are also presented in Table 6. For the SP_S2 vs. Field_LAI combination, the selected VIs (TBSI b , mDI c , TBSI c ) present very similar statistics and all were constructed with the same spectral bands, centred at red (665 nm), the red edge (783 nm), and near infrared (865 nm). However, the TBSI b fitted to a second degree polynomial function outperforms the other VIs with an RMSE = 0.32, R 2 = 0.79 and RMSE = 0.35, R 2 = 0.74 for calibration and cross validation, respectively.
The SP_S2 vs. LAI_S2 combination presents excellent statistics for the three selected VIs (TBSI b , SR , TBSI a ), with the TBSI b the outperforming one with an RMSE = 0.14, R 2 = 0.83 and RMSE = 0.18, R 2 = 0.76, respectively, for calibration and cross validation. This was expected because the SP_S2 was used as input data to derive the LAI_S2 products. However, careful interpretation is recommended given the lower correlation between the Field_LAI and the LAI_S2 (r = 0.69). The TBSI b was constructed with green (560 nm), red edge (705 nm), and near infrared (865 nm) bands. Figure 8 depicts the agreement between measured and estimated LAI for the cross validation with the two combinations: SP_S2 vs. Field_LAI (upper panel) and SP_S2 vs. LAI_S2 (lower panel). Clearly, the SP_S2 vs. LAI_S2 combination presents better agreement with the measured and estimated LAI.

Field Spectral Data Resampled to 10 nm (FSP_10)
The LAI estimation based on MLRA was preceded by a band analysis tool in order to select the most sensitive bands among the initial number of 70 bands included in the FSP_10 dataset. The results indicated that increasing the number of bands in the models reduces the LAI estimation accuracy. Indeed, the maximum average R 2 (0.77) was obtained with three band models while the minimum average R 2 (0.70) was acquired with the total number of bands that comprises the hyperspectral dataset (70 bands) ( Figure A1, Appendix A). The three selected bands are centered at 565 nm, 675 nm, and 775 nm, with a 10 nm spectral resolution, corresponding to green, red, and red edge regions of the spectra, respectively. Table 7 presents the goodness-of-fit statistics of the three best performing MLRA for LAI estimation using the three above mentioned bands. The relevance vector machine (RVM) algorithm shows better consistency for LAI estimation in both cross validation and external validation datasets.
With the cross validation dataset, the RVM has an RMSE = 0.54, which is higher than the values

Field Spectral Data Resampled to 10 nm (FSP_10)
The LAI estimation based on MLRA was preceded by a band analysis tool in order to select the most sensitive bands among the initial number of 70 bands included in the FSP_10 dataset. The results indicated that increasing the number of bands in the models reduces the LAI estimation accuracy. Indeed, the maximum average R 2 (0.77) was obtained with three band models while the minimum average R 2 (0.70) was acquired with the total number of bands that comprises the hyperspectral dataset (70 bands) ( Figure A1, Appendix A). The three selected bands are centered at 565 nm, 675 nm, and 775 nm, with a 10 nm spectral resolution, corresponding to green, red, and red edge regions of the spectra, respectively. Table 7 presents the goodness-of-fit statistics of the three best performing MLRA for LAI estimation using the three above mentioned bands. The relevance vector machine (RVM) algorithm shows better consistency for LAI estimation in both cross validation and external validation datasets. Table 7. Goodness-of-fit statistics for the LAI models of the best performing machine learning regression algorithms (MLRA) applied to different combinations of spectral and leaf area index. --CV-cross validation; EV-external validation; FSP_10-field spectral data resampled to 10 nm; FSP_S2-field spectral data resampled to Sentinel-2 bands; SP_S2-Sentinel-2 spectral data (SP_S2); Field_LAI-field leaf area index; LAI_S2-Sentinel-2 leaf area index product data.

MLRA/Type of
With the cross validation dataset, the RVM has an RMSE = 0.54, which is higher than the values of RMSE of other MLRA; nevertheless, the RVM strongly outperform the other algorithms when used in the independent dataset with an RMSE = 0.50, which is clearly better compared to, for example, the SVM (RMSE = 0.67). The very high Gaussian process regression (VHGPR) also presented a good performance for the LAI estimation, both with the cross validation and external validation, although with a higher RMSE (0.63), when compared with the RVM algorithm for the external validation dataset. Figure 9 depicts the comparison between measured and estimated LAI values for cross validation and external validation for the two best performing MLRA with the FSP_10 dataset. For cross validation (Figure 9 left panel), the two algorithms achieved good agreement in the entire LAI range, but there is a slight underestimation for LAI values above two. For the external validation dataset (Figure 9 right panel), the two algorithms exhibit good estimation in the full range of LAI values, although the VHGPR has lower accuracy for values exceeding two. CV-cross validation; EV-external validation; FSP_10-field spectral data resampled to 10 nm; FSP_S2-field spectral data resampled to Sentinel-2 bands; SP_S2-Sentinel-2 spectral data (SP_S2); Field_LAI-field leaf area index; LAI_S2-Sentinel-2 leaf area index product data. Figure 9. Comparisons of the LAI measured and estimated with the two best performing machine learning regression algorithms (MLRA) using the field spectra resampled to 10 nm (FSP_10). CVcross validation, EV-external validation. RVM-relevance vector machine; VH-GPR-VH. Gaussian

Field Spectral Data Resampled to Sentinel-2 Bands (FSP_S2)
A band analysis was also performed prior to the MLRA assessment to select the best bands of the eight comprising the FSP_S2, for the LAI estimation. The results indicated that the LAI estimation accuracy increases from a minimum average of R 2 = 0.51, when only one band is involved, to a maximum average of R 2 = 0.77 when three bands are involved in the modelling ( Figure A2, Appendix A). There was no accuracy improvement with additional band inclusion in the models. The selected bands match to the red and red edge Sentinel-2 bands, centered at 665 nm, 705 nm, and 783 nm.
The goodness-of-fit statistics of the three best performing MLRA when using the FSP_S2 dataset for LAI estimation are presented in Table 7. As with FSP_10, the RVM algorithm shows better consistency for both cross validation and external validation datasets. In the cross validation, the RVM has an RMSE = 0.52 and R 2 = 0.73, which is slightly poorer than the SVM values of RMSE = 0.48 and R 2 = 0.78. However, in the independent dataset, the RVM presents higher accuracy with RMSE = 0.53, compared to the SVM values of RMSE = 0.90. The bagging trees (BaT) algorithm is equally consistent, but with relatively higher RMSE (0.63 and 0.64) for the cross validation and external validation datasets, respectively. Figure 10 presents the comparison between LAI measured and estimated for both cross validation and external validation datasets of the two best performing MLRA for the FSP_S2 dataset. For the cross validation ( Figure 10 left panel), the two algorithms reveal good agreement between the measured and estimated LAI throughout the whole range, but the estimation accuracy decreases as the LAI values increase, resulting in a spread-out pattern for LAI values above two. For the external validation, the two algorithms present an overestimation trend in the full LAI range. validation and external validation datasets of the two best performing MLRA for the FSP_S2 dataset. For the cross validation ( Figure 10 left panel), the two algorithms reveal good agreement between the measured and estimated LAI throughout the whole range, but the estimation accuracy decreases as the LAI values increase, resulting in a spread-out pattern for LAI values above two. For the external validation, the two algorithms present an overestimation trend in the full LAI range.

Sentinel-2 Spectral Data (SP_S2)
A band analysis performed preceding the MLRA assessment revealed that two bands centered at blue (490 nm) and red (665 nm) are the most appropriate for LAI estimation using the SP_S2 dataset. Inclusion of more bands decreased the estimation accuracy ( Figure A3, Appendix A).
(RMSE = 0.23, R 2 = 0.61), and SVM (RMSE = 0.23, R 2 = 0.60). Figure 11 presents the comparison between LAI measured and estimated for cross validation for the two best performing MLRA with the combinations, SP_S2 vs. Field_LAI (upper panel) and SP_S2 vs. LAI_S2 (lower panel). For the SP_S2 vs. Field_LAI, the SVM algorithm shows better agreement while the RFF presents the best agreement for the SP_S2 vs. LAI_S2. Figure 11. Comparisons of the LAI measured and estimated with the best performing machine learning regression algorithms (MLRA) using the Sentinel-2 spectral data (SP_S2) with the field leaf area index (Field_LAI) and the Sentinel-2 leaf area index product (LAI_S2). GPR-Gaussian process regression; SVM-support vector regression; RFF-random forest Fitensemble; RVM-relevance vector machine. The full line depicts the linear regression and the dashed line the linear regression through the origin.

LUT Inversion Based LAI Estimation and Validation
In this section, we present the goodness-of-fit statistics of LAI estimation by the LUT inversion approach using three different combinations of spectral and LAI datasets: (i) FSP_S2 and the field LAI (FSP_S2 vs. Field LAI); (ii) spectral data of Sentinel-2 images and the LAI product derived from Sentinel-2 (SP_S2 vs. LAI_S2); and (iii) SP_S2 combined with Field LAI (SP_S2 vs. Field_LAI). Table 8 presents the goodness-of-fit statistics of LAI retrieval through the LUT inversion using the three combinations of datasets. For the three datasets, three cost functions, namely the K(x) = (log(x)) 2 , Figure 11. Comparisons of the LAI measured and estimated with the best performing machine learning regression algorithms (MLRA) using the Sentinel-2 spectral data (SP_S2) with the field leaf area index (Field_LAI) and the Sentinel-2 leaf area index product (LAI_S2). GPR-Gaussian process regression; SVM-support vector regression; RFF-random forest Fitensemble; RVM-relevance vector machine. The full line depicts the linear regression and the dashed line the linear regression through the origin.

LUT Inversion Based LAI Estimation and Validation
In this section, we present the goodness-of-fit statistics of LAI estimation by the LUT inversion approach using three different combinations of spectral and LAI datasets: (i) FSP_S2 and the field LAI (FSP_S2 vs. Field LAI); (ii) spectral data of Sentinel-2 images and the LAI product derived from Sentinel-2 (SP_S2 vs. LAI_S2); and (iii) SP_S2 combined with Field LAI (SP_S2 vs. Field_LAI). Table 8 presents the goodness-of-fit statistics of LAI retrieval through the LUT inversion using the three combinations of datasets. For the three datasets, three cost functions, namely the K(x) = (log(x)) 2 , the K(x) = x(log(x)) -x, and the Bhattacharyya divergence, evidenced good LAI estimation. Nevertheless, the highest LAI estimation accuracy was achieved applying the Bhattacharyya divergence cost function to the SP_S2 vs. LAI_S2 dataset, resulting in an RMSE = 0.20 and R 2 = 0.70. Furthermore, this dataset combination showed better performance than the others with all the evaluated cost function. However, the SP_S2 vs. LAI_S2 dataset produced relatively lower association measure values, 0.65 ≤ R 2 ≤ 0.72, compared to the range of 0.82 ≤ R 2 ≤ 0.86 for FSP_S2 vs. Field LAI. Figure 12 illustrates the agreement between the measured and estimated LAI values with the two best performing cost functions and using the three datasets. With the FSP_S2 vs. Field LAI dataset, the two cost functions present a strong underestimation of the LAI for values higher than one, while the estimation for the LAI values below one is more accurate (Figure 11 first two plots). The underestimation is also evident when the SP_S2 vs. Field LAI is used with all the cost functions ( Figure 12, third and fourth plots). On the other hand, for the SP_S2 vs. LAI_S2, the two cost functions accurately estimate the LAI throughout the entire range ( Figure 12, last two plots). Table 8. Goodness-of-fit statistics of look-up-table (LUT) inversion using the field spectral data aggregated according to Sentinel-2 bands and the field LAI (FSP_S2 vs. Field LAI); the spectral data of Sentinel-2 images and the LAI product derived from Sentinel-2 (SP_S2 vs. LAI_S2); the spectral data of Sentinel-2 images and the field LAI (SP_S2 vs. Field LAI), while considering different cost function (CF) algorithms. aggregated according to Sentinel-2 bands and the field LAI (FSP_S2 vs. Field LAI); the spectral data of Sentinel-2 images and the LAI product derived from Sentinel-2 (SP_S2 vs. LAI_S2); the spectral data of Sentinel-2 images and the field LAI (SP_S2 vs. Field LAI), while considering different cost function (CF) algorithms.  Figure 12 illustrates the agreement between the measured and estimated LAI values with the two best performing cost functions and using the three datasets. With the FSP_S2 vs. Field LAI dataset, the two cost functions present a strong underestimation of the LAI for values higher than one, while the estimation for the LAI values below one is more accurate (Figure 11 first two plots). The underestimation is also evident when the SP_S2 vs. Field LAI is used with all the cost functions ( Figure 12, third and fourth plots). On the other hand, for the SP_S2 vs. LAI_S2, the two cost functions accurately estimate the LAI throughout the entire range ( Figure 12, last two plots). Figure 12. Comparisons of the measured and estimated LAI of the best LUT inversion cost functions using the field spectral data resampled to Sentinel-2 bands with the field leaf area index (FSP_S2 vs. F_LAI-first two plots), sentinel-2 spectral data with field leaf area index data (SP_S2 vs. F_LAI- Figure 12. Comparisons of the measured and estimated LAI of the best LUT inversion cost functions using the field spectral data resampled to Sentinel-2 bands with the field leaf area index (FSP_S2 vs. F_LAI-first two plots), Sentinel-2 spectral data with field leaf area index data (SP_S2 vs. F_LAI-third and fourth plots), and Sentinel-2 spectral data with Sentinel-2 leaf area index product data (SP_S2 vs. LAI_S2-last two plots). The full line depicts the linear regression and the dashed line the linear regression through origin. The Breusch-Pegan (BP) test p-value is reasonably high for almost all the datasets and modelling approaches, excluding the MLRA (BaT), for FSP_S2 and VI (mNDc and TBSIc) for SP_S2 vs. Field_LAI, all for cross validation, confirming the presence of homoscedasticity. The residuals' normality was confirmed from the Jarque-Bera (JB) test, showing higher values than p-values for all the developed models. The results of BP and JB tests are included in Table A2 (Appendix A).
Generally, there was no substantial difference in LAI estimation accuracy between the two types of datasets, hyperspectral (FSP_10) and multispectral (FSP_S2) for both VIs and MLRA approaches. In fact, the advantage of hyperspectral data over multispectral data in estimating LAI remains a matter of debate [78]. Theoretically, the high spectral resolution of hyperspectral data disclose the spectral details obscured with multispectral data for LAI estimation [79]. However, the LAI insensitive bands included in hyperspectral data require additional computational time and distort the accuracy of LAI retrieval, with this the reason why there is a need for dimensionality reduction of hyperspectral data [78].
The performance of LUT inversion for LAI retrieval is in accordance with the findings of other studies for maize using the same and other retrieval methods, with RMSE values of 0.40-0.43 [2], 0.46-1.21 [59], 0.73 [80], 0.41-0.76 [81], and 0.63 [82]. The underestimation trend observed while using the FSP_S2 dataset for the inversion was also found in maize by [2]. The reason could be related with the row planting pattern of maize, which diverges from the turbid assumption of the PROSAIL model [83]. In the case of the present study, this could be exacerbated by the low input cropping system of the field area and important differences in the planting geometry (Table 1), which resulted in high heterogeneity within and between the sampled fields.
Concerning the applied cost functions, our findings show that the commonly used RMSE was not the best performing cost function, but instead the contrast function, (K(x) = (log(x)) 2 ), when inverting through the field spectral data (FSP_S2), and the Bhattacharyya divergence cost function, when using the SP_S2 dataset for the inversion. Similar findings were reported by [68] for LAI. Table 9 summarizes the RMSE and slope (b) of the best models identified in the different combinations of spectral and LAI data and modelling approaches tested.
For all the modelling approaches, the combination, SP_S2 vs. LAI_S2, yielded the best accuracy, which was, in fact, expected because the SP_S2 is an input for the derivation of LAI_S2 products. The comparison of combinations involving Field_LAI shows that the statistical approach based in VIs using the SP_S2 vs. Field_LAI yielded the most accurate LAI estimation (RMSE = 0.35 and b = 0.82), outperforming the physically-based approach of LUT inversion that is often considered more robust. These results may be due to the small geographical scale (field level) of our study and involving a single crop type. According to [78], statistical predictive LAI models are based on numerical relationships, and thus rely greatly on the specific location, including the crop's condition and soil background reflectance, therefore, they are more suitable for small-scale studies. Additionally, the parameterization of PROSAIL for the application of LUT inversion may have been hindered by the large heterogeneity within and between crop fields in the study area. FSP_10-field spectral data resampled to 10 nm; FSP_S2-Field spectral data resampled to Sentinel-2 bands; Sentinel-2 spectral data (SP_S2); Field_LAI-field leaf area index; LAI_S2-Sentinel-2 leaf area index product data; VI-vegetation index; MLRA-machine learning regression algorithm; * b coefficient of the linear regression between measured and predicted LAI.

Conclusions
In this paper, we successfully calibrated and validated models for maize LAI estimation in low-input crop systems based on statistical and physical approaches. Hyperspectral and multispectral data, obtained both from field and satellite sensors, were tested for the LAI retrieval and further compared with field and satellite LAI data. The robustness of the developed models was indicated by the consistency of the selected electromagnetic band regions, whichever the modelling approach or dataset combination was applied. Additionally, the models reasonably copedwith transferability issues by adjusting to relatively different cropping systems (Fields 1-4) and different weather conditions (2015,2016,2018), as suggested by the results of the external validation process. The most accurate model involves the TBSI b spectral index. This VI is built with three bands centered in the red, red edge, and near infrared regions of the electromagnetic spectrum with widely known biophysical significance for LAI estimation. The hyperspectral data (aggregated to 10 nm) failed to improve the LAI estimation accuracy comparatively to Sentinel-2 multispectral data. This finding is of particular relevance for the operational application of spectral data in crop monitoring, though Sentinel-2 data is freely available and presents good spatial and spectral resolutions. However, future research should consider using field LAI data acquired with high precision equipment, including other crop types and extensive sampling, in order to increase the ground truth data and, as a result, improve the accuracy of LAI retrieval.

Conflicts of Interest:
The authors declare no conflict of interest.
Appendix A Figure A1. Variation of the models' R 2 as a function of the number of bands involved in machine learning modelling using the field spectral data resampled to 10 nm (FSP_10). The highest R2 was achieved using three bands centered at 565 nm, 675 nm, and 775 nm. Figure A2. Variation of the models' R 2 as a function of the number of bands involved in machine learning modelling using the field spectral data resampled to Sentinel-2 bands (FSP_S2). The highest R2 was achieved using three bands: Centered at 665 nm, 705 nm, and 783 nm.

Conflicts of Interest:
The authors declare no conflict of interest.
Appendix A Figure A1. Variation of the models' R 2 as a function of the number of bands involved in machine learning modelling using the field spectral data resampled to 10 nm (FSP_10). The highest R2 was achieved using three bands centered at 565 nm, 675 nm, and 775 nm. Figure A2. Variation of the models' R 2 as a function of the number of bands involved in machine learning modelling using the field spectral data resampled to Sentinel-2 bands (FSP_S2). The highest R2 was achieved using three bands: Centered at 665 nm, 705 nm, and 783 nm. Yielded R 2 Number of bands Figure A2. Variation of the models' R 2 as a function of the number of bands involved in machine learning modelling using the field spectral data resampled to Sentinel-2 bands (FSP_S2). The highest R 2 was achieved using three bands: Centered at 665 nm, 705 nm, and 783 nm.
Remote Sens. 2018, 10, x FOR PEER REVIEW 25 of 31 Figure A3. Variation of the models' R 2 as a function of the number of bands involved in machine learning modelling using the spectral data from Sentinel-2 images (SP_S2). The highest R2 was achieved using two bands: Centered at 490 nm and 665 nm.  Figure A3. Variation of the models' R 2 as a function of the number of bands involved in machine learning modelling using the spectral data from Sentinel-2 images (SP_S2). The highest R 2 was achieved using two bands: Centered at 490 nm and 665 nm.