Sensitivity of Common Vegetation Indices to the Canopy Structure of Field Crops

Leaf inclination angle distribution is an important canopy structure characteristic which directly impacts the fraction of the intercepted solar radiation. Together with the leaf area index (LAI) it determines the structure and fractional cover of a homogeneous crop canopy. Unfortunately, this key canopy parameter has usually been ignored when applying common multispectral vegetation indices to the mapping of LAI, although its impact is known from model simulations. Therefore, we measured leaf angles and determined their distribution (quantified using the leaf mean tilt angle, MTA) for six crop species with different structures growing on 162 plots with a broad range of LAI (1.1–5.0) and leaf chlorophyll content (26–94 μg cm−2). Next, we calculated six vegetation indices widely used for LAI monitoring—the normalized difference vegetation index (NDVI), the enhanced vegetation index (EVI), the two band enhanced vegetation index (EVI2), the modified triangular vegetation index (MTVI2), the optimized soil adjusted vegetation index (OSAVI) and the wide dynamic range vegetation index (WDRVI)—from airborne imaging spectroscopy data. We calculated the Spearman’s correlation coefficient Rs, a non-parametric statistic chosen because of the non-normal distribution of canopy parameters. All studied indices depended on the LAI (0.50 ≤ Rs ≤ 0.71), but the dependence on the MTA was of similar magnitude (−0.83 ≤ Rs ≤ −0.53) with EVI, EVI2, OSAVI and MTVI2 depending more strongly on MTA than on LAI. All studied indices were good proxies (0.78 ≤ Rs ≤ 0.88) for vegetation fractional cover (Fcover) which, for homogeneous crop canopies, is a nonlinear function of LAI and MTA. EVI2 and MTVI2 were the most strongly correlated with Fcover, although the difference to the other studied indices was small. This first study involving a large range of crop structures confirms the results from canopy reflectance simulations and emphasizes the necessity of leaf angle information for the successful mapping of LAI with Earth observation data.


