Toward Landsat and Sentinel-2 BRDF Normalization and Albedo Estimation: A Case Study in the Peruvian Amazon Forest

The Amazon Remote sensing data analysis is the way to study such a large geographical extent during an extended period of the its increase in temporal resolution the combination with Landsat a strong emphasis has been put on exploiting these data. Though these satellites provide near nadir observations, surface reﬂectance time series are affected by illumination variability throughout the year. These effects can be corrected using a Bidirectional Reﬂectance developed a methodology to derive Landsat surface albedo and BRDF. It is based on the BRDF parameters from the MODerate Resolution Imaging Spectroradiometer (MODIS) which are disaggregated at Landsat spatial resolution (30 m). In this work, we apply the Franch et al. (2014a) method to normalize the surface reﬂectance for BRDF effects using the NASA’s Harmonized Landsat Sentinel-2 (HLS) product. We apply this method to the Tambopata region in Peru from 2013 to 2017 and validate it using ground-based albedometer measurements. The results show that the near infrared reﬂectance can increase up to 0.06 (20%) for low solar angles while the impact on the red range and the NDVI is minor ( < 0.01). The evaluation of the surface albedo against ﬁeld measurements shows an error of 0.01.

The Amazon forest has been the focus of study by the science community during the last few decades. Remote sensing data analysis is the only way to study such a large geographical extent during an extended period of time. Since the launch in 2015 of Sentinel 2 and its increase in temporal resolution through the combination with Landsat sensors, a strong emphasis has been put on exploiting these data. Though these satellites provide near nadir observations, surface reflectance time series are affected by illumination variability throughout the year. These effects can be corrected using a Bidirectional Reflectance Distribution Function (BRDF) model. Franch et al. (2014a) developed a methodology to derive Landsat surface albedo and BRDF. It is based on the BRDF parameters from the MODerate Resolution Imaging Spectroradiometer (MODIS) which are disaggregated at Landsat spatial resolution (30 m). In this work, we apply the Franch et al. (2014a) method to normalize the surface reflectance for BRDF effects using the NASA's Harmonized Landsat Sentinel-2 (HLS) product. We apply this method to the Tambopata region in Peru from 2013 to 2017 and validate it using ground-based albedometer measurements. The results show that the near infrared reflectance can increase up to 0.06 (20%) for low solar angles while the impact on the red range and the NDVI is minor (<0.01). The evaluation of the surface albedo against field measurements shows an error of 0.01.

INTRODUCTION
Surface albedo is an essential parameter which is important not only for developing climatic models, but also for Earth's energy balance studies. In the physical climate system, albedo determines the radiation balance of the surface and affects the surface temperature and boundary-layer structure of the atmosphere. Particularly, in the Amazon forest, deriving an accurate surface albedo is critical to study the impact of anthropogenic land use change and to successfully identify the vegetation changes and climatic consequences associated with them (Bala et al., 2007). An analysis of these changes requires a high temporal resolution, traditionally only achievable through coarse spatial resolution sensors. However, in areas with such a high frequency of cloud cover, high to moderate spatial resolution albedo products are becoming increasingly necessary. While there is an already established surface albedo product for MODIS at coarse spatial resolution (Schaaf et al., 2002), there is no such a global albedo product at high to moderate spatial resolution (Landsat-Sentinel-2 at 10-30 m). The reason behind this lies in the fact that derivation of albedo, estimated through the angular integration of the Bidirectional Reflectance Distribution Function (BRDF), usually requires an appropriate angular sampling of the surface reflectance.
The Harmonized Landsat/Sentinel-2 (HLS) project (Claverie et al., 2018) provides a surface reflectance product that combines observations from USGS/NASA's Landsat 8 and ESA's Sentinel-2 satellites at moderate spatial resolution. These satellites' sampling characteristics provide nearly constant observation geometry and low illumination variation within each scene, making the albedo computation a challenging exercise. Additionally, when a surface reflectance time series, combining measurements from both sensors, is computed, variations due to differences in view geometry between them arise. In extreme cases, the differences in the view angle for a ground target can increase up to 20 • . Variations in the seasonal illumination also impact the surface reflectance value. In order to correct such effects, the official HLS product provides a Nadir BRDF-Adjusted Reflectance (NBAR) product based on the BRDF normalization proposed by Roy et al. (2016). However, this correction is not optimal because the considered coefficients are constant and do not depend on the surface type or condition.
Previous studies have worked toward generating a medium to high spatial resolution albedo product. (Shuai et al., 2011) used the MODIS Bidirectional Reflectance Distribution Function (BRDF) MDC43 product at 500 m (Schaaf et al., 2002) to generate a Landsat surface albedo product that later (Shuai et al., 2011) was implemented back to the 80 s to the "pre-MODIS era." Their method presents a Root Mean Square Error (RMSE) generally lower than 0.03 when compared to in situ data, but it is mostly determined by a negative bias. Furthermore, this method does not provide a BRDF product at Landsat scale. Franch et al. (2014a), presented an alternative algorithm to calculate the Landsat surface albedo based on the BRDF parameters estimated from the MODerate Resolution Imaging Spectroradiometer (MODIS) Climate Modeling Grid (CMG) surface reflectance product (M{O,Y}D09) using the VJB method (Vermote et al., 2009). The validation of the results against field measurements showed an RMSE of 0.015 (7%). This work just analyzed the accuracy of the model for surface albedo estimation. However, this algorithm also calculates the BRDF parameters for each pixel, allowing the correction of directional effects through the computation of the normalized surface reflectance.
In this work, we apply the Franch et al. (2014a) method to the HLS product in the Tambopata region (Peru) from 2013 to 2017 and analyze its feasibility for BRDF normalization purposes. The method is validated using a ground-based albedometer that has been retrieving continuous measurements since 2013.

