Exploring the use of vegetation indices to sense canopy nitrogen to phosphorous ratio in grasses

Reduced availability of plant nutrients such as nitrogen (N) and phosphorous (P) has detrimental effects on plant growth. Plant N:P ratio, calculated as the quotient of N and P concentrations, is an ecological indicator of relative N and P limitation. Remote sensing has already been widely used to detect plant traits in foliage, particularly canopy N and P concentrations and could be used to detect canopy N:P faster and at lower cost than traditional destructive methods. Despite the potential opportunity of applying remote sensing techniques to detect canopy N:P, studies investigating canopy N:P remote detection are scarce. In this study, we examined if vegetation indices developed for canopy N or P detection can also be used for canopy N:P detection. Using in situ spec- trometry, we measured the reflectance of a common grass species, Yorkshire fog ( Holcus lanatus L.), grown under different nutrient ratios and levels. We calculated 60 VIs found in literature and compared them to optimized VIs developed specifically for this study. The VIs were calculated using both the original narrow band spectra and the spectra resampled to the band properties of six satellite sensors (MSI – Sentinel 2, OLCI – Sentinel 3, MODIS – Terra/Aqua, OLI – Landsat 8, WorldView 4 and RapidEye) to investigate the influence of bandwidths and band positions. The results showed that canopy N:P was significantly related to both existing VIs (r 2 =0.16 - 0.48) and optimized VIs (r 2 =0.59 – 0.72) with correlations similar to what was observed for canopy N or canopy P. Existing VIs calculated with MSI and OLI sensors bands showed higher correlation with canopy N:P compared to the other sensors while the correlation with optimized VIs was not affected by the differences in sensors’ bands. This study might lead to future practical applications using in situ reflectance measurements to sense canopy N:P in grasslands.