Introduction
Canopy structure, or its architecture, is understood as the spatial distribution and orientation of vegetation components within a canopy.It directly affects the interception of solar radiation and thus also how radiation is reflected or intercepted by the canopy.Vegetation structure has an effect on land surface energy balance and it influences leaf temperature, photosynthetic activity and plant growth [1][2][3][4].Vegetation structure thus plays a dual role in remote sensing.First, it is a canopy characteristic which needs to be retrieved to understand biophysical processes.Second, due to its direct effect on canopy reflectance, it needs to be taken into account-and corrected for-when other canopy biophysical or biochemical variables are determined from remote sensing data.
The complexity of canopy structure varies with biome and plant functional type.Crop canopy is usually assumed to be the simplest possible case-horizontally homogeneous and with a constant distribution of leaf normals, i.e., the vertical profile of the canopy is the same at all locations.Its structure can thus be described by the amount of leaves above the unit horizontal ground area and their directional distribution, which can be quantified using the leaf area index (LAI) and leaf angle distribution (LAD), respectively.LAI is a key variable for understanding the biophysical characteristics of the vegetated areas of the planet [5][6][7].For agricultural management, crop LAI measured over large areas can be used to estimate productivity on national to global scales.The other structural characteristic, LAD, quantifies the distribution of the angles between leaf normals and the vertical direction.If the azimuth angles of leaves are uniformly distributed, LAD and LAI completely determine the transmission of radiation through the canopy [8].Compared with LAI, much less attention has been paid to LAD when describing the biophysical characteristics of vegetation.
Spectral vegetation indices (VIs), calculated from canopy reflectance in a small number of selected spectral bands mostly in the visible and near infrared (NIR) spectral domains [9][10][11][12], are commonly used for the remote estimation of vegetation biophysical characteristics.Visible bands are utilized for determining the absorption by plant pigments active in this spectral region, whereas the reflectance in the NIR domain, where green leaves have minimal scattering, contains information on canopy cover and thickness.However, in addition to the number of leaves, canopy reflectance is also affected by many factors, such as soil background, atmospheric effects, canopy element optical properties and canopy architecture.
Despite its central role, there is a noticeable lack of empirical evidence on the effect of LAD on even the most common VIs.The main reason is likely a lack of in-situ measurements of canopy LAD, or the central moment of the distribution, the mean tilt angle (MTA).A traditional direct approach for LAD measurement is using clinometers in contact with the leaf surface [13].This method is laborious and time consuming, and direct measurement is inevitable to disturb the leaf natural orientation.Indirect leaf mean tilt angle estimation relies on the inversion of the gap fraction below the canopy [14].Some specialized techniques were applied to characterize 3D canopy structure, e.g., portable laser canning [15,16], which were resource demanding and rigid to operate.A digital photography method was developed to estimate the LAD of tree species [17][18][19] and extended and validated for field crops [20].
The effect of MTA is obvious across the whole optical spectrum as this parameter affects canopy transmittance in both the view and illumination directions and thus the brightness and visibility of the vegetation background.An especially strong dependence of reflectance on MTA has been found in the far-red edge spectral region, around 750 nm [21,22].However, this spectral region is not included in most satellite instruments for which data is commonly available, and hence also in the most widely used VIs.We therefore point our focus on the effects of LAD on the most common wavelength combinations in the visible spectrum and NIR used in vegetation remote sensing.
A number of VIs have been developed for estimating LAI by establishing statistical relationships between VIs and measured LAI for various vegetation types.The most obvious first choice is the normalized difference vegetation index (NDVI) [23] which reflects the fractional absorbed photosynthetically active radiation and exhibits a strong relationship with vegetation characteristics [24][25][26].This classic VI based on reflectance in red and NIR bands has numerous modifications aimed at improving robustness under varying atmospheric and background conditions, as well as linearity with specified biophysical variables.We included in this study the enhanced vegetation index (EVI) [27], which includes also reflectance in the blue band, the two-band version of EVI (EVI2, [28]), OSAVI [29], the modified triangular vegetation index (MTVI2) [30], which also includes a green band, and WDRVI [31,32].All these indices have reported dependencies on MTA at low to middle LAI [33], based on simulated reflectance spectra.To our knowledge, the relationships have not been tested using empirical data covering a wide range of canopy structures.
The purpose of our study was to determine the role of two key canopy structural variables, MTA and LAI, in determining the values of commonly used and robust vegetation indices under actual field conditions.We measured the leaf angle distributions and reflectance properties of 162 study plots growing six different species with a large variation in both LAI and MTA.To generalize and validate the accuracy of our field-measured data, we also applied a common canopy reflectance model with inputs taken from spectral and structural measurements at the test site.