METHODOLOGY
We implement the Franch et al. (2014a) method to HLS data on the Amazon forest. We extract a 9 × 9 km area centered on the albedometer tower location (12.818 S, 69.281 W). All images are then masked for cloud mask and cloud shadow effects using the product's cloud mask. Next, we run an unsupervised classification over every scene using the Iterative Self-Organizing Data Analysis Technique Algorithm (ISODATA) developed by Ball and Hall (1965). This classification is used to unmix the BRDF parameters retrieved from MODIS at HLS spatial resolution (30 m). Following the Vermote et al. (2009) notations, the surface reflectance (ρ) is written as: where θ s is the sun zenith angle, θ v is the view zenith angle, φ is the relative azimuth angle, F 1 is the volume scattering kernel, based on the Rossthick function derived by Roujean et al. (1992) and F 2 is the geometric kernel, based on the Li-sparse model (Li and Strahler, 1986) but considering the reciprocal form given by Lucht (1998) and corrected for the Hot-Spot process proposed by Maignan et al. (2004). F 1 and F 2 are fixed functions of the observation geometry, but k 0 , k 1 , and k 2 are free parameters. Following this notation, we define V as k 1 /k 0 and R for k 2 /k 0 . In order to invert the MODIS BRDF parameters (V and R) we use the VJB method (Vermote et al., 2009;Franch et al., 2014b). We unmix the BRDF parameters to the HLS 30 m spatial resolution, by assuming that the surface reflectance of a MODIS pixel can be written as a weighted sum of n Landsat classes: where A x , B x and N x represent proportions of each the n classes within the x MODIS pixel. BRDF parameters for each class, which are unknowns in Equation (2), can be derived through matrix inversion. For further details on this methodology, we refer the reader to Franch et al. (2014a). We only invert k 1 and k 2 parameters since they describe the directional shape, while k 0 describes the reflectance magnitude. Considering the HLS surface With the BRDF parameters at HLS resolution, we apply Equation (4) for deriving the normalized surface reflectance (ρ N ) at 45 degrees of solar zenith angle and nadir observation (Vermote et al., 2009): Finally, using the BRDF parameters at HLS resolution we derive the surface albedo which is validated against field measurements. The surface albedo is estimated following the same methodology as the official MCD43 MODIS product (Schaaf et al., 2002). The albedo can be written by integrals of the BRDF model through the black-sky albedo and the white-sky albedo (Schaaf et al., 2002).