Introduction
Nutrients play an essential role in plant growth and foliage nutrient concentration is linked to several physiological and ecosystem processes. Chlorophyll content, photosynthetic capacity, leaf life span, light use efficiency and biomass primary productivity have all been related to foliar nitrogen (N) concentration at the leaf or canopy level (Bakker et al., 2011;Evans, 1989;Green et al., 2003;Kergoat et al., 2008;Reich, 2012;Reich et al., 1999;Wright et al., 2004). Similarly, foliage phosphorus (P) has been correlated to leaf life span and photosynthetic capacity (Wright et al., 2004) and influences the photosynthetic rate (Vcmax) (Walker et al., 2014). Non-optimal levels of foliar N and P will often result in reduced plant growth (Zhao and Zeng, 2009). vegetation level, the N:P ratio is an informative variable that not only indicates potential limitation of N and P but which is also related to species composition, species richness, productivity and functional trait composition (Fujita et al., 2014;Roeling et al., 2018;Wassen et al., 2005).
Foliar N and P concentration can be measured in the laboratory, but this is labor intensive and expensive, especially for large sample sizes. An alternative approach consists of using remote sensing to estimate plant traits, among which foliage nutrient concentration (Homolová et al., 2013). Foliage N concentration influences the reflectance spectra through specific absorption features attributed to N-bonds and protein absorption. These absorption features are located at 1020 nm, 1510 nm, 1940 nm, 2060 nm, 2180 nm, 2300 nm and 2350 nm (Kumar et al., 2006). However, due to the overlapping absorption features of multiple compounds as well as the strong absorption by water in the shortwaveinfrared region (SWIR, 1400 -3000 nm), the interpretation of the reflectance spectra is difficult (Kumar et al., 2006). For this reason, foliage N detection by remote sensing is mainly linked to chlorophyll absorption and often includes regions of the spectrum associated with chlorophyll detection, i.e. the red-edge and near-infrared (NIR) regions (Schlemmer et al., 2013). Field spectrometry has been extensively applied to estimate canopy N in a variety of crops using so called spectral vegetation indices (VIs), which consist of a combination of spectral reflectance bands (Hansen and Schjoerring, 2003;Li et al., 2014;Schlemmer et al., 2013;Tian et al., 2011). Other ecosystems have been investigated for canopy N estimation using VIs computed from airborne or spaceborne sensors at different spatial resolutions, including temperate forests (HySpex airborne sensor at 3 m, Wang et al. (2016b)), tropical forests (RapidEye satellite sensor at 5 m, Cho et al. (2013)), Mediterranean forests (MERIS satellite sensor at 1 km, Loozen et al. (2018)), as well as chaparral vegetation (AVIRIS airborne sensor at 18 m, Serrano et al. (2002)) and savannah (RapidEye satellite sensor at 5 m, Ramoelo et al. (2012)).
Similar to canopy N, canopy P has been estimated using spectral indices, though to a lesser extent. Several studies have aimed to develop VIs for canopy P estimation in agricultural lands (Kawamura et al., 2011;Mahajan et al., 2014;Pimstein et al., 2011) and in a Carex dominated grassland (Wang et al., 2016a), all using field spectrometry.
VIs calculated from spectrometry measurements have already been extensively studied for canopy N detection (Pacheco-Labrador et al., 2014) and high accuracy results have been obtained (Tian et al., 2011). Different categories of VIs exist, the most common ones being the simple ratio (SR), the normalized difference (ND) and the simple difference (SD) (Tian et al., 2011;Wang et al., 2016a). A VI can be classified as either narrowband or broadband (Thenkabail et al., 2002). Narrowband VIs are computed from narrow reflectance bands measured through imaging spectrometry, which includes in situ reflectance measurement using a spectroradiometer as well as reflectance measured by specific airborne or spaceborne sensors (Tian et al., 2011;Wang et al., 2015). Broadband VIs are traditionally obtained from multispectral satellite sensors although they can also be computed from resampled narrow reflectance bands obtained from other sources (Clevers and Gitelson, 2013). A common approach to develop VIs is by selecting the optimal indices from all pair-wise combinations of wavelengths in the visible, NIR and shortwave-infrared (SWIR) regions (400-2500 nm), i.e. a band combination analysis (Hansen and Schjoerring, 2003;Tian et al., 2011).
Canopy N:P is relevant to ecological studies as it is a nutrient limitation indicator. It is thus somewhat surprising that remote sensing of canopy N:P using VIs has rarely been investigated. One study investigated canopy N:P detection in savanna grasslands using field spectrometry coupled with partial least square regression with significant results (r 2 = 0.69 -0.85, Ramoelo et al., 2013). Canopy N:P detection was also investigated in a boreal forest using vegetation indices calculated from airborne (Compact Airborne Spectrographic Imager, CASI) and spaceborne (Hyperion EO-1) imaging spectrometry (Gökkaya et al., 2015). The results showed that the VIs could be related to canopy N:P with r 2 ranging from 0.34 to 0.70. In this context, there is a need to examine in more detail how spectral properties, either from in situ measurement or satellite sensors, influence the performance of VIs regarding canopy N:P detection.
In this study, we aim at identifying if existing VIs that have been used for the estimation of canopy N or canopy P can also be used to remotely sense canopy N:P. To do so, we will evaluate three sub questions. How do VIs perform for canopy N:P estimation compared to canopy N and canopy P estimation? How do existing VIs perform compared to optimized VIs developed using the data collected in this study? How is the correlation between canopy N:P and the VIs influenced by different sensors' bandwidths and band positions?
To create a data set representing a large range of canopy N:P values, we chose to execute the study under controlled conditions. We have grown a common temperate grassland species, Holcus lanatus L. using six nutrient treatments that reflect N limited, P limited, and N and P colimited conditions, under high and low nutrient availability, respectively. Existing VIs found in the literature were computed using the original narrow band reflectance spectra measured by in situ spectrometry as well as the resampled spectra corresponding to six different satellites sensors (MSI aboard Sentinel 2, OLCI aboard Sentinel 3, MODIS aboard Terra-Aqua, OLI aboard Landsat 8, WorldView 4 and RapidEye). We compared the existing VIs to optimized VIs, obtained using a band combination analysis carried out on both the original narrow band spectra and the resampled broadband spectra.

