Adjustment of Sentinel-2 Multi-Spectral Instrument ( MSI ) Red-Edge Band Reflectance to Nadir BRDF Adjusted Reflectance ( NBAR ) and Quantification of Red-Edge Band BRDF Effects

Optical wavelength satellite data have directional reflectance effects over non-Lambertian surfaces, described by the bidirectional reflectance distribution function (BRDF). The Sentinel-2 multi-spectral instrument (MSI) acquires data over a 20.6◦ field of view that have been shown to have non-negligible BRDF effects in the visible, near-infrared, and short wave infrared bands. MSI red-edge BRDF effects have not been investigated. In this study, they are quantified by an examination of 6.6 million (January 2016) and 10.7 million (April 2016) pairs of forward and back scatter reflectance observations extracted over approximately 20◦ × 10◦ of southern Africa. Non-negligible MSI red-edge BRDF effects up to 0.08 (reflectance units) across the 290 km wide MSI swath are documented. A recently published MODIS BRDF parameter c-factor approach to adjust MSI visible, near-infrared, and short wave infrared reflectance to nadir BRDF-adjusted reflectance (NBAR) is adapted for application to the MSI red-edge bands. The red-edge band BRDF parameters needed to implement the algorithm are provided. The parameters are derived by a linear wavelength interpolation of fixed global MODIS red and NIR BRDF model parameters. The efficacy of the interpolation is investigated using POLDER red, red-edge, and NIR BRDF model parameters, and is shown to be appropriate for the c-factor NBAR generation approach. After adjustment to NBAR, red-edge MSI BRDF effects were reduced for the January data (acquired close to the solar principal where BRDF effects are maximal) and the April data (acquired close to the orthogonal plane) for all the MSI red-edge bands.


Introduction
The Sentinel-2A and Sentinel-2B satellites carry the Multi Spectral Instrument (MSI) that sense medium-resolution multi-spectral data designed for global land surface monitoring [1].Combined, the Sentinel-2A and -2B systems have a global median average revisit interval of five days [2].The MSI has 13 spectral bands ranging from the visible to the shortwave infrared (0.433 µm to 2.19 µm) with four 10 m visible and near-infrared bands, six 20 m red-edge, near-infrared (NIR), and short wave infrared bands, and three 60 m bands [1].The MSI red-edge bands are centered at 705 nm, 740 nm, and 783 nm, and are not present on the MSI heritage Landsat and Satellite Pour l'Observation de la Terre (SPOT) sensors.The red-edge encompasses the rapid change in vegetation reflectance from the red to the near infrared and is of interest for studies of chlorophyll and nitrogen content [3][4][5], leaf area index [6][7][8], and vegetation stress [9], and also may be useful for certain classification applications such as for burned area mapping [10,11].
The Sentinel-2A and -2B MSI remote sensing systems are in sun-synchronous polar orbits and acquire images over a 20.6 • field of view that results in directional reflectance effects over non-Lambertian surfaces.These effects, often described by the bidirectional reflectance distribution function (BRDF), can be non-trivial.In a recent study of 17 million pairs of forward and back scatter reflectance observations, pronounced visible, near-infrared, and short wave infrared BRDF effects were documented across the 290 km wide MSI swath [12].In this study, MSI red-edge BRDF effects are quantified using the same data that were extracted over 20 • × 10 • of southern Africa for two 10 day periods in January 2016 (acquired close to the solar principal plane when BRDF effects are maximal) and in April 2016 (acquired close to the orthogonal plane when BRDF effects are minimal).
The nadir BRDF-adjusted reflectance (NBAR), i.e., the reflectance at nadir (0 • view zenith) for a specified solar zenith, nominally has no BRDF effects and therefore is well suited for applications that require consistent reflectance in space and time.The Moderate Resolution Imaging Spectral radiometer (MODIS) NBAR product is generated systematically on a global basis using BRDF model parameters that are derived by inversion against 16 days of MODIS observations [13].This approach will not work for Sentinel-2 MSI data because there are insufficient cloud-free MSI observations available every 16 days for BRDF model inversion [12] and over longer periods the surface usually changes, precluding reliable BRDF model parameter estimations [14,15].Recently, an empirical c-factor approach was published to adjust Landsat reflectance to NBAR [16].The approach is based on multiplying the reflectance by the ratio of the reflectance modeled using MODIS BRDF spectral model parameters for the Landsat observed view zenith and for a nadir view.Global fixed BRDF model parameters were used and were defined by averaging, for each MODIS reflective wavelength band, the MODIS 500 m BRDF product (MCD43A1) spectral model parameters [13] over the entire global land surface for 12 months [16].It was developed based on observations that BRDF shapes are sufficiently similar over the Landsat 15 • field of view that a single global fixed set of MODIS BRDF spectral model parameters can be used regardless of the land cover or condition [16].Consequently, Landsat NBAR can be derived using the c-factor approach for all the Landsat global long-term records and prior to the year 2000 availability of the MODIS BRDF product.The approach was subsequently demonstrated and used to correct Sentinel-2 MSI 10 m blue (490 nm), green (560 nm), red (665 nm), and near-infrared (NIR) (842 nm), in addition to the MSI 20 m short wave infrared (SWIR) (1610 nm and 2190 nm), reflectance to NBAR [12].The three MSI 20 m red-edge bands were not considered because there are no equivalent red-edge MODIS BRDF spectral model parameters.
In this study, the c-factor approach is applied to the three MSI red-edge bands and used to adjust the red-edge reflectance to NBAR.MSI 705 nm, 740 nm, and 783 nm red-edge BRDF parameters are derived by linear wavelength-based interpolation of the red (645 nm) and NIR (858 nm) fixed global MODIS BRDF spectral model parameters.The efficacy of the interpolation is investigated by comparison with Polarization and Directionality of the Earth's Reflectances (POLDER) red, red-edge and NIR BRDF spectral model parameters derived from a recently available POLDER BRDF database [17].The c-factor NBAR approach is demonstrated considering the same MSI data used in [12] and the view zenith variation in MSI red-edge reflectance quantified before (surface reflectance) and after BRDF correction (NBAR) are compared.The paper concludes with a discussion and implications of the study findings.