in the reference [20]).
Canopy leaf-angle distribution was determined using the photographic approach [17], validated and extended to field crops [20] on 6 July 2012 in the same growth stage as during the spectroscopy measurement.Canopy photographs were taken using a Nikon D1X digital camera mounted on a tripod and leveled with bubble level outside of the plots approximately 1 m away from the edge of the field.The camera height was between 30 and 50 cm depending on crop height.Five to six photographs were taken for each species from one or several plots.Finally, the inclination angles of 75-100 leaves or leaf fragments oriented perpendicularly to the viewing direction were measured using ImageJ 1.47 software, which is an open source image processing program [34], and used for the calculation of the leaf mean tilt angle (MTA).We assumed LAD to be species-specific, as it has been shown to vary more among species than within species [8,13,35].A more detailed description of the test site and the photographic method was given by Zou et al. [20].
Airborne spectroscopy imagery covering the study area was acquired using an AISA Eagle II imaging spectrometer (Specim Ltd., Oulu, Finland) on 25 July 2011.The data were collected between 09:36 and 10:00 local time at approximately 600 m above ground.The flight direction was parallel with the sunrays.The solar zenith angle varied between 50.4 • and 48.1 • and the measurement zenith angle between 0 and 18.9 • was determined by the swath width of the sensor and aircraft orientation changes.Reflected radiation was measured in 64 channels in the visible and NIR spectrum (400-1000 nm) with a spectral resolution of 9-10 nm and a spatial resolution of 0.4 m.The spectroscopic imagery was radiometrically calibrated using the CaliGeo 4.9.9 software package (Specim Ltd., Oulu, Finland) and georectifed using Parge (ReSe Applications Schläpfer, Wil, Switzerland).Finally, the radiometric data were converted to top-of-canopy hemispherical-directional reflectance factors using ATCOR-4 (ReSe Applications Schläpfer, Wil, Switzerland).The spectroscopic data were extracted from the imagery and averaged as average plot spectra.
A SunScan ceptometer (Delta-T Devices, Cambridge, UK) was used to measure LAI in the field.The data used in this study were collected within five days of the airborne campaign.The SunScan ceptometer bar was entered from the edge of the plot at approximately 45 • to the row direction to minimize the row effects on LAI measurement as suggested in the instrument manual.LAI was calculated from the measured transmittance following the model implemented in the SunScan hardware.The possibility of light transmitted through the canopy at polar angle θ can be modeled for a random spatial distribution element canopy using an exponential model: where a is the leaf absorptivity for light set to 0.85 in this study (SunScan SSI user manual) and k(θ, χ) is the canopy directional light extinction coefficient depending on LAD.The canopy transmittance model used by SunScan uses an ellipsoidal LAD model depending on the single parameter χ to characterize LAD.Commonly, χ is assumed to have equal unity corresponding to random distribution of leaf normals due to a lack of information on the actual LAD.In this study, we calculated χ from the measured MTA using the following function derived from Equation ( 16) in [36]: For an ellipsoidal leaf angle distribution model, k(θ, χ) can now be calculated as [13]: The calculations were performed as follows.First, we determined MTA from the measured leaf angles.Next, we calculated the species-specific ellipsoidal leaf angle distribution parameter χ from Equation ( 2) and two extinction coefficient values from Equation (3): k(θ S , χ) corresponding to the solar zenith angle θ S at the time of the SunScan measurements, and k(0, χ) for the nadir direction.To obtain the LAI, we solved Equation (1) for LAI and used the SunScan-measured P(θ S ) and k(θ S , χ).Finally, the LAI for each plot was averaged from 4-5 sequential SunScan readings.We used the LAI value for each plot to calculate the fraction of the canopy cover as: where the light transmission in the nadir direction, P(0), was calculated from Equation ( 2) with θ = 0. Fcover can also be directly related to SunScan measurements.After making all the substitutions listed above, we can write an equation relating Fcover and SunScan-measured transmittance measurements: A portable spectroradiometer (ASD Handheld, Boulder, CO, USA) with spectral range between 325 and 1075 nm and a spectral resolution of 3 nm was used to measure the bare soil reflectance in four sample plots in a harvested area on 7 October 2011.Bare soil reflected radiation was collected at a distance approximately 40 cm above the soil surface and a field of view of 25 • .Soil reflectance was calculated as the ratio of soil-reflected radiance to that measured above a calibrated white Spectralon panel (Labsphere Inc., Sutton, NH, USA).Due to the different solar illumination conditions between soil spectral and airborne spectral measurements, the measured soil reflectance was corrected with a bidirectional reflectance distribution function (BRDF) model developed by Walthall et al. [37] and modified by Nilson and Kuusk [38].The soil spectrum used for PROSAIL model simulation is presented in Figure 1.

Vegetation Indices
The VIs (Table 1) used in the study (NDVI, EVI, EVI2, MTVI2 and OSAVI) make use of a total of four spectral bands (blue, green, red, and NIR) available in most modern satellite remote sensing instruments.The corresponding wavelengths used in the study are blue: 470 nm, green: 550 nm, red: 680 nm and NIR: 800 nm.As these exact wavelengths were not available in the AISA data, we used the closest available wavelength.The actual wavelengths used in the analysis were 470 nm, 551 nm, 682 nm, and 795 nm, respectively.

PROSAIL Simulations
We studied the sensitivity of the VIs to canopy structure (MTA and LAI) using the PROSAIL canopy reflectance model.The model is a combination of the leaf optical reflectance model PROSPECT-5 [39] and the canopy reflectance model SAILH [40,41].The PROSPECT-5 model simulates leaf reflectance and transmittance from six input parameters: leaf chlorophyll a and b content ( ), leaf carotenoid content ( ), leaf dry matter content ( ), leaf water content ( ), leaf brown pigment content ( ), and the mesophyll structure parameter .The SAILH sub-model simulates canopy reflectance using LAI, MTA, hot spot size parameter, soil reflectance, solar zenith angle, sensor view angle, sensor azimuth angle and the fraction of incident diffuse sky radiation.The leaf inclination angle distribution is characterized from the leaf mean tilt angle to generate an ellipsoidal leaf angle distribution [42].

