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

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₁ is the volume scattering kernel, based on the Rossthick function derived by Roujean et al. (1992) and F₂ 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₁ and F₂ are fixed functions of the observation geometry, but k₀, k₁, and k₂ are free parameters. Following this notation, we define V as k₁/k₀ and R for k₂/k₀. 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₁ and k₂ parameters since they describe the directional shape, while k₀ describes the reflectance magnitude. Considering the HLS surface reflectance and the classification image, we derive k${}_0^{\text{Landsat}/\text{Sentinel}2}$ image from this equation:

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).

Materials

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.

Results

Figure 1 shows the directional (red), BRDF-normalized using Roy et al. (2016) (green) and BRDF-normalized using Franch et al. (2014a) (black) surface reflectance for each day of the year (DOY) considered in the time series (2013–2017). Note that we have represented differently Landsat 8 (dots) and Sentinel 2 (triangle) data. We evaluate separately (a) the red band (band 4 for both Landsat 8 and Sentinel 2), (b) the near infrared (NIR) band (band 5 for Landsat 8 and band 8a for Sentinel 2), (c) the Normalized Difference Vegetation Index (NDVI) and (d) the Difference Vegetation Index (DVI=NIR-RED). The figure also displays in the secondary y-axis the solar zenith angle (SZA) for each observation.

Figure 1. Landsat 8 (dots) and Sentinel 2 (triangles) surface reflectance in the NIR (A) and RED (B), NDVI (C) and DVI (D) with no normalization (red), HLS BRDF normalization (green) and BRDF-normalized (black). The error bars displayed represent the uncertainty of the Landsat 8 surface reflectance product (Vermote et al., 2016), assuming the same error for Sentinel 2. The solar zenith angle (SZA) is displayed in blue in the secondary y axis.

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.

Table 1. Validation of the HLS albedo with field measurements.

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.

Conflict of Interest Statement

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Acknowledgments

We want to thank the Pontificia Universidad Católica del Peru for sharing the albedometer data. This work was supported by the NASA grant Support for the HLS (Harmonized Landsat-Sentinel-2) Project (no. NNX16AN88G).

References