Sentinel-2A Data
The Sentinel-2A data used in [12] were used in this study.Specifically, two ten day periods of Sentinel-2A MSI data acquired close to the solar principal (January 2016) and orthogonal planes (April 2016) over approximately 20 • × 10 • of southern Africa were analyzed.In this way, MSI data sensed where BRDF effects are maximal (solar principal plane) and minimal (solar orthogonal plane) [18,19] were examined.In each 10-day period, adjacent laterally overlapping MSI orbit swaths sensed in the forward and backward scatter directions were examined.More than 6.6 million (January 2016) and 10.6 million (April 2016) pairs of reflectance observations sensed three or seven days apart over the same location in the forward and backscatter directions in overlapping Sentinel-2A orbit swaths were considered.The maximum view zenith difference was 23.86 • .The Sentinel-2A MSI data were projected into the MODIS sinusoidal projection, but first had to be registered due to misregistration that is evident between Sentinel-2A overlapping orbits [20,21].The top of atmosphere reflectance data were corrected to surface reflectance using the SEN2COR atmospheric correction software [22].Only pairs of forward and backward reflectance values that were cloud-and snow-free, unsaturated, and had no significant change in their three or seven day separation derived by the application of a blue band and normalized difference vegetation index filter [12] were used.
In this study, different to what was reported in [12], the three MSI red-edge bands were used and are illustrated as black lines in Figure 1.The MSI red-edge bands centered at 705 nm and 740 nm are 15 nm wide and the MSI red-edge band centered at 783 nm is 20 nm wide.Figure 1 also shows the MODIS red and NIR spectral response functions (dashed lines).
Remote Sens. 2017, 9, x FOR PEER REVIEW 3 of 17 in the forward and backward scatter directions were examined.More than 6.6 million (January 2016) and 10.6 million (April 2016) pairs of reflectance observations sensed three or seven days apart over the same location in the forward and backscatter directions in overlapping Sentinel-2A orbit swaths were considered.The maximum view zenith difference was 23.86°.The Sentinel-2A MSI data were projected into the MODIS sinusoidal projection, but first had to be registered due to misregistration that is evident between Sentinel-2A overlapping orbits [20,21].The top of atmosphere reflectance data were corrected to surface reflectance using the SEN2COR atmospheric correction software [22].Only pairs of forward and backward reflectance values that were cloud-and snow-free, unsaturated, and had no significant change in their three or seven day separation derived by the application of a blue band and normalized difference vegetation index filter [12] were used.
In this study, different to what was reported in [12], the three MSI red-edge bands were used and are illustrated as black lines in Figure 1.The MSI red-edge bands centered at 705 nm and 740 nm are 15 nm wide and the MSI red-edge band centered at 783 nm is 20 nm wide.Figure 1 also shows the MODIS red and NIR spectral response functions (dashed lines).

Fixed Global Annual MODIS Red and NIR BRDF Model Parameters
The fixed global annual MODIS red and NIR BRDF model parameters defined in [16] were used in this study (Table 1).They were defined by averaging the MODIS 500 m BRDF product spectral model parameters (MCD43A1) [13] over the entire global land surface for 12 months [16].Only the highest quality and snow-free MCD43A1 data were used.The MODIS red and NIR BRDF model parameters are used because they bound the MSI red-edge wavelengths.The MODIS red and NIR bands are 50 nm and 35 nm wide, respectively, and their spectral response functions do not overlap with the MSI red-edge spectral response functions (Figure 1).