Vegetation Indices
The VIs (Table 1) used in the study (NDVI, EVI, EVI2, MTVI2 and OSAVI) make use of a total of four spectral bands (blue, green, red, and NIR) available in most modern satellite remote sensing instruments.The corresponding wavelengths used in the study are blue: 470 nm, green: 550 nm, red: 680 nm and NIR: 800 nm.As these exact wavelengths were not available in the AISA data, we used the closest available wavelength.The actual wavelengths used in the analysis were 470 nm, 551 nm, 682 nm, and 795 nm, respectively.

Index
Formula References

PROSAIL Simulations
We studied the sensitivity of the VIs to canopy structure (MTA and LAI) using the PROSAIL canopy reflectance model.The model is a combination of the leaf optical reflectance model PROSPECT-5 [39] and the canopy reflectance model SAILH [40,41].The PROSPECT-5 model simulates leaf reflectance and transmittance from six input parameters: leaf chlorophyll a and b content (C ab ), leaf carotenoid content (C Car ), leaf dry matter content (C m ), leaf water content (C w ), leaf brown pigment content (C bp ), and the mesophyll structure parameter N. The SAILH sub-model simulates canopy reflectance using LAI, MTA, hot spot size parameter, soil reflectance, solar zenith angle, sensor view angle, sensor azimuth angle and the fraction of incident diffuse sky radiation.The leaf inclination angle distribution is characterized from the leaf mean tilt angle to generate an ellipsoidal leaf angle distribution [42].
First, to illustrate the effects of LAI and MTA on VIs, we varied the MTA between 15 and 75 degrees, set the LAI value to 1, 3 and 5, and set all other parameters to values characteristic of healthy green vegetation: C ab = 40 µg cm −2 , C Car = 8 µg cm −2 , C w = 0.001 − 0.02 cm.Other parameters were fixed to values taken from the literature: C m = 0.005 g cm −2 (the average value of the six species [43][44][45][46]), N = 1.55 (the average value of various crops [27]), C bp = 0 (assuming the leaves were green during the measurements), and the hot spot parameter was set to 0.01 (this parameter has no practical influence on the results as the view direction was far from the hot spot because of low sun).The fraction of diffuse radiation was calculated using the 6S atmosphere radiative transfer model [47] from actual observation geometry and nearby sun photometer measurements.The parameters describing the illumination and view geometry in the model were set to coincide with the airborne spectroscopy imagery acquisition: solar zenith angle 49.4 • , sensor zenith angle 9 • , and azimuth angle 90 • .
Next, we generated 100,000 parameter combinations by uniformly sampling the variation ranges of C ab , C w , LAI and MTA within the measured (for field-measured parameters) or natural variation range (for parameters for which this information was available).According to measurements [20], C ab varied between 25 and 100 µg cm −2 , LAI between one and five, and MTA between 15 and 70 degrees (MTA = 15, 30, 40, 50, 60 and 70 • ).C w varied in a wide natural range of the parameter, between 0.001 and 0.020 cm; C car was linked to C ab as 1:5 based on the LOPEX93 dataset [48].

Statistical Analysis
We analyzed the correlation between VIs and canopy parameters (LAI, MTA and Fcover) for model simulations and empirical data using the Spearman's correlation coefficient R s (−1 ≤ R s ≤ 1) which is a statistical measure of the strength of the relationship between paired data.Its interpretation is similar to that of Pearson's correlation coefficient: the closer the R s is to ±1 the stronger the correlation.Unlike Pearson's R, there is no requirement of normality and assumption of a linear relationship between the two variables.