HLS Data
The harmonization process of the HLS product is based on a set of algorithms for atmospheric correction, cloud and cloud-shadow screening and co-registration of data from the Operational Land Imager (OLI) onboard Landsat-8 and the Multi-Spectral Instrument (MSI) onboard Sentinel-2. Though the HLS product is also corrected for spectral adjustment and BRDF normalization. However, in this work we considered uncorrected data to explore the feasibility of the proposed algorithm for BRDF normalization. For this study, we used all HLS cloud free data available until December 2017.

MODIS Data
We also processed MODIS surface reflectance Collection 6 data at 1 km spatial resolution (M{OY}D09GA; Vermote, 2015) over the same areas and time period to derive the BRDF parameters using the VJB method (Vermote et al., 2009;Franch et al., 2014b).

In-situ Data
The in-situ albedo was derived from a validation station (12.818 S, 69.281 W) located on a homogeneous area of dense tropical forest in the Tambopata region of the Peruvian Amazon. At the station, a Campbell Scientific NR01 4-component net radiometer installed on a 44-meter tower does continuous measurements of upward and downward shortwave radiation over the canopy and the 10-min average value is recorded. More details of the validation station and instruments can be found in Vihermaa et al. (2016). Data was acquired during the time period studied (2013-2017). However, due to technical problems with data storage devices, there was just simultaneous field data (to the satellite overpass) during 2017. All parameters show good interconsistency between both satellites' products. The figure also shows that these sensors provide data during the wet season. This is generally challenging at coarser spatial resolutions caused generally by the presence of scattered clouds, as we noticed while processing this dataset. The surface reflectance not normalized for BRDF effects (red) shows a clear dependency with the SZA, showing higher values for low solar angles and lower values for high SZA. This effect is not fully corrected with the BRDF-normalization of the HLS product (green), which after the correction still shows a "smile" effect. However, it is corrected on the BRDF-normalized data using Franch et al. (2014a) (black), which normalize all observations to SZA=45 • and nadir observation. Comparing the normalized and not normalized data, the angular effect is stronger in the NIR and the DVI, where the difference can get up to 0.06, than the RED band while the NDVI barely shows an impact, with differences up to 0.01. Finally, we evaluate the BRDF normalization by using the BRDF parameters to derive the HLS surface albedo. This product is averaged using 3 × 3 pixels window around the location of the albedometer. Table 1 shows the comparison of the satellite estimations with the field measurements during 2017 considering a total of 7 different days. The product shows an RMSE of 0.011, which is mostly determined by a positive bias. Despite this bias, this result shows a better agreement than the results in Franch et al. (2014a) that reported an RMSE of 0.015 when analyzing several land covers.

DISCUSSION AND CONCLUSIONS
The main purpose of this study is evaluating the feasibility of Franch et al. (2014a) method to normalize the BRDF effects on HLS data. This method was presented as a surface albedo algorithm for medium resolution data (Landsat). However, it also provides pixel-based BRDF parameters that in the original manuscript was presented as an intermediate result and was not further analyzed. The data analyzed in this work shows that BRDF effects, caused mainly by the solar zenith angle yearly evolution, impacts more strongly the NIR region and the DVI by increasing them up to 0.06 for low solar zenith angles. The RED spectral region and the NDVI show a minor effect with differences lower than 0.01. These effects are corrected when applying Franch et al. (2014a) method on the particular area analyzed, while the current HLS BRDF normalization shows a "smile" effect. This results demonstrates the potential application of this method to correct BRDF effects. Note that we apply the method on a particular area of the Peruvian Amazon forest which is representative of the Southwest Amazon moist forest ecoregion (Olson et al., 2001). Therefore, in future works we will apply this method to different land cover types. Additionally, this method shows potential to be applied in the operational context as the processing of a single HLS tile takes approximately 15 min using a computing machine with an Intel(R) Xeon(R) CPU E5-2650 0 @ 2.00 GHz processor. The evaluation of the BRDF parameters using the surface albedo as a reference shows an error of 0.01, which is mainly caused by a positive bias. However, this result lies within the error reported by Franch et al. (2014a).

AUTHOR CONTRIBUTIONS
BF leaded the work, applied the methodology and wrote the manuscript. EV contributed to the design of the work and supervised findings. SS prepared all the HLS data, worked on the preprocessing and revised the manuscript. J-CR supervised the findings and reviewed the manuscript. AS-A, extracted the field measurements. JV-N worked on the manuscript preparation. JM supervised the work.