Fixed Global Annual MODIS Red and NIR BRDF Model Parameters
The fixed global annual MODIS red and NIR BRDF model parameters defined in [16] were used in this study (Table 1).They were defined by averaging the MODIS 500 m BRDF product spectral model parameters (MCD43A1) [13] over the entire global land surface for 12 months [16].Only the highest quality and snow-free MCD43A1 data were used.The MODIS red and NIR BRDF model parameters are used because they bound the MSI red-edge wavelengths.The MODIS red and NIR bands are 50 nm and 35 nm wide, respectively, and their spectral response functions do not overlap with the MSI red-edge spectral response functions (Figure 1).The MCD43A1 MODIS spectral BRDF model parameters f iso (λ), f vol (λ), and f geo (λ) are defined for each MODIS band by the semi-empirical Ross-Thick/Li-Sparse-Reciprocal BRDF model [13,24] and may be used to predict reflectance for any view and solar geometry as: where ρ λ Ω, Ω is the predicted reflectance for the viewing vector Ω (view zenith and azimuth angles) and solar illumination vector Ω' (solar zenith and azimuth angles); K vol (Ω,Ω') and K geo (Ω,Ω') are volumetric scattering and geometric-optical model kernels, respectively, which only depend on the sun-view geometry; and f iso (λ), f vol (λ), and f geo (λ) are the MODIS spectral BRDF model parameters.

Spectral BRDF Parameters Derived from the POLDER BRDF Database
MODIS has no red-edge bands suitable for land surface monitoring.In order to check the methodology to interpolate MODIS red-edge BRDF parameters (Section 3.2), an independent analysis was undertaken using different sensor data.There are only a small number of global coverage satellite sensors with a dedicated multi-angular sensing capability, such as the Multi-angle Imaging SpectroRadiometer (MISR) [25] and the Polarization and Directionality of the Earth's Reflectances (POLDER) [26].POLDER was selected as it has red-edge bands and the multi-angular data can be used to invert semi-empirical BRDF models [27].POLDER measures the radiance in nine spectral bands (443 nm to 1020 nm) with nadir pixel dimensions of 6 × 7 km and provides up to 16 observations of the same target per orbit that are sensed over a range of view zenith angles up to nearly 70 • [26].
Recently, a POLDER BRDF database was released that provides high-quality POLDER observations sensed over globally distributed sites that are free from significant aerosol and cloud contamination [17].The database has an interface that enables interactive analysis of the POLDER observations and the inversion of BRDF models.All of the available snow-, ice-, and water-free observations in the database for 2008, i.e., more than 900,000 observations covering 8654 globally distributed sites, were considered.At each site, the snow-, ice-, and water-free POLDER observations for each calendar month of 2008 (typically at least 40 observations per site per month) were inverted against the Ross-Thick/Li-Sparse-Reciprocal model [13,24] to derive spectral BRDF model parameters f iso (λ), f vol (λ), and f geo (λ) for each of the six POLDER bands in the database.Only site observations that resulted in positive BRDF model parameter values were retained following suggestions in previous MODIS BRDF research [13,28].In addition, to quality check the BRDF model parameters, the white sky albedo, WSA(l), was derived for each band to ensure legal 0 ≤ WSA(l) ≤ 1 values [28,29].Then, the mean value of each BRDF model parameter for all the globally distributed POLDER sites was derived, i.e., in the same way that we derived the global fixed MODIS BRDF parameters [16].The results are shown in Table 2.Only the BRDF model parameters for the POLDER red (670 nm), red-edge (765 nm), and NIR (865 nm) bands were used in this study.The values for the other three POLDER bands are shown to provide context and in case researchers wish to use them for other purposes.The POLDER red (670 nm), red-edge (765 nm), and NIR (865 nm) bands have 15 nm, 38 nm, and 33.5 nm band widths and non-overlapping spectral response functions [30].
Table 2. Fixed global annual spectral BRDF model parameters for the six POLDER bands derived from the POLDER database [17].The 765 nm band is a red-edge band.

Sentinel-2 MSI Red-Edge Band NBAR Derivation
The nadir BRDF-adjusted reflectance (NBAR) was derived for each MSI red-edge band (705 nm, 740 nm, or 783 nm) using the c-factor method as [12]: where ρ λ is the Sentinel-2A MSI reflectance sensed with view zenith θ Sentinel v and solar zenith θ Sentinel s for red-edge band wavelength λ, and NBAR λ is the NBAR equivalent estimated for a solar zenith angle (θ s ) and nadir view zenith (θ ν = 0) and k is the local solar zenith (set as the average solar zenith of the pair of forward and backward scattering observations).The ρ λ values in (2) are defined as (1) using red-edge band interpolated MODIS spectral BRDF model parameters (Section 3.2).

Derivation of Sentinel-2 MSI Red-Edge Band BRDF Spectral Model Parameters by Linear Interpolation between the Red and NIR MODIS BRDF Spectral Model Parameters
In order to compute MSI red-edge NBAR, as (2), red-edge BRDF spectral model parameters are required.They were derived by linear interpolation between the red and NIR MODIS BRDF spectral model parameters.The BRDF shapes at red and NIR wavelengths are often quite different and typically red wavelengths have higher relative reflectance variation with respect to view and solar angle than the NIR [31,32].The wavelength dependence of BRDF is not well studied, and only a small number of studies have considered red-edge wavelengths [33][34][35].Typically, the BRDF of terrestrial surfaces can be described by dome or bowl shapes with greater reflectance in the back scatter direction [36][37][38] that is primarily due to shadow-hiding that occurs at the sand-sized microscale over soil surfaces [39] and at leaf and crown scales over vegetated surfaces [40].
In this study, we implemented a simple linear wavelength dependent interpolation of the semi-empirical Ross-Thick/Li-Sparse-Reciprocal red and NIR BRDF model parameters.The BRDF model implicitly models surface heterogeneity (assuming linear scaling of reflectance) [13,24,41] and it is meaningful to interpolate the BRDF model parameters and then predict reflectance from the interpolated parameter values, e.g., [42,43].
The Sentinel-2 MSI red-edge band (705 nm, 740 nm, 783 nm) BRDF model parameters values were derived by linear interpolation between the bounding MODIS red and NIR band fixed global annual BRDF model parameters as: where f * i (λ) is the interpolated BRDF model parameter i for the MSI red-edge band centered on wavelength λ; f i (645) and f i (858) are the MODIS red and NIR fixed global annual BRDF model parameters (Table 1), respectively; and i is one of three parameters (iso, geo, or vol).
To gain insights into the efficacy of the above rather simple interpolation, the same approach was applied to the POLDER red and NIR bands as: where f * i (λ) is the interpolated BRDF model parameter i for the POLDER red-edge band centered at 765 nm; f i (670) and f i (865) are the POLDER red and NIR fixed global annual BRDF model parameters (Table 2), respectively; and i is one of three parameters (iso, geo, or vol).Plots of predicted reflectance derived as (1), and c-factors derived as (2), were generated using the fixed global annual POLDER red-edge BRDF model parameters (Table 2) and using the interpolated red-edge BRDF model parameters derived as (4).In addition, plots of predicted reflectance and c-factors for the POLDER red and NIR BRDF model parameters (Table 2) were generated.The plots were generated as a function of view zenith angle and for different fixed solar zenith angles (0 • , 30 • , or 45 • ).
3.3.Quantification of Sentinel-2 MSI Red-Edge Directional Reflectance Effects and Reduction of Directional Reflectance Effects in MSI Red-Edge NBAR More than 6.6 million (January 2016) and 10.6 million (April 2016) pairs of reflectance observations sensed three or seven days apart over the same location in the forward and backscatter directions in overlapping Sentinel-2A orbit swaths were considered.Red-edge band reflectance differences between adjacent overlapping Sentinel-2A swaths in the forward and backward scatter directions were quantified for both before (surface reflectance) and after BRDF correction (NBAR).Following the same approach as [12], the mean absolute difference, and also the relative absolute percentage difference, were derived for every pair of overlapping observations from forward and backward reflectance values as: where ∆ρ λ and ∆ρ λ * are the mean absolute and the relative absolute percentage reflectance differences, respectively; and ρ f orward,λ i and ρ backward,λ i are a pair of forward and backward reflectance values, respectively, and there are n pairs for Sentinel-2A wavelength λ.The same equations were applied for the pairs of reflectance adjusted to NBAR.Ordinary least squares (OLS) regression fits of the differences between the observed reflectance pairs and also the NBAR pairs as a function of the view zenith angle were generated independently for each Sentinel-2A red-edge band.The OLS slope term multiplied by the maximum observed MSI view zenith range (23.86 • ) was derived to quantify the average difference between surface reflectance in the forward and backward scatter directions at the MSI scan edges and, for brevity, is termed the B-F difference [12,16].

Sentinel-2 MSI Red-Edge Band NBAR Derivation
Table 3 shows the Sentinel-2A MSI red-edge band BRDF model parameters values derived by linear interpolation, as (3), of the MODIS red and NIR BRDF model parameters.The interpolated red-edge band parameter values are intermediate between the MODIS red and NIR BRDF parameter values and increase with wavelength because the MODIS BRDF parameter values are greater at the NIR than at the red wavelengths (Table 1).The POLDER red-edge (765 nm) parameters derived by linear interpolation, as (4), of the POLDER red and NIR fixed global annual spectral BRDF model parameters (Table 2) were f_iso = 0.2040, f_geo = 0.0299, and f_vol = 0.1094.The POLDER and Sentinel-2 MSI red-edge band BRDF spectral model parameters are expected to be different because they were defined for different red-edge wavelengths and because of POLDER and MODIS sensor differences.
Figure 2 shows predicted POLDER reflectance, derived as (1) over a ±60 • view zenith angle range and for three fixed solar zenith angles, generated using the fixed global annual POLDER red, red-edge, and NIR BRDF model parameters (Table 2) and using the interpolated POLDER red-edge parameters derived as (4) from the red and NIR POLDER parameters.The reflectance values are symmetrical for solar zenith = 0 • around nadir (view zenith = 0 • ) and are bowl shaped with increasing reflectance in the backscatter (positive view zenith) at higher solar zenith angles, with an evident hot-spot when the view and solar zenith angles coincide.Similar predicted reflectance anisotropy was observed for reflectance predicted from the MODIS red and NIR BRDF parameters (Table 1) reported in [16].Of most relevance here is that the predicted POLDER 765 nm reflectance (black) and the interpolated equivalent reflectance (orange) are dissimilar.The interpolated red-edge reflectance (orange) is close to the mean of the red (red) and NIR (blue) predicted reflectances because the POLDER red-edge band central wavelength (765 nm) is about mid-way between the POLDER red (670 nm) and NIR (865 nm) central wavelengths.However, the POLDER predicted red-edge reflectance (black) is much closer to the predicted NIR reflectance (blue).Evidently, the linear interpolation approach is not particularly appropriate, i.e., the POLDER wavelength dependence of the BRDF in the red-edge cannot be modelled reliably by linear interpolation of red and NIR BRDF parameters, and therefore interpolation of the MODIS red and NIR BRDF parameters to Sentinel-2 MSI red-edge wavelengths will also be unreliable.Figure 2 shows predicted POLDER reflectance, derived as (1) over a ±60° view zenith angle range and for three fixed solar zenith angles, generated using the fixed global annual POLDER red, red-edge, and NIR BRDF model parameters (Table 2) and using the interpolated POLDER red-edge parameters derived as (4) from the red and NIR POLDER parameters.The reflectance values are symmetrical for solar zenith = 0° around nadir (view zenith = 0°) and are bowl shaped with increasing reflectance in the backscatter (positive view zenith) at higher solar zenith angles, with an evident hot-spot when the view and solar zenith angles coincide.Similar predicted reflectance anisotropy was observed for reflectance predicted from the MODIS red and NIR BRDF parameters (Table 1) reported in [16].Of most relevance here is that the predicted POLDER 765 nm reflectance (black) and the interpolated equivalent reflectance (orange) are dissimilar.The interpolated red-edge reflectance (orange) is close to the mean of the red (red) and NIR (blue) predicted reflectances because the POLDER red-edge band central wavelength (765 nm) is about mid-way between the POLDER red (670 nm) and NIR (865 nm) central wavelengths.However, the POLDER predicted red-edge reflectance (black) is much closer to the predicted NIR reflectance (blue).Evidently, the linear interpolation approach is not particularly appropriate, i.e., the POLDER wavelength dependence of the BRDF in the red-edge cannot be modelled reliably by linear interpolation of red and NIR BRDF parameters, and therefore interpolation of the MODIS red and NIR BRDF parameters to Sentinel-2 MSI red-edge wavelengths will also be unreliable.Figure 3 shows the c-factor values derived as ( 2) for the same wavelengths, BRDF model spectral parameters, and solar zenith angles as seen in Figure 2. The c-factor values are plotted over the Sentinel-2 view zenith angle range.They vary linearly with view zenith and are similar to each other.This is because the BRDF shapes are similar over the Sentinel-2 field of view (Figure 2) and the c-factor provides a multiplicative scaling factor that is dependent on the shape but not the magnitude of the BRDF.The red-edge c-factors derived using the fixed global annual POLDER red-edge parameters and the interpolated parameters are almost the same, with less than a 1% difference at the greatest Sentinel-2 view zenith angle.Evidently, the difference in the magnitude of the predicted POLDER red-edge (black) and interpolated red-edge (orange) reflectance evident in Figure 2 is unimportant for red-edge BRDF normalization using the c-factor approach.
Remote Sens. 2017, 9, x FOR PEER REVIEW 8 of 17 Sentinel-2 view zenith angle range.They vary linearly with view zenith and are similar to each other.This is because the BRDF shapes are similar over the Sentinel-2 field of view (Figure 2) and the c-factor provides a multiplicative scaling factor that is dependent on the shape but not the magnitude of the BRDF.The red-edge c-factors derived using the fixed global annual POLDER red-edge parameters and the interpolated parameters are almost the same, with less than a 1% difference at the greatest Sentinel-2 view zenith angle.Evidently, the difference in the magnitude of the predicted POLDER red-edge (black) and interpolated red-edge (orange) reflectance evident in Figure 2 is unimportant for red-edge BRDF normalization using the c-factor approach.2) and using the interpolated POLDER red-edge 765 nm (orange) BRDF model parameters.Shown for Sentinel-2 MSI ± 10.3° view zenith angle range and for three fixed solar zenith angles (0°, 30°, 45°).The orange and blue lines are almost identical and the blue lines are plotted over the orange lines.