Results
In the modeled data for typical healthy green vegetation, all indices showed a clear and mostly negative dependence on MTA (Figure 2).The curves were most similar for the lowest LAI value, LAI = 1.At medium LAI values (LAI = 3) of the three LAI levels, the indices showed a smaller dependence on MTA.Also, the curves in Figure 1 become nonlinear and indicate lower dependence on MTA at low MTA values.The latter is especially evident for NDVI (Figure 2a) and WDRVI (Figure 2f).At LAI = 5, the nonlinearity was even larger compared with lower LAI values.The dependence of NDVI and WDRVI on MTA was slightly positive for the smallest MTA values.However, the dependence was very small and for practical purposes, NDVI and WDRVI were insensitive to the variation of MTA between below MTA = 50 • for large LAI values.The other indices, EVI, EVI2, OSAVI and MTVI2, were still highly sensitive to MTA (R s > 0.6) even at LAI = 5.
All indices showed a monotonic dependence on LAI when all other inputs in the PROSAIL model (see Section 2.3) were fixed to the values of green healthy vegetation (Figure 2).The dependence, measured as the distance between the lines in Figure 1, generally decreased with LAI and MTA.For NDVI, the lines for LAI = 3 and LAI = 5 nearly intersect at MTA = 15 • indicating a lack of dependence of the index on LAI.This well-known saturation with LAI, however, is not present in EVI, EVI2 and MTVI2, for the LAI intervals investigated here.
In the 100,000 simulations (Figures S1-S3) and empirical data, the results visualized in Figure 2 were corroborated: VIs showed generally stronger correlations with MTA than LAI (Table 2).The only exceptions were NDVI and WDRVI, for which the correlation with LAI was stronger than with MTA in both datasets; for NDVI, the correlation with LAI was stronger than with MTA in the field-measured data.However, even the strongest correlations with LAI, (simulated dataset: R s = 0.71 and measured dataset: R s = 0.58) for NDVI and WDRVI were weaker than the strongest correlation (simulated dataset: R s < −0.81 and measured dataset: R s < −0.71) for EVI, EVI2 and MTVI2.The strongest correlations in both the simulated and measured datasets were found between all studied indices and Fcover (Table 2).In the simulated dataset, was consistently above 0.93 with = 0.95, indicating a near-perfect correlation for NDVI and WRDVI.In the measured data, the values of were somewhat smaller, around 0.8, with the lowest value ( = 0.78) produced by NDVI and WRDVI.
The correlations in the empirical dataset are presented graphically in Figures 3-5.The species-specific MTA values cause the clustering of points into vertical stripes in these figures.The length of each stripe depends on the variation of a VI on parameters other than MTA.The variation of these parameters in the measured dataset is provided in Table 1 in the reference [20].Despite considerable variation, the strong effect of MTA on EVI, EVI2, OSAVI and MTVI2 is clear from Figure 3b-e  The strongest correlations in both the simulated and measured datasets were found between all studied indices and Fcover (Table 2).In the simulated dataset, R s was consistently above 0.93 with R s = 0.95, indicating a near-perfect correlation for NDVI and WRDVI.In the measured data, the values of R s were somewhat smaller, around 0.8, with the lowest value (R s = 0.78) produced by NDVI and WRDVI.
The correlations in the empirical dataset are presented graphically in Figures 3-5.The species-specific MTA values cause the clustering of points into vertical stripes in these figures.The length of each stripe depends on the variation of a VI on parameters other than MTA.The variation of these parameters in the measured dataset is provided in Table 1   The correlation between the VIs and LAI were all positive (Figure 4).No obvious saturation at high LAI values (LAI > 3) was visible for any of the indices used in this study.However, the point clouds were divided into two subclouds with a small set of measurements with LAI = 2-3 (corresponding to faba bean and lupin, the two non-cereal crops with the lowest recorded MTAs) with VI values higher than the rest of the data.The correlation between the VIs and LAI were all positive (Figure 4).No obvious saturation at high LAI values (LAI > 3) was visible for any of the indices used in this study.However, the point clouds were divided into two subclouds with a small set of measurements with LAI = 2-3 (corresponding to faba bean and lupin, the two non-cereal crops with the lowest recorded MTAs) with VI values higher than the rest of the data.The large R s values between the VIs and Fcover in Table 2 are accompanied by the close to linear point clouds with little scatter in Figure 5. Especially the VIs with the largest R s values (EVI, EVI2, OSAVI, MTVI2) are also characterized by the smallest scatter; however, the performance differences between the indices are small.A small non-linearity with Fcover can be visually detected for EVI, EVI2 and MTVI2 (Figure 5b-e), but not for OSAVI (Figure 5d).The large values between the VIs and Fcover in Table 2 are accompanied by the close to linear point clouds with little scatter in Figure 5. Especially the VIs with the largest values (EVI, EVI2, OSAVI, MTVI2) are also characterized by the smallest scatter; however, the performance differences between the indices are small.A small non-linearity with Fcover can be visually detected for EVI, EVI2 and MTVI2 (Figure 5b-e), but not for OSAVI (Figure 5d).

