Assessment of Red-Edge Position Extraction Techniques: A Case Study for Norway Spruce Forests Using HyMap and Simulated Sentinel-2 Data

Systematic quantification and monitoring of forest biophysical and biochemical variables is required to assess the response of ecosystems to climate change and gain a deeper understanding of the carbon cycle. Red-Edge Position (REP) is a hyperspectrally detectable parameter, which is sensitive to Chlorophyll (Chl) content. In the current study, REP was modelled for Norway spruce Forest canopy Reflectance and Transmittance (FRT) using Radiative Transfer Modelling (RTM) (resampled to HyMap and Sentinel-2 spectral resolution) as well as calculated from the real HyMap and simulated Sentinel-2 image data. Different REP extraction methods (PF, LE, 4PLI and its optimized versions for HyMap and Sentinel-2 spectral resolution) were assessed. The lowest differences in REP values calculated from image-extracted spectra and from the theoretical RTM simulations were found for the 4PLI method including its HyMap and Sentinel-2 optimized versions (4PLIH and 4PLIS). Despite its simplicity, the 4PLI REP extraction technique demonstrated its potential usefulness for estimating canopy chlorophyll (Chl × LAI) content using both airborne hyperspectral (HyMap) data as well as space-borne Sentinel-2 image data.


Introduction
Leaf pigments have been identified as important indicators of physiological status, senescence and stress [1][2][3], while Chlorophyll (Chl) is the leaf pigment responsible for photosynthesis.Estimates of foliar biochemicals such as the level of Chl provide us with information on the plant's development stage, productivity and stress [2,4,5]; therefore, Chl is one of the key biophysical variables used for vegetation monitoring.Direct field methods using standard approaches for estimating Chl content frequently require destructive harvesting and they are labour-intensive and costly.
The biochemical characteristics of vegetation, including Chl, affect its spectral properties and remote sensing, especially hyperpspectral remote sensing, offers an alternative set of techniques in different domains, such as laboratory measurements, field and imaging spectroscopy; therefore, it has opened new ways of monitoring Chl content over a wide range of scales [6,7].
The relationship between the spectral properties of vegetation and the biochemical and biophysical properties can be studied using empirical correlations.Such approaches are relatively straightforward but unfortunately lack robustness as the empirical correlations of vegetation's spectral properties and biophysical/biochemical parameters are site and time specific.To overcome this, physical models have been developed.Physical methods employ inverted Radiative Transfer (RT) models (e.g., PROSPECT [8,9], LEAFMOD [10], SAIL [11] and SCOPE [12]) to estimate the biochemical content at the leaf level [13].RT modeling simulates the transfer of radiation in the canopy by computing the interaction between a plant and solar radiation [14,15].In comparison with the direct extrapolation and the canopy-integrated approaches, inversion models offer the potential of a more generic approach to quantify the biochemical parameters of vegetation based on spectral data.However, numerous in-put parameters are requested for the RT models, which are still difficult and costly to obtain for regional studies.
The wavelength of the maximum slope in the region of the red-edge (680 to 780 nm), termed the Red-Edge Position (REP) [16,17], is a spectral measure that is less sensitive to the effect of variable leaf and canopy Chl content and environmental conditions on Chl content estimation.REP is an important index that indicates the abrupt leaf reflectance change between 680 nm and 780 nm of the vegetation spectra caused by the combined effects of strong Chl absorption and leaf internal scattering [18].Increases in the amount of Chl result in a broadening of the major Chl absorption feature centred around 680 nm [18,19], causing a shift in the slope and REP towards longer wavelengths [16,[20][21][22].Therefore, REP, due to it being sensitive to Chl variation, provides a valuable remote sensing means for estimating leaf or canopy Chl content [20,23,24].Various techniques have been developed for the extraction of REP parameters from different sources of spectral data, among the most widely used REP extraction methods are such methods as the Maximum First Derivative (MFD) [22, [25][26][27][28], the Polynomial Fitting (PF) technique [29], the Inverted Gaussian (IG) technique [30] and Linear Extrapolation (LE) technique [7].
While laboratory, field and imaging spectroscopy have been used to extract REP and estimate Chl content, new opportunities are approaching since the new state-of-the-art satellite, Sentinel-2A, was launched on 23 June 2015.Sentinel-2A (S2), a superspectral instrument, is the first satellite of the S2 constellation carrying 13 bands, covering the Visible-Near Infrared (VNIR) and Shortwave-Infrared (SWIR) domains with a spatial resolution ranging from 10 to 60 m.Effective utilization of S2 data, which is now in orbit, provides a very good opportunity for vegetation studies on a regional and continental level especially due to the high revisit time of S2 satellites [31][32][33][34].In addition, S2 data contain more spectral bands within the red-edge than data provided by e.g., Landsat, ASTER and other multispectral space-borne instruments.This allows such techniques as REP calculation, which were originally developed for hyperspectral data, to be employed for S2 data.This is specially related to the bands centred at 665 (band 4), 705 (band 5), 740 (band 6) and 783 (band 7) nm [35].
The potential of S2 data has been investigated for diverse vegetation parameter estimations such as LAI [36][37][38], Chl and nitrogen [39], biophysical parameters [31,40] and REP [39,41].However, only a limited number of forest-oriented applications have been published so far.This may be due to the fact that forests are generally more challenging targets than agricultural fields, as it is also much more time and cost-demanding to collect the ground-truth data.In addition, the architectural structure, clumping of optically active surfaces at multiple scales and spatial-temporal foliage dynamics contribute to this problem.Assessment of forest biophysical and biochemical variables faces further complexity by the contribution of vegetation, litter, soil, bark as well as plant and relief shadow, all of which influence the radiometric signal [42].As a result, all these factors together make forest-oriented remote sensing applications more challenging than compared to agriculture ones.Only three studies have been published so far on this topic: one deals with discrimination of tropical forests types using hyperspectral and simulated S2 data [33].In the second one, the potential of EnMap and S2 data was assessed for detecting drought stress phenomena in deciduous forests [32].The third study, which is the only one so far that utilizes the real S2 data rather than simulated datasets, deals with the suitability of the S2 red-edge wavelengths for burn severity discrimination [43].
The main goal of this paper was to test REP extraction methods while focusing on those suitable for both the broader-band hyperspectral data used in this study (e.g., HyMap) as well as for space-borne S2 data.This study was forest-oriented and intends to answer the following questions: (i) What are Forests 2016, 7, 226 3 of 17 the differences in the REP values calculated from the theoretical Radiative Transfer Model (RTM) simulations and the real (hyperspectral-image-extracted) spectra, and which of the studied REP calculation methods provides the best fit between the theoretical RTM simulations and the real hyperspectral image data (HyMap)?(ii) What is the influence of spectral sampling on the REP calculation reliability (comparison of theoretical RTM simulations, real hyperspectral image data (HyMap) and simulated S2 data)?(iii) What do the relationships between the calculated REP values and canopy Chl (Chl × LAI) content look like when using the theoretical RTM simulations, real hyperspectral image data (HyMap) and simulated S2 data?

Study Area
The study was performed in the Sokolov lignite basin situated in the western part of the Czech Republic.This region has been affected by long-term lignite mining and heavy industrial pollution.The central part of the basin is surrounded by boreal forests with Norway spruce (Picea abies L. Karst.) as the predominant species, which is the main object of interest in this study.The current research was conducted for four homogeneous (i.e., single-species and evenly aged) Norway spruce forests stands: Erika (E), Habartov (H), Mezihorská (M) and Studenec (S) located in the neighbourhood of the lignite mining area.The trees did not show any visible symptoms of damage at the time of data acquisition [44], although these forests are situated in an area with a historically high level of environmental pollution.The location of the Sokolov lignite basin with the positions of the studied Norway spruce stands can be seen in Figure 1.simulations and the real (hyperspectral-image-extracted) spectra, and which of the studied REP calculation methods provides the best fit between the theoretical RTM simulations and the real hyperspectral image data (HyMap)?(ii) What is the influence of spectral sampling on the REP calculation reliability (comparison of theoretical RTM simulations, real hyperspectral image data (HyMap) and simulated S2 data)?(iii) What do the relationships between the calculated REP values and canopy Chl (Chl × LAI) content look like when using the theoretical RTM simulations, real hyperspectral image data (HyMap) and simulated S2 data?

Study Area
The study was performed in the Sokolov lignite basin situated in the western part of the Czech Republic.This region has been affected by long-term lignite mining and heavy industrial pollution.The central part of the basin is surrounded by boreal forests with Norway spruce (Picea abies L. Karst.) as the predominant species, which is the main object of interest in this study.The current research was conducted for four homogeneous (i.e., single-species and evenly aged) Norway spruce forests stands: Erika (E), Habartov (H), Mezihorská (M) and Studenec (S) located in the neighbourhood of the lignite mining area.The trees did not show any visible symptoms of damage at the time of data acquisition [44], although these forests are situated in an area with a historically high level of environmental pollution.The location of the Sokolov lignite basin with the positions of the studied Norway spruce stands can be seen in Figure 1.

Airborne Hyperspectral Data Acquisition
Two HyMap airborne hyperspectral datasets were acquired over the area of interest.On 27 July 2009, the instrument was flown at an average altitude of 2570 m and nine cloud-free flight lines in a NE-SW direction were obtained with a ground spatial resolution of 5 m.The second HyMap dataset was acquired on 21 August 2010, obtaining seven cloud-free flight lines at an average altitude of 2475

Airborne Hyperspectral Data Acquisition
Two HyMap airborne hyperspectral datasets were acquired over the area of interest.On 27 July 2009, the instrument was flown at an average altitude of 2570 m and nine cloud-free flight lines in a NE-SW direction were obtained with a ground spatial resolution of 5 m.The second HyMap dataset was acquired on 21 August 2010, obtaining seven cloud-free flight lines at an average altitude of 2475 m with 4 m ground spatial resolution.Both datasets covered the spectral range between ca.450-2500 nm with a FWHM of ca.10-15 nm.
The atmospheric correction was performed using the ATCOR-4 software with the aerosol optical thickness obtained from in-situ measurements by a Microtops II Sunphotometer.For validation purposes, the Hemispherical-Conical Reflectance Factors (HCRF) of several near-Lambertian surfaces (asphalt, concrete, artificial grass and a water pool) were measured.The data were acquired using an ASD FieldSpec 3 spectroradiometer, equipped with a standard pistol grip, simultaneously with the airborne data acquisition.Since the HyMap data showed strong Bidirectional Reflectance Distribution Function (BRDF) effects due to the wide field of view, semi-empirical nadir normalization was performed using the Ross-Li kernel approach [45,46].Ortho-georectification of the 2009 dataset was performed in PARGE software, while the 2010 data were corrected using the ORTHO module.More details about the HyMap 2009 and 2010 data pre-processing can be found in Kopačková et al. [44].

Simulated Sentinel-2 Image Data
Sentinel-2 data was simulated using the S2 end-to-end simulation tool (S2eteS), which is described in detail in Segl et al. [47].It consists of a forward simulation module and a backward simulation module.
The forward simulation consists of a spatial module, a spectral module, an atmospheric module and a spectral module.The spatial module creates a virtual S2 sensor overpass over the image scene, taking into account parameters such as the specific sensor observation geometry Segl et al. [47].The atmospheric module converts the surface reflectance data to Top Of Atmosphere (TOA) reflectance taking into account atmospheric parameters such as e.g., Aerosol Optical Thickness (AOT), variable Columnar Water Vapour (CWV), followed by a radiance calculation.The spectral module performs the spectral resampling of the resulting data.The resultant at-sensor radiance is converted to DN in the radiometric module considering e.g., sensor noise and detector nonlinearity Segl et al. [47].
The backward simulation then uses the S2 pre-processing chain to retrieve simulated-at-ground reflectance data, which include the aforementioned sensor specific parameters.More detailed information on the whole simulation process can be found in Segl et al. [47].
This allows a close-to-real sensor simulation of the S2 data, rather than just a simple spectral and spatial resampling of the HyMAP data to S2 data.Therefore, this approach allows close-to-real-world tests and demonstrations of S2 applications such as, for example, mapping of iron-bearing surface cover types [48] or the REP calculation prior to the launch of the sensor presented here.

In-Situ and Laboratory Measurement of Vegetation Properties
The needle sampling scheme was described in detail in Kopačková et al. [44].In this study, the Chl content determined for the Upper (U 1 ) portion of the sunlit crown of 50 mature Norway spruce trees at four study sites was used: Erika (10 trees), Habartov (15 trees), Mezihorská (15 trees) and Studenec (10 trees).The samples were collected from the same trees in both years.The total Chl was extracted in DiMethylForamide (DMF) based on Porra et al. [49].The foliar pigment content was then calculated using the equations of Wellburn [50].The samples were scanned by a table scanner to determine their projection area (LAP), which was then transformed to the total needle area (LAT) using the LAP/LAT empirical correction factors [51].LAT was then used to calculate the Specific Leaf Area (SLA) describing the ratio of needle total area and dry weight.The SLA values were finally used to normalize the values of Chl by sample area.Another part of the collected samples was used for the spectral measurement of the needles using ASD FieldSpec 3 [52].In addition, the spectra of Norway spruce bark and dry litter were measured as inputs for the canopy RTM.These measurements were performed using the SR-2500 spectroradiometer (Spectral Evolution Inc., Lawrence, MA, USA) combined with the contact probe.A set of Digital Hemispherical Photographs (DHP) was acquired at each site in order to estimate the canopy biophysical and structural parameters using a Canon EOS 400D camera (Canon Inc., Tokyo, Japan) equipped with a Sigma 4.5 mm Circular fisheye lens.The acquired DHPs were pre-processed in GIMP 2.8 software (GIMP Development Team).CanEye software (Canon Inc., Tokyo, Japan) was finally used to estimate the true and effective Leaf Area Index (LAI), CAnopy Closure (CAC) and CRown Cover (CRC).

Leaf and Canopy Radiative Transfer Modelling
The ASD reflectance acquired for the needle sample calibration dataset (2 × 30 samples) was used for the parameterization of the PROSPECT-5 leaf RTM [8,53].This model was introduced by Jacquemoud and Baret [8] and currently it is one of the most frequently used leaf-level RTM due to its simplicity and very good inversion ability [9], which were the main reasons for using this model.A validation dataset (2 × 20 samples) was used for an independent assessment of the needle simulations.The main outputs of PROSPECT-5 forward simulations are single-leaf Reflectance (R) and Transmittance (T) at particular wavelengths.PROSPECT-5 was run in the forward mode with known values of Chl on the input with a structural parameter N ranging from 1.0 to 5.0 (0.1 step).The simulated R and T were transformed to infinite reflectance (R ∞ ) using the equations of Zarco-Tejada and Miller [54] as the ASD outputs (when using a contact probe to measure needle reflectance) have more the character of R ∞ : where: R ∞ is infinite reflectance at the wavelengths lambda; R is single leaf reflectance at wavelength lambda; T is single leaf transmittance at wavelength lambda.Infinite reflectance (R ∞ ) is defined as the reflectance of a leaf stack where the influence of multiple radiation scattering is taken into account.The simulated R ∞ spectra were then compared with the ASD reflectance to find the optimal value of the parameter N using the Root Mean Square Error (RMSE) merit function minimization technique.The optimal N value was used together with the Chl values of the validation dataset samples to assess the accuracy of the PROSPECT-5 simulations.The simulated spectra were again transformed from R and T to R ∞ and then compared with the corresponding measured ASD spectra of the validation dataset used for calculating the RMSE of simulation.
The simulated R and T were up-scaled to the shoot level using the spectral invariants theory following Rautiainen et al. [55].This theory is based on an assumption that the canopy albedo at the given wavelength λ is related to leaf/needle albedo on the same wavelength λ trough, the photon recollission probability p [56,57]) and represents an alternative approach for modelling the radiation scattering occurring in a vegetation canopy.U 1 shoot-level spectra were simply averaged to obtain the single-shoot-level spectrum for each sampled tree.Mean shoot-level spectra were then determined for each particular stand by averaging the spectra of the corresponding sampled trees, which were finally used as the input of the canopy RTM.
The FRT model was used for up-scaling the needle/shoot-level spectra to the canopy [58].One single mock-up (a file definition of the properties of the vegetation stand whose spectral signature is going to be simulated by the model) was prepared for all the studied sites.First, Diameter Breast Height (DBH) and crown radius (R c ) were determined using a combination of in-situ and two allometric equations [59] as below: where: h is tree height; DBH is diameter breast height; R c is crown radius.Crown Cover (CRC) was calculated from estimated CAC by processing the DHPs and the given fixed value of the tree distribution parameter (TDP = 1.5) [59].Then, the tree density was estimated Forests 2016, 7, 226 6 of 17 using the CRC and R c values.The Dry Leaf Weight (DLW) was finally calculated from the Density (DEN) (determined in the previous step), Specific Leaf Weight (SLW) and LAI (estimated from the DHPs).The measured Norway spruce bark spectrum was used to define the spectral properties of the wooden parts of the trees considered in the FRT model and the measured dry litter spectrum was utilized as the forest floor understory background spectrum.The parameterized FRT model was run in the forward mode to provide a canopy spectrum for each considered forest stand.The simulated spectra were resampled from the original 5 nm sampling to the spectral resolution of the HyMap and S2 data.Finally, they were compared with the corresponding spectra extracted from the HyMap and simulated S2 image data.

Red-Edge Position Calculation
The REP, resampled to HyMap and S2 spectral resolution, was calculated from the FRT canopy RTM as well as from the real HyMap (5-m spatial resolution) and simulated S2 imagery (20-m spatial resolution).The REP values extracted from the FRT simulations are signed as REP SIM and the values calculated from the image datasets are signed as REP IMG in the following text.Each Norway spruce stand was determined by Regions Of Interest (ROI) (Figure 2).The spatial extend of the ROI was defined with respect to the S2 data spatial resolution (20 m) In order to obtain representative statistical inputs in terms of HyMap/S2 spatial resolutions and the spatial distributions of sampling points used for the Chl × LAI content estimation, the image pixels were grouped to ROIs.The same ROI was used in both HyMap (5-m spatial resolution) and simulated S2 imagery (20-m spatial resolution).The spectra of the pixels falling within the ROI were used in the following analyses (e.g., REP assessment).The averaged REP values derived for each ROI were then correlated with the per-ROI-averaged Chl × LAI content (see more in Section 2.8).
Forests 2016, 7, 226 6 of 17 (DEN) (determined in the previous step), Specific Leaf Weight (SLW) and LAI (estimated from the DHPs).The measured Norway spruce bark spectrum was used to define the spectral properties of the wooden parts of the trees considered in the FRT model and the measured dry litter spectrum was utilized as the forest floor understory background spectrum.The parameterized FRT model was run in the forward mode to provide a canopy spectrum for each considered forest stand.The simulated spectra were resampled from the original 5 nm sampling to the spectral resolution of the HyMap and S2 data.Finally, they were compared with the corresponding spectra extracted from the HyMap and simulated S2 image data.

Red-Edge Position Calculation
The REP, resampled to HyMap and S2 spectral resolution, was calculated from the FRT canopy RTM as well as from the real HyMap (5-m spatial resolution) and simulated S2 imagery (20-m spatial resolution).The REP values extracted from the FRT simulations are signed as REPSIM and the values calculated from the image datasets are signed as REPIMG in the following text.Each Norway spruce stand was determined by Regions Of Interest (ROI) (Figure 2).The spatial extend of the ROI was defined with respect to the S2 data spatial resolution (20 m) In order to obtain representative statistical inputs in terms of HyMap/S2 spatial resolutions and the spatial distributions of sampling points used for the Chl × LAI content estimation, the image pixels were grouped to ROIs.The same ROI was used in both HyMap (5-m spatial resolution) and simulated S2 imagery (20-m spatial resolution).The spectra of the pixels falling within the ROI were used in the following analyses (e.g., REP assessment).The averaged REP values derived for each ROI were then correlated with the per-ROI-averaged Chl × LAI content (see more in Section 2.8).To derive the REP from the HyMap, the polynomial fitting (PF) method and Linear Extrapolation (LE) were employed as well as the 4-Point Linear Interpolation (4PLI) in its general form (4PLI) and in the form tuned for the spectral resolution of HyMap (4PLIH) [31].To derive the REP from the S2, the 4-Point Linear Interpolation method adjusted to the spectral resolution of S2 data was employed (4PLIS).Each algorithm is designated based on the following formulae: 2.6.1.Polynomial Fitting (PF) A polynomial function was least-squares-fitted to the reflectance spectrum between the wavelengths, corresponding to the minimum reflectance in the red and the maximum NIR (shoulder) reflectance: where: λ is the band between 670 nm to 780 nm.The REP is determined using the first derivative of the function used for approximating the source spectrum.
The advantage of this method is that it is only necessary to set the spectral range.

Linear Extrapolation (LE)
This method is based on the supposition of the REP location at the intersection of the Far Red Line (FRL) and Near Infrared Line (NIL) defined by the first derivative of the reflectance at 680 and 700 nm (FRL) or 725 and 760 nm (NIL).For the HyMap spectral resolution can be defined as: NIL is defined as: NIL = a 2 λ + b 2 (7) where: a 1 and a 2 are the gains of FRL resp.NIL; b 1 and b 2 are offsets of FRL resp.NIL; D 680 , D 700 , D 725 and D 760 are first derivatives of the reflectance at 680, 700, 725 and 760 nm, respectively.REP is then defined as:

4-Point Linear Interpolation (4PLI)
It is assumed that the reflectance curve at the red-edge can be approximated by a straight line centred near the mid-point between the reflectance at 780 nm and 670 nm.As such, it is, in a general form, defined as: The form adjusted to the HyMap data resolution is defined as: Forests 2016, 7, 226 8 of 17 The form adjusted to the S2 data resolution is defined as: where: R 670 , R 780 , R 700 and R 740 are the reflectance values at 670, 780, 700 and 740 nm, respectively.

Comparison between REP IMG and REP SIM
REP SIM and REP IMG were then compared for each individual REP extraction method (4PLI, PF, LE, 4PLIH and 4PLIS).The fits of the REP SIM and REP IMG values were expressed in terms of RMSE and coefficient of determination (R 2 ).The RMSE was calculated based on the following equation: where: REP IMGi is the REP value extracted from the image data for the forest stand i; REP SIMi is the REP value extracted from the FRT simulated spectra for the forest stand i; n is the number of forest stands (in this case n = 4).The calculated REP SIM and REP IMG for HyMap and S2 spectral resolution were compared with the REP values considered as the comparative standard (REP STD ).This comparative standard was defined as the REP value extracted from the original FRT simulations using the PF method because: (i) the spectral resolution of the original FRT simulations (5 nm) is higher than the spectral resolution of the HyMap data and much higher than the spectral resolution of the S2 data; (ii) the reflectance curves can be approximated by a polynomial function with very high accuracy; (iii) the PF is the only fully continuous REP extraction method in contrast to other methods based on discrete approaches using various types of inter or extrapolations.
Comparing the REP SIM values extracted from the HyMap and S2 level spectra with the standard REP STD , we can observe the influence of the spectral resolution on the REP calculation accuracy: the REP SIM values are calculated from the simulated spectra that differs from the spectra used for calculating the standard value only by the spectral resolution.Comparing the REP IMG values with the standard enables a study of precisely how it is possible to find the real REP value using the real imagery and which REP extraction method is optimal in each case.

Vegetation Properties and REP Correlation
As described in Section 2.6, the REP values derived for each ROI were then correlated with the per-ROI-averaged Chl × LAI content (Figure 2).In this specific case, for each of 4 sites one representative Chl × LAI value was obtained, an average value of 10 to 15 sampled trees (Erika: 10 trees, Habartov: 15 trees, Mezihorská: 15 trees and Studenec: 10 trees).Due to the low number of forest stands studied, the relationship and the strength of the correlations performed are described using R 2 and RMSE, respectively.First, the correlation was performed between the REP SIM values of the particular forest stands and Chl × LAI, which shows how the correlations should theoretically look.Then, the REP derived from the image data (REP IMG for HyMap and simulated S2 data respectively) were correlated with the corresponding Chl × LAI values.

Comparison of REP STD , REP IMG and REP SIM Values
The REP STD was compared with REP SIM and the corresponding REP IMG at both the HyMap (5-m spatial resolution, PF, LE, 4PLI and 4PLIH methods) and S2 (20-m spatial resolution, 4PLIS method) levels (Figure 3).The relative fitting of the REP IMG and REP SIM can be considered satisfactory.A considerable difference can be seen between the REP SIM and REP IMG spectra when using the LE extraction method followed by the PF algorithm at all the study sites for both years (2009 and 2010).

Comparison of REPSTD, REPIMG and REPSIM Values
The REPSTD was compared with REPSIM and the corresponding REPIMG at both the HyMap (5-m spatial resolution, PF, LE, 4PLI and 4PLIH methods) and S2 (20-m spatial resolution, 4PLIS method) levels (Figure 3).The relative fitting of the REPIMG and REPSIM can be considered satisfactory.A considerable difference can be seen between the REPSIM and REPIMG spectra when using the LE extraction method followed by the PF algorithm at all the study sites for both years (2009 and 2010).Table 1 lists the RMSE values which describe the difference between the REP SIM and REP IMG values as well as R 2 , which explains the correlation between the REP IMG and REP SIM values obtained by applying various REP extraction methods in both years (2009 and 2010).Comparison of REP SIM /RE IMG values with REP STD yielded following results.For the REP SIM/STD values, the best fit was found for 4PLI and PLIH (the same R 2 ; however RMSE are generally higher for the 2010 dataset).For the PF and LE methods, lower R 2 and higher RMSE were achieved; in addition, much higher deviations between REP SIM/STD and REP IMG/STD were found.Furthermore, comparable results with those derived for HyMap (4PLI and PLIH) were achieved for simulated S2 data when employing the (4PLIS) method.The RMSE values obtained for the 2010 dataset were generally higher compared with the ones obtained from the 2009 dataset.This can be explained by a higher signal-to-noise-ratio (SNR), which characterized the 2009 HyMap dataset when compared to the 2010 HyMap data [60].
The REP values were then extracted from the HyMap imagery as well as from the simulated S2 image data using the 4PLIH and 4PLIS methods (Figure 4).The results displayed in Figure 4 show that the REP spatial gradients are in good agreement when using HyMap (5-m spatial resolution) and simulated S2 image data (20-m spatial resolution).In addition, The HyMap REP image was resampled to 20-m spatial resolution to be spatially comparable with the S2 level data and to be able to display a scatterplot between the REP values extracted from HyMap and simulated S2 datasets (Figure 5).
The REP values were then extracted from the HyMap imagery as well as from the simulated S2 image data using the 4PLIH and 4PLIS methods (Figure 4).The results displayed in Figure show that the REP spatial gradients are in good agreement when using HyMap (5-m spatial resolution) and simulated S2 image data (20-m spatial resolution).In addition, The HyMap REP image was resampled to 20-m spatial resolution to be spatially comparable with the S2 level data and to be able to display a scatterplot between the REP values extracted from HyMap and simulated S2 datasets (Figure 5).The scatterplot of the REP values extracted from the HyMap and simulated S2 datasets suggested that the difference between these two datasets had a linear character with a systematic offset (Figure 5).The scatterplot of the REP values extracted from the HyMap and simulated S2 datasets suggested that the difference between these two datasets had a linear character with a systematic offset (Figure 5).

Correlation of the REP Calculated by Different Methods with Chl × LAI Values
The REP values obtained by the described extraction methods were correlated with the corresponding Chl × LAI content.The R 2 of REPSIM-Chl × LAI and REPIMG-Chl × LAI are listed in Table 3.

Correlation of the REP Calculated by Different Methods with Chl × LAI Values
The REP values obtained by the described extraction methods were correlated with the corresponding Chl × LAI content.The R 2 of REP SIM -Chl × LAI and REP IMG -Chl × LAI are listed in Table 3.

Discussion
This study modelled the REP for a Norway spruce Forest canopy Reflectance and Transmittance (FRT) using RTM (resampled to HyMap and S2 spectral resolution) as well as calculated the REP from the real HyMap and simulated S2 data.The S2 data were simulated from the HyMap 2009 and 2010 datasets using the forward and backward simulation modules (S2eteS) [47]; during this transformation, the effect of the real atmosphere was added to the image data as well as the real property of the sensor (e.g., sensor noise and detector nonlinearity).Therefore, this approach allowed a close-to-real sensor simulation of the S2 data, rather than just simple spectral and spatial resampling.
The REP was calculated using the following methods: PF, LE and 4PLI and its optimized versions for the spectral setting of the HyMap (4PLIH) and S2 data (4PLIS).The REP values calculated from image-extracted spectra (REP IMG ) were compared with the REP values extracted from canopy spectra simulated by couplings of PROSPECT-5 and FRT RTMs (REP SIM ) as well as with the comparative standard (REP STD ).The best match between HyMap REP IMG and both REP SIM and REP STD was found in the case of the 4PLI method and the optimized version (4PLIH).Furthermore, comparable results with those derived for HyMap, were achieved for simulated S2 data when employing the (4PLIS) method to extract the REP.The deviations between the REP SIM and REP IMG can be caused by the following factors.Unlike the PROSPECT-5-simulated data, the image datasets, HyMap and simulated Sentinel-2, both contain noise (HyMap because it is a real aerial hyperspectral dataset and the simulated S2 data because sensor noise was added as well as the parameters of the atmosphere).Moreover, the leaf internal structure considered by this model does not match the real structure of coniferous needles [59].Consequently, there might be some differences between the measured needle spectra and the model simulated ones.These differences are then propagated through the shoot and canopy up-scaling process, which can result in a mismatch between the model-simulated and image-extracted vegetation spectra.In general, a better fit between REP IMG and REP SIM and REP STD (higher R 2 and lower RMSE) for both HyMap and the simulated S2 image data was found for the 2009 dataset when compared to the 2010 dataset.This can be explained by the higher signal-to-noise-ratio (SNR), which characterized the 2009 HyMap dataset [60].
When correlating the REP extracted from hyMap data using the above-mentioned methods with the Chl × LAI content, the highest R 2 were obtained using the 4PLI method together with its optimized versions (4PLIH).Furthermore, comparable correlation levels were obtained between the Chl × LAI and REP derived from the HyMap employing 4PLIH and the simulated S2 data when employing the 4PLIS.A similar trend was demonstrated by Delegido et al. [35] in the study employed over diverse agriculture sites.The Chl × LAI maps were derived for simulated S2 data and aerial Airborne Hyperspectral System (AHS) data employing the Normalized Area Over the reflectance Curve (NAOC) method.It was demonstrated that the estimated chlorophyll from simulated S2 image was practically the same as the one estimated from the original AHS image.Similarly, Clevers and Gitelson [39] demonstrated the significance of the red-edge bands of S2 sensor for estimation Chl content in crops and grasslands.These results show that due to the inclusion of red-edge bands, S2 has a great potential for delivering high-quality products (e.g., estimates of chlorophyll content).
The spatial resolution of either the HyMap (5-m spatial resolution) or S2 image data (20-m spatial resolution) did not allow a direct comparison between the sampling tree-tops and corresponding pixel values.Instead, it was possible to compare the tree groups/clusters with corresponding ROI's pixel groups of HyMap and simulated S2 data.To be able to obtain comparable statistical results, the same ROIs were used for both HyMap and S2 statistics.This clustering decreased the number of the available input data for the empirical modelling of Chl × LAI.As a result, the Chl sampling points were averaged per site as well as the REP values from the pixels falling within the ROI's (in total 4 sites were evaluated).On the other hand, the authors believe that it helped to minimize the errors resulting from false identifications between sampled tree-tops and corresponding pixels.Moreover, the values were more representative due to them being an average of 10-15 sampled trees.No strong statistical conclusions can be made using this limited dataset; however, the results demonstrate that the S2 REP is a promising tool for mapping Chl gradients and their spatial trends when assessing forests.However, further research focusing on the utilization of S2 REP to estimate forest biochemical parameters, which would also include an extensive validation campaign, is still needed and remains a challenge for the future.
The next step should also include testing new retrieval algorithms.As suggested by Delegido et al. [35], it would be interesting to apply and evaluate advanced non-parametric statistical models to S2 data, more sophisticated regression methods (e.g., stepwise multiple linear regression, principal component regression, partial least square regression [61,62] and new machine learning techniques such as Gaussian processes regression [63][64][65][66].However, to employ, evaluate and validate these techniques much more extensive input data and ground truth datasets are required, especially for forests.

Conclusions
In this forest-oriented study, it was demonstrated that it is possible to calculate the red-edge-position (REP) using the simulated Sentinel-2 (S2) image data with comparable accuracy when using HyMap hyperspectral data.The lowest differences in REP values calculated from image-extracted spectra and from theoretical RTM simulations were found for the 4PLI method including its HyMap (4PLIH) and S2 optimized versions (4PLIS).When assessing correlation levels between REP and the Chl × LAI contents the optimal solutions, in the case of HyMap, were found to be the 4PLI and 4PLIH methods and for simulated S2 image data the 4PLIS method.A comparable correlation level between REP and Chl × LAI was obtained when using simulated S2 image data despite their lower spatial and spectral resolution compared to airborne hyperspectral imagery, such as the HyMap data.Despite its simplicity, the 4PLI REP extraction technique demonstrated its potential usefulness for forest assessment based on both airborne hyperspectral (e.g., HyMap) as well as the space-borne S2 image data.Although S2 data have neither such high spatial nor spectral resolution as the data provided by current airborne hyperspectral sensors, they will make a major contribution to remote-sensing-based vegetation studies, especially because of their remarkably high temporal resolution and spatial coverage.
So far there are still very few studies devoted to the use of S2 data for estimating forest biochemistry/biophysics.These results demonstrate that the S2 REP has a potential to estimate chlorophyll gradients and to map spatial trends when assessing forests.However, further research focusing on the use of S2 REP to estimate forest biochemical parameters, which would also include an extensive validation campaign, is still needed and remains a challenge for the future.

Figure 2 .
Figure 2. Example (the Mezihorská site) showing the Regions Of Interest (ROI) and the sampling points overlaid on the hyMap (left) and simulated Sentinel 2 (right).At the top: per-ROI-averaged reflectance derived for HyMap (left) and simulated S2 image data (right) which was used in the following analyses.

Figure 2 .
Figure 2. Example (the Mezihorská site) showing the Regions Of Interest (ROI) and the sampling points overlaid on the hyMap (left) and simulated Sentinel 2 (right).At the top: per-ROI-averaged reflectance derived for HyMap (left) and simulated S2 image data (right) which was used in the following analyses.

Figure 3 .
Figure 3. REP averages per sites: The REP STD compared with REP SIM and the corresponding REP IMG at both the HyMap (5-m spatial resolution, PF, LE, 4PLI and PLIH methods) and S2 (20-m spatial resolution, 4PLIS method) levels.

Table 1 lists
the RMSE values which describe the difference between the REPSIM and REPIMG values as well as R 2 , which explains the correlation between the REPIMG and REPSIM values obtained by applying various REP extraction methods in both years (2009 and 2010).

Table 1 .
Comparison of RMSE and R 2 values of REP SIM and REP IMG (Calculated by 4PLI, PF, LE, 4PLIH and 4PLIS methods).The obtained results showed that the best fit (characterized by the lowest RMSE together with the highest R 2 ) between the REP IMG and REP SIM values was found when using the 4PLI extraction method (RMSE 2009 = 0.45 nm, R 2 2009 = 0.66; RMSE 2010 = 0.43 nm, R 2 2010 = 0.74) and its 4PLIH version optimized specifically for the HyMap dataset (RMSE 2009 = 0.42 nm, R 2 2009 = 0.66; RMSE 2010 = 0.38 nm, R 2 2010 = 0.74).The 4PLIS method optimized for S2 data yielded results, which are generally comparable with its hyperspectral versions 4PLI or 4PLIH, although the RMSE value was slightly higher in this case (RMSE 2009 = 0.75 nm, R 2 2009 = 0.67; RMSE 2010 = 1.11 nm, R 2 2010 = 0.69).Both the REP IMG and REP SIM values were compared with the REP STD values.The differences between the REP IMG (or REP SIM ) values and the corresponding REP STD values (REP SIM/STD and REP IMG/STD respectively) were expressed by RMSE and R 2 .Both RMSE and R 2 values for 2009 and 2010 can be found in Table 2.

Table 2 .
Comparison of RMSE and R 2 values of REP SIM/STD and REP IMG/STD (For HyMap: 4PLI, PF, LE, 4PLIH and for simulated S2 data: 4PLIS methods).

Table 3 .
R 2 and RMSE describing the relationship of the REPSIM/REPIMG and Chl × LAI values.

Table 3 .
R 2 and RMSE describing the relationship of the REP SIM /REP IMG and Chl × LAI values.The R 2 obtained for REP IMG -Chl × LAI ranged from 0.22 to 0.67.The most consistent results when compared REP IMG -Chl × LAI and REP SIM -Chl × LAI as well as the highest correlations between HyMap REP IMG and the Chl × LAI values were obtained using the 4PLI extraction method or its derivatives 4PLIH (Chl × LAI: RMSE 2009 = 44.05mg/m 2 , RMSE 2010 = 32.04 mg/m 2 ).Employing the 4PLIS method, comparable correlation levels with those derived for HyMap were obtained between the simulated S2 REP IMG and Chl × LAI (Chl × LAI: RMSE 2009 = 49.28 mg/m 2 , RMSE 2010 = 32.11mg/m 2 ).