Quantification of Sentinel-2 MSI Red-Edge Band Bi-Directional Reflectance Effects
Figures 4 and 5 show scatterplots of the MSI red-edge band surface reflectance differences as a function of view zenith for the 10-days of January and April data, respectively.In each period, pairs of MSI observations sensed in the forward and backward scatter directions over the same location were compared.The plots were generated using the same method as [12], but are shown here for the MSI red-edge bands.The figures show the surface reflectance differences for each pair plotted as a function of the view zenith angle.The view zenith BRDF effect across the Sentinel-2 field of view is evident for all three red-edges bands and in both months.The view zenith BRDF effect is greater in January than in April because the January data were sensed close to the solar principal plane when BRDF effects are pronounced and the April data were sensed closer to the orthogonal plane where

Quantification of Sentinel-2 MSI Red-Edge Band Bi-Directional Reflectance Effects
Figures 4 and 5 show scatterplots of the MSI red-edge band surface reflectance differences as a function of view zenith for the 10-days of January and April data, respectively.In each period, pairs of MSI observations sensed in the forward and backward scatter directions over the same location were compared.The plots were generated using the same method as [12], but are shown here for the MSI red-edge bands.The figures show the surface reflectance differences for each pair plotted as a function of the view zenith angle.The view zenith BRDF effect across the Sentinel-2 field of view is evident for all three red-edges bands and in both months.The view zenith BRDF effect is greater in January than in April because the January data were sensed close to the solar principal plane when BRDF effects are pronounced and the April data were sensed closer to the orthogonal plane where BRDF effects are less pronounced.This is reflected in the slopes of the OLS regression lines that are significant for all bands and both months (p-values < 0.0001) and have higher r 2 values (>0.7) for January than for April (>0.3) (Table 4).The B-F difference, i.e., the OLS slope term multiplied by the maximum observed view zenith range (23.86 • ), quantifies the average difference between surface reflectance in the forward and backward scatter directions at the MSI scan edges.The B-F differences for all three red-edge bands are about 0.07 to 0.08 (January) and 0.02 to 0.03 (April) (Table 4).
Remote Sens. 2017, 9, x FOR PEER REVIEW 9 of 17 BRDF effects are less pronounced.This is reflected in the slopes of the OLS regression lines that are significant for all bands and both months (p-values < 0.0001) and have higher r 2 values (>0.7) for January than for April (>0.3) (Table 4).The B-F difference, i.e., the OLS slope term multiplied by the maximum observed view zenith range (23.86°), quantifies the average difference between surface reflectance in the forward and backward scatter directions at the MSI scan edges.The B-F differences for all three red-edge bands are about 0.07 to 0.08 (January) and 0.02 to 0.03 (April) (Table 4).         4 and 5; OLS regression coefficient of determination (r 2 ), the OLS regression F-test p-value, and the B-F difference.A total of 6,600,685 and 10,656,197 pairs of forward and back scatter surface reflectance values were considered in January and April, respectively.4).4).and 7 show the surface NBAR differences for the January and April data, respectively.Clearly, the normalization of the surface reflectance to NBAR reduces the BRDF effects that are less apparent in the NBAR scatterplots than the corresponding surface reflectance differences (Figures 4  and 5).However, as previously reported for the other MSI bands [12], the BRDF normalization is not perfect.Although the April NBAR exhibit almost no residual BRDF effects with OLS regression lines that are near horizontal (Figure 7), the January NBAR data have evident small residual BRDF effects (Figure 6).The OLS regressions all have low r 2 values (<0.21), but are statistically significant (p < 0.0001) (Table 6).The NBAR B-F differences are small, ~0.02 for the January data and <0.007 for the April data.Table 6.Summary of the ordinary least squares (OLS) linear regressions of the surface NBAR differences illustrated in Figures 6 and 7; OLS regression coefficient of determination (r 2 ), the OLS regression F-test p-value, and the B-F difference.A total of 6,600,685 and 10,656,197 pairs of forward and back scatter surface reflectance values were considered in January and April, respectively.The OLS slope magnitudes and the B-F NBAR differences (Table 6) are always smaller than the corresponding values for the surface reflectance data (Table 4).Table 7 summarizes the mean absolute (∆ρ λ ) and relative absolute percentage (∆ρ λ * ) differences between the pairs of forward and backward scatter surface NBAR values.A comparison of these difference statistics with the corresponding surface reflectance difference statistics (Table 5) reveals that, for all bands and for both months, the NBAR differences are smaller than the corresponding surface reflectance differences.The ∆ρ λ * surface NBAR values are about half the ∆ρ λ * surface reflectance values for January and vary among the three red-edge bands from 6.0% to 7.5%.In April, when the BRDF effects were less pronounced, the ∆ρ λ * surface NBAR values vary among the three red-edge bands from 8.6% to 9.7% (Table 7) and are smaller than the ∆ρ λ * surface reflectance values (Table 5) by less than 1%.Evidently, the MODIS BRDF parameter c-factor approach systematically reduces Sentinel-2 red-edge BRDF effects.6).