Discussion
Clearly, the observed correlation with the leaf inclination angle in the field crop canopy (quantified by MTA) cannot be ignored when estimating other canopy biophysical parameters.Although it is possible that the spectral effect is caused by some other parameter strongly correlated with MTA, no other parameter has been identified in earlier modeling studies [29].Thus far, the influence of LAD on VIs has not been adequately investigated with empirical data.In this study, we attempted to bridge this gap for six commonly used multispectral VIs characterizing canopy greenness.
The field-measured set of vegetation parameters covers a wide range of values covering a large majority of natural variations in both LAI and MTA.This allowed us to use a simple comparison of correlation coefficients, which depend on the level of variation included in each of the studied variables.Since the measured MTA and LAI ranges are representative of their natural variation, a comparison of the values of R s tells us about the role of each vegetation parameter in determining the studied VIs.
According to both PROSAIL model simulations and empirical data, all studied VIs depended on MTA with NDVI being the least affected.This supports the findings of previous studies by Liu et al. [33].WDRVI is a nonlinear transformation of NDVI performed identically with NDVI when measured with Spearman's R, however, the more common Pearson's R 2 -value was slightly higher for WDRVI (0.25 and 0.30 for measured NDVI and WDRVI, respectively).NDVI and WDRVI were also the two indices with the strongest correlations with LAI in both the modeled and measured data (Table 2).However, the relatively strong (compared with other indices studied here) correlation between LAI and NDVI is compromised by the well-known saturation of the VI's sensitivity to LAI evident in the variation of the distance between the lines in Figure 1a.The same non-linearity is present in WDRVI (Figure 1f).EVI, EVI2 and MTVI2 were designed to be more sensitive to the variations inside the canopy by minimizing background influence, and working at large LAI values without saturation.In our measured dataset, the indices saturate less, with the lines in Figure 2b-e nearly parallel and more closely equidistant.Unfortunately, these VIs were more strongly correlated with MTA (R s < −0.71) than NDVI and WDRVI.The specific strength of the correlation between MTA and NDVI and WDRVI may somewhat depend on the spectral properties of soil.Using the sample curves in Figure 1b, a narrow-leafed lupin (MTA = 18 • ) canopy with LAI = 1 would have EVI ≈ 0.5.A similar EVI value is produced by an oat (MTA = 58 • ) canopy of LAI = 4.Our results demonstrate empirically that it is very difficult to quantify LAI using remotely sensed VIs when the MTA is unknown [49].
EVI and EVI2 are widely used for tracing the dynamics of vegetation gross primary production (GPP) and phenology inter-or inner-annually [50][51][52][53].As MTA had a stronger effect on these VIs than LAI, caution should be used when applying them to estimate LAI across structurally potentially different canopies, e.g., seasonality studies or application over multiple crop species.Depending on cover type and LAI values, it may be more optimal to use NDVI or WDRVI, which produced relatively high R s with LAI in both model simulations and empirical data: the R s between the LAI and NDVI was 0.71 and 0.58 in simulated and measured data, respectively.This also agrees quantitatively with the previous finding by Colombo et al. [54] of R 2 = 0.33 for data pooled over five vegetation types (soybean, corn, vineyards, poplar plantation and deciduous forest).
In our study, the correlation between Fcover and WDRVI was the strongest with the R s values close to the value of 0.96 reported by Gitelson [55] in the simulated dataset.They also proposed the following equation for calculating Fcover from remote sensing data: Fcover = 80.838WDRVI + 34.021 (6) (Equation ( 2) in [55]).The proposed equation was more accurate for empirical (RMSE = 0.22, bias of mean = 0.18) than model-simulated (RMSE = 0.38, bias of mean = 0.37) data with some overestimation in both datasets.
For measured data, the best-performing indexes for estimating Fcover were EVI2 and MTVI2, with R s = 0.88.In agreement with previous studies, we found that a new independent variable Fcover, a nonlinear combination of LAI and MTA, is the key driver of variation in the commonly used VIs.Hence, all the VIs studied here can be efficiently used to predict this structural parameter.The slightly better performance of EVI2 and MTVI2 can be attributed to the capability of the Spearman's R to ignore the nonlinearity with Fcover.When the more traditional Pearson's R is used (despite the non-normality of the data), the index most strongly correlated with Fcover is OSAVI (R = 0.85) with a marginal difference to EVI2 and MTVI2 (R = 0.83).Similarly to this study, however, some nonlinearity and dependence on illumination conditions is expected [56] under varying conditions.
The results corroborate earlier findings, especially those by Liu et al. [33] and Gitelson [55].As we have demonstrated empirically, even for homogeneous crop canopies, NDVI or its more sophisticated counterparts, EVI or EVI2, depend not only on LAI, but also on the second structural parameter, MTA.Furthermore, the dependence of the most widely used VIs on MTA can be stronger than on LAI with the actual effects of the two canopy parameters depending on their actual range of variation.Therefore, in order to map LAI from multispectral data, effective methods need to be developed to remotely measure leaf angles.This can be achieved, for example, using canopy reflectance in the red edge [22].Although Earth observation data in this spectral region is scarce, it is covered by some systems, most notably Sentinel 2. We hope that future investigations into the application of this sensor into MTA estimations will improve the mapping of MTA, and thus also LAI.