Culture of the plants
Holcus lanatus L. (Yorkshire fog) is a perennial common grass species, generally found on nutrient-rich soils throughout Europe. Seeds were collected from the Middenduin nature reserve (52˚24΄N 4˚35΄E) located in the western Netherlands and cultivated under controlled conditions in a greenhouse of Utrecht University. Seeds were sown on moist quartz sand and were transplanted after 22 days into pots containing a mixture of quartz sand (Carlo Bernasconi, Zürich, CH, 0.1-0.7 mm) and sand collected close to the area where seed were harvested under a ratio of 11:1. This was done to introduce soil fauna, bacteria and fungi for a complete soil community. Each pot was planted with four seedlings. Plants were grown further under nutrient treatments that lasted from July 2015 until June 2016 in the greenhouse with 400 Wm −2 light from 9:00 am to 4:00 pm and a temperature regime that mimicked the temperate conditions of the Netherlands. Six nutrient treatments were applied with three N:P ratios, each in high and low nutrient supply levels. The nutrient ratios were N:P = 5 (N-limitation), N:P = 15 (balanced supply) and N:P = 45 (P-limitation). Each nutrient treatment had 8 replicates (n = 48). The nutrient treatments were applied following Güsewell (2005). N was supplied as KNO 3 and Ca(NO 3 ) 2 , P was provided as KH 2 PO 4 . KNO 3 and KH 2 PO 4 supplied part of the potassium, the rest was added as KCl. All the other essential macro-and micro-nutrients were provided in non-limiting supply using a standard Hoagland solution.

Reflectance measurements
The canopy reflectance of the grasses was measured inside the greenhouse under controlled conditions when plants were fully grown on February 19 th 2016 using a FieldSpec Pro Fr spectroradiometer (Analytical Spectral Device, Boulder, GO, USA). This device measures reflectance between 350 and 2500 nm with a resampled spectral resolution of 1 nm. We performed the spectral measurements with several lamps positioned above the plants on each side of the spectroradiometer to ensure constant light conditions. The spectroradiometer was held by a tripod at nadir position approximately 20 cm above the canopy. The field of view was 8°, ensuring a ground field of view of approximately 12 cm 2 . Each measurement of the spectroradiometer was the average of 50 successive scans. We measured each pot four times, after turning the pot by 90 degrees between each measurement to reduce geometrical effects. We averaged these four measurements to produce one canopy measurement per pot. The spectroradiometer was internally calibrated using a Spectralon white reference panel (Labsphere, North Sutton, NH, USA).

Leaf chemical measurements
Two leaves were sampled from each pot (n = 48) on February 25 th 2016, collected in paper bags and dried in an oven at 60°C for 48 h. The leaf samples were then ground using mixer mill (MM400, Retsch). The N concentration (g N 100 g −1 dry matter, %N) was measured using an elemental CN analyzer (Fisons NA 1500 NCS). The P concentration (g P 100 g −1 dry matter, %P) was measured using the total reflection X-ray fluorescence spectroscopy method (TXRF, S2 Picofox, Bruker, Germany). N:P ratios were calculated as the ratio between the weight based concentrations of N and P.

Spectral bands considered: original and resampled to satellite bands
Prior to analysis, the spectral range was reduced from 350-2500 nm to 400-2450 nm because the signal to noise ratio was too low for the minimum and maximum spectral ranges. This spectral range matches with the spectral coverage of earth observation sensors. We considered two sets of spectral bands: the original narrow band reflectance spectra and the spectra resampled to the properties of satellite sensor bands. We included the spectra resampled to satellite sensor bands to assess the possibility of detecting canopy N:P with VIs derived from actual satellite measurements. We included six different satellite sensors (MSI-Sentinel-2, OLCI-Sentinel-3, MODIS-Terra/Aqua, OLI-Landsat 8, WorldView 4, and RapidEye). We chose these sensors because they present a variety of band properties, from 5 to 21 bands and from 2.5 nm to 350 nm spectral resolution (Table 1, Fig. 1). They are all operational and have a long measurement record in the case of MODIS and Landsat. We obtained the resampled spectra by using the spectral response function corresponding to each sensor band. The spectral response functions were downloaded from the website of the sensor producer. For the MODIS sensor, the measured spectral response function could not be found and was approximated using a standard normal distribution (μ = 0, σ = 1).