Discussion
BRDF effects are known to be present in most optical wavelength data including the red-edge [33][34][35].A previous study, undertaken using the same data and analytical techniques as this one, documented pronounced visible, near-infrared, and short wave infrared BRDF effects across the Sentinel-2A MSI field of view [12].In this study, we found similar BRDF effects in the three MSI red-edge bands, i.e., in MSI bands 5 (705 nm), 6 (740 nm), and 7 (783 nm).The B-F difference (that quantifies the average difference between surface reflectance in the forward and backward scatter directions at the MSI scan edges) varied among the three MSI red-edge bands from about 0.07 to 0.08 (January) and from about 0.02 to 0.03 (April).To provide context, the median red-edge surface reflectance for the Southern Africa study data was 0.22 (band 5), 0.25 (band 6), and 0.27 (band 7) for the January data and 0.18 (band 5), 0.25 (band 6), and 0.28 (band 7) for the April data.Further, the relative absolute percentage differences (Δρ λ ̅̅̅̅̅ * ) between the pairs of forward and backward scatter surface reflectance values varied among the red-edge bands from 13.3% to 14.9% (January) and 9.4% to 10.5% (April).The magnitude of these red-edge BRDF effects are not insignificant, and are comparable in magnitude to Landsat atmospheric correction errors [44], and will likely be an issue for some Sentinel-2 MSI red-edge band applications.For example, the Sentinel-2 red-edge bands are of interest for studies of chlorophyll and nitrogen content, but statistical [45,46] and model inversion [47] approaches will require minimization or incorporation of BRDF variations into the estimation approach.6).