Conclusions
The results of this study reveal that the test six common multispectral vegetation indices showed strong correlations with MTA, especially for EVI, EVI2 and MTVI2, with R s < −0.75 for both PROSAIL simulations and field-measured data.All studied VIs were proxies for the fractional canopy cover determined by two key elementary canopy parameters (leaf area index and leaf mean tilt angle).For the studied crops, EVI2 and MTVI2 were the best indicator for Fcover with R s = 0.88 for measured data and R s = 0.93 for model simulations.

Figure 1 .
Figure 1.The soil reflectance used for PROSAIL simulations.

Figure 1 .
Figure 1.The soil reflectance used for PROSAIL simulations.
in the reference [20].Despite considerable variation, the strong effect of MTA on EVI, EVI2, OSAVI and MTVI2 is clear from Figure 3b-e.For the other indices (Figure 3a,f) R s was lower (R s = −0.53)and the index values at MTA = 54 • and MTA = 58 • (wheat and oat, respectively) could be as high as for MTA = 18 • and MTA = 27 • (lupin and faba bean, respectively).

Figure 3 .
Figure 3. Empirical relationships between MTA and VIs calculated from Spearman's correlation coefficient Rs: (a) NDVI, (b) EVI, (c) EVI2, (d) OSAVI, (e) MTVI2, (f) WDRVI.A simple trend line calculated from linear regression was added to visualize the trend of the variation.

Figure 3 .
Figure 3. Empirical relationships between MTA and VIs calculated from Spearman's correlation coefficient Rs: (a) NDVI, (b) EVI, (c) EVI2, (d) OSAVI, (e) MTVI2, (f) WDRVI.A simple trend line calculated from linear regression was added to visualize the trend of the variation.

Figure 4 .
Figure 4. Empirical relationships between LAI and VIs calculated from Spearman's correlation coefficient Rs: (a) NDVI, (b) EVI, (c) EVI2, (d) OSAVI, (e) MTVI2, (f) WDRVI.A simple trend line calculated from linear regression was added to visualize the trend of the variation.

Figure 4 .
Figure 4. Empirical relationships between LAI and VIs calculated from Spearman's correlation coefficient R s : (a) NDVI, (b) EVI, (c) EVI2, (d) OSAVI, (e) MTVI2, (f) WDRVI.A simple trend line calculated from linear regression was added to visualize the trend of the variation.

Figure 5 .
Figure 5. Empirical relationships between Fcover and VIs calculated from Spearman's correlation coefficient R s : (a) NDVI, (b) EVI, (c) EVI2, (d) OSAVI, (e) MTVI2, (f) WDRVI.A simple trend line calculated from linear regression was added to visualize the trend of the variation.

Table 1 .
Vegetation indices used in this study.

Table 1 .
Vegetation indices used in this study.

Table 2 .
Spearman's correlation coefficient (R s ) between spectral indices and canopy structural variables.