Descriptive statistics and vegetation indices
Descriptive statistics and boxplots of canopy N:P, N and P, hereafter designated as canopy traits, were produced and the Pearson correlation coefficients between each of these canopy traits were calculated.

Existing VIs
We evaluated existing VIs, which we selected through a literature search. We found 60 VIs to evaluate (Table A1), including VIs that were either developed for N or P detection or, when developed for another purpose, for instance VIs correlated with chlorophyll, photosynthesis or the presence of vegetation (see the "developed for" column of the Table  A1 for a detailed list), have been used for N or P detection.
Based on their equations, we categorized the VIs as either two bands simple ratio (SR, band band 1 2 ), two bands normalized difference (ND, ) or three bands VI (TB). When they did not belong to any of the aforementioned categories, they were labelled other VI (OVI, see Table A1 for equations associated with TB or OVI).
We calculated linear models between canopy N:P, canopy N or canopy P and both the VIs and the natural log transformed VIs calculated using the original narrow band spectra.
The ten narrow band VIs that obtained the best correlation (r 2 ) with canopy N:P were also calculated using the spectra resampled to the six satellite sensor bands mentioned in section 2.4. The indices were computed using the sensor bands closest to the VI nominal wavelengths. Several VIs could not be computed because the sensors' band locations were too distant (> 100 nm) from the VI nominal wavelength.

Optimized VIs
From the spectral signature we took from our experimental data, we calculated optimal VIs using combinations of all available wavelengths, i.e. a band combination analysis. Following Tian et al. (2011), every combination of two, 1 and 2 , or three wavelengths, 1 , 2 and 3 , was used to compute two bands or three bands VIs. The two bands VIs were categorized as simple ratio (SR) 1 2 , single difference (SD) . In order to decrease computation time, the three bands VIs were computed using only 1 out of 10 narrow bands (Pacheco-Labrador et al., 2014). Subsequently, linear regressions between each of the obtained VIs and canopy N:P, canopy N and canopy P were calculated. Specific wavelength regions of high correlation between the canopy traits and the three categories of two bands VI were identified using heatmaps. The optimized VIs were obtained by selecting the band combination with the highest r 2 , for both the narrow band spectra and the spectra resampled to satellite sensor bands.

Regression models
All regression models were assessed using the determination coefficient (r 2 ) and the cross validated Relative Root Mean Square Error (RRMSE cv ) values, obtained using the leave-one-out cross validation method (Clevers and Gitelson, 2013) and calculated following Eq. (1) (Yao et al., 2010): , is a measurement, with n the total number of measurements, P i represents the predicted value, O i , the observed value, and Ō i the mean of all observed values. The significance level was at a pvalue < 0.01. All the statistical analyses were performed in the R environment (R Development Core Team, 2014).

Descriptive statistics of canopy N:P, canopy N and canopy P
The results of the descriptive statistics showed that the range of canopy N:P values (6.1-75.0, Fig. 2) corresponds to the range usually observed in natural environments (Roeling et al., 2018). As expected, canopy N:P was correlated with both canopy N and canopy P (Table 2).

Existing vegetation indices
Among the 60 existing VIs tested, 31 VIs showed significant relationships with at least one canopy trait, while 29 VIs showed no significant relationship with any of the canopy traits (Table A1). The VI showing the highest correlation with canopy N:P (r 2 = 0.48), was the TB (R 498 , R 413 , R 442 , Fig. 3) VI developed for N detection (Tian et al., 2011). The results obtained for canopy N and canopy P were similar, with the highest r 2 equal to 0.44 for canopy N with ND (R 1220 , R 710 ) VI. For canopy P, the highest correlation (r 2 = 0.52) was obtained with TB (R 498 , R 413 , R 442 ), hence the same VI as for canopy N:P. The two VIs showing the highest correlation with canopy N:P, TB (R 498 , R 413 , R 442 ) and TB (R 434 , R 496 , R 401 ) were both based on three wavelengths located in the blue region of the spectrum. Other VIs that showed a significant   (r2) between the canopy trait and each combination of two bands (λ1 and λ2, nm) between 400-2450 nm for canopy N:P, canopy N and canopy P and each VI category investigated. SR, Simple Ratio; SD, Simple Difference, ND, Normalized Difference. relationship with canopy N:P were based on the NIR and red-edge regions of the spectrum.