Discussion
BRDF effects are known to be present in most optical wavelength data including the red-edge [33][34][35].A previous study, undertaken using the same data and analytical techniques as this one, documented pronounced visible, near-infrared, and short wave infrared BRDF effects across the Sentinel-2A MSI field of view [12].In this study, we found similar BRDF effects in the three MSI red-edge bands, i.e., in MSI bands 5 (705 nm), 6 (740 nm), and 7 (783 nm).The B-F difference (that quantifies the average difference between surface reflectance in the forward and backward scatter directions at the MSI scan edges) varied among the three MSI red-edge bands from about 0.07 to 0.08 (January) and from about 0.02 to 0.03 (April).To provide context, the median red-edge surface reflectance for the Southern Africa study data was 0.22 (band 5), 0.25 (band 6), and 0.27 (band 7) for the January data and 0.18 (band 5), 0.25 (band 6), and 0.28 (band 7) for the April data.Further, the relative absolute percentage differences (∆ρ λ * ) between the pairs of forward and backward scatter surface reflectance values varied among the red-edge bands from 13.3% to 14.9% (January) and 9.4% to 10.5% (April).The magnitude of these red-edge BRDF effects are not insignificant, and are comparable in magnitude to Landsat atmospheric correction errors [44], and will likely be an issue for some Sentinel-2 MSI red-edge band applications.For example, the Sentinel-2 red-edge bands are of interest for studies of chlorophyll and nitrogen content, but statistical [45,46] and model inversion [47] approaches will require minimization or incorporation of BRDF variations into the estimation approach.For Sentinel-2 red-edge applications that consider subtle changes in reflectance, there is a need for the correction of BRDF effects.Recently, an empirical c-factor approach was demonstrated to adjust MSI reflectance to a nadir view and specified solar zenith, i.e., to estimate nadir BRDF-adjusted reflectance (NBAR) [12].The approach was applied to the MSI visible, near-infrared, and short wave infrared reflectance and shown to reduce BRDF effects.It was developed based on observations that BRDF shapes are sufficiently similar over the narrow field of view that a single global fixed set of MODIS BRDF spectral model parameters can be used to provide a first-order BRDF correction with little sensitivity to the land cover type, condition, or surface disturbance [16].As there are no MODIS red-edge BRDF model parameters, a linear wavelength interpolation of the MODIS red and NIR global fixed BRDF parameters to the MSI red-edge band wavelengths was implemented.The developed method worked successfully.For all the MSI red-edge bands and for both months of data considered, the NBAR forward and backscatter differences were smaller than the corresponding surface reflectance differences.The approach can be applied in a computationally efficient manner globally.As noted in [16], the approach should not be used to derive surface albedo, which would require a more accurate characterization of the surface BRDF than provided by fixed BRDF parameters.
The red-edge BRDF normalization to NBAR was not perfect because of factors including atmospheric correction and sensor calibration errors [12].In addition, a concern for this study is the appropriateness of the red-edge BRDF parameter generation by linear wavelength interpolation of the MODIS red and NIR BRDF parameters.To check that this simple interpolation approach was sufficient, an independent evaluation was undertaken using POLDER red, red-edge, and NIR fixed BRDF parameters.POLDER data were used because POLDER has a red-edge band with a multi-angular observation capability and the recent POLDER BRDF database makes access to POLDER data straightforward [17].The POLDER BRDF parameters derived from a year of globally distributed sites were used to compute global fixed red, red-edge, and NIR BRDF parameters.The linear wavelength interpolation approach was applied to the POLDER red and NIR fixed global BRDF parameters to derive interpolated POLDER red-edge BRDF parameters.Predicted reflectance values for a range of view zenith and solar zenith angles derived using the linearly interpolated and the fixed red-edge parameters were compared.The reflectances had similar BRDF shapes but quite different magnitudes, indicating that the POLDER wavelength dependence of the BRDF in the red-edge cannot be modelled reliably by the linear interpolation of red and NIR BRDF parameters.However, the corresponding POLDER derived c-factors were linear and similar over the Sentinel-2 field of view, i.e., indicating that the interpolation approach is sufficient for the purposes of this research.More sophisticated interpolation approaches, for example, by fitting a spline or sigmoid curve to the MODIS visible, NIR, and SWIR BRDF parameters could be developed.However, in the absence of MODIS resolution multi-angular red-edge wavelength data required to test the interpolation, future research will need to proceed by modelling and/or by consideration of multi-spectral and multi-angular airborne and field measurements.
The majority of available medium resolution atmospheric correction algorithms assume a Lambertian surface and so the c-factor approach can be meaningfully applied to the resulting atmospherically corrected data.We note, however, that surface BRDF and atmospheric corrections can be coupled but that this requires knowledge of the atmospheric constituents and/or sufficient multi-angular observations for model inversion [37, [48][49][50], which are not available from Sentinel-2 data.Despite this, research on the use of more physically-based correction approaches, and ones that include corrections for topographic and adjacency effects, that are particularly evident in medium spatial resolution data [51,52], is recommended.
Finally, as observed in [12,16], the spatial scale difference between medium resolution Sentinel-2 or Landsat data and coarse resolution MODIS data means that the fixed global mean 12 month MODIS BRDF spectral model parameters provide an average representation of the local surface reflectance anisotropy.Certain geographic locations and times may have greater red-edge reflectance anisotropy than captured by the MODIS interpolated red-edge BRDF spectral model parameters.In addition, although more than 17 million pairs of forward and back scatter reflectance observations sensed close to the solar principal and orthogonal planes were considered, they did not include evergreen needleleaf and snow and ice land covers or observations sensed with solar zenith >50 • .Therefore, the efficacy of the red-edge normalization described here may not be globally representative.

Conclusions
Sentinel-2 MSI red-edge band BRDF effects were observed in more than 6.6 million (January 2016) and 10.6 million (April 2016) pairs of forward and back scatter reflectance observations extracted over approximately 20 • × 10 • of southern Africa.The BRDF effects were quantified and caused MSI red-edge surface reflectance variations up to 0.08 (reflectance units) across the 290 km wide swath.This variation is not negligible and may be an issue for MSI red-edge applications.The recently published MODIS BRDF parameter c-factor approach, developed to derive Landsat NBAR [16], and shown to also work for the generation of Sentinel-2 visible, near-infrared, and short wave NBAR data [12], was adapted in this study for MSI red-edge band implementation.The algorithm reduced BRDF effects for all three MSI red-edge bands and for both months of study data.Residual BRDF effects remained in the NBAR data, about 0.02 for the January data (that had more pronounced surface reflectance BRDF effects) and <0.007 for the April data.The MSI red-edge band BRDF parameters needed to implement the red-edge c-factor NBAR algorithm are provided (Table 3) so that users may implement the algorithm.

Figure 3
Figure3shows the c-factor values derived as (2) for the same wavelengths, BRDF model spectral parameters, and solar zenith angles as seen in Figure2.The c-factor values are plotted over the

Figure 3 .
Figure 3. c-factors (derived as 2) using the fixed global annual POLDER NIR 865 nm (blue), red-edge 765 nm (black), and red 670 nm (red) BRDF model parameters (Table2) and using the interpolated POLDER red-edge 765 nm (orange) BRDF model parameters.Shown for Sentinel-2 MSI ± 10.3° view zenith angle range and for three fixed solar zenith angles (0°, 30°, 45°).The orange and blue lines are almost identical and the blue lines are plotted over the orange lines.

Figure 3 .
Figure 3. c-factors (derived as 2) using the fixed global annual POLDER NIR 865 nm (blue), red-edge 765 nm (black), and red 670 nm (red) BRDF model parameters (Table 2) and using the interpolated POLDER red-edge 765 nm (orange) BRDF model parameters.Shown for Sentinel-2 MSI ± 10.3 • view zenith angle range and for three fixed solar zenith angles (0 • , 30 • , 45 • ).The orange and blue lines are almost identical and the blue lines are plotted over the orange lines.

Figure 4 .
Figure 4. Sentinel-2A MSI red-edge surface reflectance differences in the Southern Africa January swath image overlap zones plotted against view zenith for a total of 6,600,685 pairs of forward and back scatter surface reflectance values.The plot colors show the relative frequency of occurrence of similar difference values (with a log2 scale).The solid lines show ordinary least squares (OLS) linear regression fits of these data (see Table4).Results shown for MSI bands 5 (705 nm), 6 (740 nm), and 7 (783 nm).
Figure 4. Sentinel-2A MSI red-edge surface reflectance differences in the Southern Africa January swath image overlap zones plotted against view zenith for a total of 6,600,685 pairs of forward and back scatter surface reflectance values.The plot colors show the relative frequency of occurrence of similar difference values (with a log2 scale).The solid lines show ordinary least squares (OLS) linear regression fits of these data (see Table4).Results shown for MSI bands 5 (705 nm), 6 (740 nm), and 7 (783 nm).

Figure 4 .
Figure 4. Sentinel-2A MSI red-edge surface reflectance differences in the Southern Africa January swath image overlap zones plotted against view zenith for a total of 6,600,685 pairs of forward and back scatter surface reflectance values.The plot colors show the relative frequency of occurrence of similar difference values (with a log2 scale).The solid lines show ordinary least squares (OLS) linear regression fits of these data (see Table4).Results shown for MSI bands 5 (705 nm), 6 (740 nm), and 7 (783 nm).
Figure 4. Sentinel-2A MSI red-edge surface reflectance differences in the Southern Africa January swath image overlap zones plotted against view zenith for a total of 6,600,685 pairs of forward and back scatter surface reflectance values.The plot colors show the relative frequency of occurrence of similar difference values (with a log2 scale).The solid lines show ordinary least squares (OLS) linear regression fits of these data (see Table4).Results shown for MSI bands 5 (705 nm), 6 (740 nm), and 7 (783 nm).

Figure 5 .
Figure 5. Same as Figure 4 but for the April data, a total of 10,656,197 pairs of forward and back scatter surface reflectance values were considered.The solid lines show ordinary least squares (OLS) linear regression fits of these data (see Table4).

Figure 5 .
Figure 5. Same as Figure 4 but for the April data, a total of 10,656,197 pairs of forward and back scatter surface reflectance values were considered.The solid lines show ordinary least squares (OLS) linear regression fits of these data (see Table4).

Figure 6 .
Figure 6.Sentinel-2A MSI red-edge surface NBAR differences in the Southern Africa January swath image overlap zones plotted against view zenith for a total of 6,600,685 pairs of forward and back scatter surface NBAR values.The plot colors show the relative frequency of occurrence of similar difference values (with a log2 scale).The solid lines show ordinary least squares (OLS) linear regression fits of these data (see Table6).Results shown for MSI bands 5 (705 nm), 6 (740 nm), and 7 (783 nm).

Figure 6 .
Figure 6.Sentinel-2A MSI red-edge surface NBAR differences in the Southern Africa January swath image overlap zones plotted against view zenith for a total of 6,600,685 pairs of forward and back scatter surface NBAR values.The plot colors show the relative frequency of occurrence of similar difference values (with a log2 scale).The solid lines show ordinary least squares (OLS) linear regression fits of these data (see Table6).Results shown for MSI bands 5 (705 nm), 6 (740 nm), and 7 (783 nm).

Figure 7 .
Figure 7. Same as Figure 6 but for the April data, a total of 10,656,197 pairs of forward and back scatter surface NBAR values were considered.The solid lines show ordinary least squares (OLS) linear regression fits of these data (see Table6).

Figure 7 .
Figure 7. Same as Figure 6 but for the April data, a total of 10,656,197 pairs of forward and back scatter surface NBAR values were considered.The solid lines show ordinary least squares (OLS) linear regression fits of these data (see Table6).

Table 1 .
[16]d global annual MODIS spectral BRDF model parameters for the MODIS red and NIR bands (reproduced from Table5of[16]) used in this study.

Table 1 .
[16]d global annual MODIS spectral BRDF model parameters for the MODIS red and NIR bands (reproduced from Table5of[16]) used in this study.

Table 3 .
Sentinel-2 MSI red-edge band BRDF spectral model parameters derived by linear interpolation between the red and NIR MODIS BRDF spectral model parameters (

Table 5
summarizes the mean absolute (Δρ λ ̅̅̅̅̅ ) and relative absolute percentage (Δρ λ ̅̅̅̅̅ * ) differences between the pairs of forward and backward scatter surface reflectance values.The mean relative absolute difference values, (Δρ λ ̅̅̅̅̅ * ) are normalized by the reflectance magnitude which changes with wavelength and surface conditions.The Δρ λ ̅̅̅̅̅ * values vary among the bands from 13.3% to 14.9% (January) and 9.4% to 10.5% (April).Clearly, these results illustrate that the MSI red-edge bands have view zenith BRDF effects.

Table 5
January) and 9.4% to 10.5% (April).Clearly, these results illustrate that the MSI red-edge bands have view zenith BRDF effects.
summarizes the mean absolute (∆ρ λ ) and relative absolute percentage (∆ρ λ * ) differences between the pairs of forward and backward scatter surface reflectance values.The mean relative absolute difference values, (∆ρ λ * ) are normalized by the reflectance magnitude which changes with wavelength and surface conditions.The ∆ρ λ * values vary among the bands from 13.3% to 14.9% (

Table 4 .
Summary of the ordinary least squares (OLS) linear regressions of the surface reflectance differences illustrated in Figures4 and 5; OLS regression coefficient of determination (r 2 ), the OLS regression F-test p-value, and the B-F difference.A total of 6,600,685 and 10,656,197 pairs of forward and back scatter surface reflectance values were considered in January and April, respectively.

Table 4 .
Summary of the ordinary least squares (OLS) linear regressions of the surface reflectance differences illustrated in Figures

Table 5 .
Mean absolute reflectance differences (∆ρ λ ) (Equation (5)) and mean absolute relative percentage differences (∆ρ λ * ) (Equation (6)) between the pairs of forward and backward scatter Sentinel-2A surface reflectance values illustrated in Figures4 and 5.A total of 6,600,685 and 10,656,197 pairs of forward and back scatter surface reflectance values were considered in January and April, respectively.

Table 7 .
Mean absolute reflectance differences (∆ρ λ ) (Equation (5)) and mean absolute relative percentage differences (∆ρ λ * ) (Equation (6)) between the pairs of forward and backward scatter Sentinel-2A surface NBAR values illustrated in Figures6 and 7. A total of 6,600,685 and 10,656,197 pairs of forward and back scatter surface NBAR values were considered in January and April, respectively.

Table 7 .
Mean absolute reflectance differences (Δρ λ ) ̅̅̅̅̅̅ (Equation (5)) and mean absolute relative percentage differences (Δρ λ ̅̅̅̅̅̅ * ) (Equation (6)) between the pairs of forward and backward scatter Sentinel-2A surface NBAR values illustrated in Figures6 and 7. A total of 6,600,685 and 10,656,197 pairs of forward and back scatter surface NBAR values were considered in January and April, respectively.