Optimized vegetation indices
We used heatmaps to investigate the wavelength regions of high correlation between the two bands VIs, obtained from the band combination analysis, and the three canopy traits (Fig. 4). A first region of high correlation with canopy N:P (r 2 between 0.4 and 0.7) was located in the blue region for 1 (400-500 nm) and in the green -red region for and also had high correlations for canopy N and canopy P.
The optimized VI TB3 (R 541 , R 421 , R 531 , Fig. 5), located in the blue and green region of the spectrum showed the highest correlation with canopy N:P and had a higher r 2 (r 2 = 0.72, Table A2) than the best performing optimized VIs for canopy N (r 2 = 0.69; TB3, R 431 , R 2431, R 561 ) and for canopy P (r 2 = 0.67; TB1, R 531 , R 541, R 421 ).
Regarding the influence of the VI category on the result, the optimized VIs from the SD category, which were composed of a combination of bands from the red-edge and SWIR regions, showed lower correlations (r 2 = 0.53 -0.59) than all the other VI categories (r 2 = 0.64 -0.72), mostly composed of a combination of blue and green bands. This result holds for all canopy traits (canopy N:P, N and P).
Compared to the results obtained for the existing VIs, the optimized VIs performed better regarding r 2 values (Fig. 6) for canopy N:P. The best optimized VI explained 24 percent more of the variation in canopy N:P values compared to the best performing existing VI.

Existing vegetation indices
The TB (R 434 , R 496 , R 401 ) VI calculated for the MSI sensor aboard Sentinel 2 showed the highest correlation with canopy N:P (r 2 = 0.63), canopy N (r 2 = 0.58) and canopy P (r 2 = 0.62) compared to the other existing VIs tested (Table A3). Regarding the influence of the sensor band properties on the result, the existing VIs calculated using the MSI and OLI sensors band showed in general higher correlation with canopy N:P than the same VIs calculated based on OLCI or MODIS sensors (Fig. 7). For the existing VIs TB (R 498 , R 413 , R 442 ), TB (R 434 , R 496 , R 401 ) and TB (R 705 , R 717 , R 491 ) calculated based on MSI sensor bands, the TB (R 498 , R 413 , R 442 ) and TB (R 434 , R 496 , R 401 ) VIs calculated based on OLI sensor bands as well as for the TB (R 705 , R 717 , R 491 ) VI calculated based on OLCI sensor properties, the obtained r 2 were higher than for the same VIs calculated with the original narrow band spectra. Similar results were observed for canopy N and canopy P. On the contrary, all the existing VIs calculated with the spectra resampled to MODIS bands showed lower r 2 with canopy N:P compared to their narrow band  spectra counterpart and three existing VIs showed an non-significant relationship with canopy N:P.

Optimized vegetation indices
Canopy N:P could be related to the optimized VIs calculated from the resampled spectra with r 2 values ranging from 0.48 to 0.72 (Table  A2). Similar results were obtained for canopy N and canopy P. The SR, ND and TB categories of optimized VIs mostly performed better with canopy N:P (r 2 between 0.55 -0.72) compared to the SD category of VI (r 2 between 0.48 -0.59). Among the bands selected for the SR, ND and TB categories of optimized VIs, the blue and green regions of the spectrum were dominant, similar to what was observed for narrow bands optimized VIs.
The difference in satellite sensors band properties did not greatly affect the results as the r 2 values were in the same range for all sensors per VI category. The satellite sensors WorldView 4 and RapidEye, both with five bands and bandwidths between 35-140 nm, showed r 2 between 0.55 and 0.67 with canopy N:P for the SR, ND and TB categories of VI. This is comparable to the results obtained for satellites sensors with a higher number of bands and narrower bands, e.g. MSI and MODIS, for which the obtained r 2 values were between 0.67 and 0.72 for the same VI categories.
When comparing the performance of existing and optimized VIs for canopy N:P detection (Fig. 6), optimized VIs generally performed better than existing VIs for all satellite sensors tested. However, the existing VIs calculated for the MSI and OLI sensors showed r 2 values in an overlapping range with their optimized counterparts.

Existing vegetation indices
Existing VIs performed similarly for canopy N:P detection as they did for canopy N and canopy P regarding obtained r 2 values. More specifically, the two best performing VIs for canopy N:P in this study, TB (R 498 , R 413 , R 442 ) and TB (R 434 , R 496 , R 401 ), both based on the blue region of the spectrum, proved to be interesting candidates for canopy N:P detection. In previous studies, these VIs were highly correlated to canopy N in rice crops (r 2 = 0.81 and r 2 = 0.84, respectively (Tian et al., 2011), but showed a weak correlation with canopy N in Holm oak leaves (Pacheco-Labrador et al., 2014).
Among the 60 VIs tested, almost 50% (n = 29) failed to show a significant relationship with canopy N:P, canopy N or canopy P while previous studies found a significant relationship. This unreliability in prediction accuracy for previously validated VIs was already observed by Pacheco-Labrador et al. (2014), who found that the majority of the published vegetation indices tested could not be correlated to canopy N in Holm oak leaves. This might be explained by differences in growth conditions and species investigated as well as by the influence of the range in canopy N or canopy P values. This highlights the influence of the dataset on the prediction accuracy.

Optimized vegetation indices
The optimized VIs developed specifically for this study performed better for canopy N:P detection compared to the existing VIs tested, for both r 2 and RRMSE cv values.
The SR, ND and TB categories of optimized VIs were based on a combination of wavelengths from the blue and green regions of the spectrum. The blue region is linked to pigment absorption, which peaks at 430 nm for chlorophyll a (chlor-a) and 460 nm for chlorophyll b (chlor-b) (Kumar et al., 2006). The SD, ND and TB 3 optimized VIs for canopy N:P included wavelengths located at 427 nm and 421 nm, hence close to one of the absorption peaks of chlor-a. Other indices located in the blue and green regions of the spectrum have been related to chlor-a and chlor-b concentration in Vitis vinifera leaves (the BGI2 index, SR, R 450 , R 550 ; Zarco-Tejada et al. (2005)) or to canopy N in rice (ND, R 573 , R 444 ; Tian et al. (2011)). The combination of wavelengths from the blue and green regions might thus be linked with chlorophyll and has also been explored for canopy N detection (Tian et al., 2011). Our results Fig. 7. Comparison between the r 2 obtained for the relationships between canopy N:P and the existing VIs calculated using either the original narrow band spectra or the spectra resample to satellite sensors'band properties (MSI, OLCI, MODIS, OLI). No existing VI was calculated for the RapidEye and WorldView 4 sensors because their band positions are too distant (> 100 nm) from the nominal VIs' wavelengths.
showed that this wavelength combination is correlated with canopy N:P, canopy N and canopy P.
The optimized SD VI for canopy N:P combined wavelengths located in the red-edge (718 nm) and in the SWIR regions (1577 nm). The SWIR region is characterized by several nitrogen and protein absorption features. The absorption features located at 1500 nm and 1510 nm might explain the second wavelength (1577 nm) being selected, because the absorption features are known to be broadened due to scattering (Kumar et al., 2006). However, the SWIR region of the spectra is influenced by the absorption features of many compounds that interfere with each other (Kumar et al., 2006) which render the signal difficult to interpret. In particular, the strong absorption due to water molecules greatly influences the reflectance. This might be the reason why the proportion of variance explained by the SD category of optimized VIs (r 2 = 0.53 -0.59) was lower than for the SR, ND and TB categories of optimized VIs (r 2 = 0.64 -0.72).

Spectra resampled to satellite sensors' bands 4.2.1. Existing vegetation indices
The correlation between canopy N:P and the existing VIs was influenced by the sensors' band properties because different sensors showed different correlation with canopy N:P for the same VIs. Two VIs (TB R 498 , R 413 , R 442 and TB R 434 , R 496 , R 401 ) calculated from MSI and OLI sensors showed higher correlation (r 2 = 0.55 -0.63) than when calculated from the OLCI and MODIS sensors (r 2 = 0.36 -0.40) or even narrow band sensors (r 2 = 0.44 -0.48) (Fig. 7). The bands from MSI and OLI sensors used for the VIs calculation were broader than their counterparts from the OLCI and MODIS sensors and hence from the original narrow band spectra. The bandwidth of the sensors thus influenced the correlation obtained for some of the existing VIs.

Optimized vegetation indices
The correlation between canopy N:P and the optimized VIs was not greatly influenced by the sensors' band properties and the correlation was stable across the different satellite sensors tested. Although this is in contradiction to what was observed for existing VIs, this might indicate that broad band sensors, such as RapidEye and WorldView 4, could be useful for canopy N:P detection. However, our results consistently showed that the blue region was important for canopy N:P detection. Although VIs based on the blue region of the spectrum do not represent a difficulty for in situ studies, including the blue region for VIs calculated from satellite sensors might lead to an interpretation problem due to atmosphere Rayleygh scattering.

Future perspectives
This analysis investigated the possibility of detecting canopy N:P using VIs compared to canopy N and canopy P detection as well as the influence of the sensors band properties on the relationships. The results obtained in this analysis might have concrete in-situ application perspectives in the field of ecology given the importance of the canopy N:P for biodiversity studies (Fujita et al., 2014;Roeling et al., 2018;Wassen et al., 2005). Spectral VIs might be a useful method to detect canopy N:P in grasslands in a non-destructive and time efficient manner that would, for example, allow to monitor the seasonal evolution of canopy N:P or trends in changing N:P ratios in response to eutrophication or other global change factors (Wassen et al., 2013). However, as with remote sensing of canopy N or canopy P, remote sensing of canopy N:P is not a direct measurement of canopy N:P values. Especially, the correlation between canopy N:P and both canopy N and canopy P (Table 2) renders it difficult to distinguish between the separate influences on the reflectance signal. Moreover, similar studies investigating canopy N:P detection should be done on more plant species and plant communities in order to validate the results obtained in this study. Further studies should also investigate the influence of the spatial resolution of the satellite sensors on canopy N:P detection as this needs to be studied before actual satellite sensors measurements can be exploited to sense canopy N:P in natural environment.

Conclusion
Canopy N remote detection has already been extensively studied using both spectroradiometer and satellite measurements. On the contrary, canopy N:P, despite being an important indicator of nutrient limitation, has seldom been studied with remote sensing techniques. In this study, we investigated the possibility of detecting canopy N:P in the common grass species Holcus lanatus using VIs developed for canopy N and canopy P detection. The results showed that using VIs for canopy N:P detection was as effective as when applied for canopy N or canopy P detection. This held for both existing and optimized VIs as well as for narrow band and broader band VIs calculated from the spectra resampled to the spectral properties of six different satellite sensors. The influence of different satellite sensors' band properties was unclear as it differed between existing and optimized VIs. Existing VIs calculated with MSI and OLI sensors bands showed higher correlation with canopy N:P compared to the other sensors tested. In contrast, the correlation with optimized VIs was not affected by the differences in sensors' bands. Satellite sensors with a limited number of broad bands, such as WorldView 4 and RapidEye, yielded similar results as sensors with multiple and narrower bands, like MSI or Sentinel 3. In the future, these results might lead to practical applications using handheld spectrometers for in situ canopy N:P detection in grasslands. The observed consistent importance of the blue region of the spectrum for canopy N:P detection might render canopy N:P detection with actual satellite sensors complicated due to the interference with the Rayleygh scattering.

Competing interest
The authors declare that they have no competing interests.

Acknowledgments
This work was supported by The Netherlands Organization for Scientific Research (NWO) [NWO ALW-GO-AO/14-12]. We would like to acknowledge Ton Markus for his help with the figures presented in this paper.

Table A1
Relationship between the existing vegetation indices (VIs) (n = 60) and canopy N:P, canopy N and canopy P. Coefficient of determination (r 2 ), p-value and Relative Root Mean Square Error of cross-validation (RRMSEcv).