A Method for Landsat and Sentinel 2 (HLS) BRDF Normalization

Abstract

The Harmonized Landsat/Sentinel-2 (HLS) project aims to generate a seamless surface reflectance product by combining observations from USGS/NASA Landsat-8 and ESA Sentinel-2 remote sensing satellites. These satellites’ sampling characteristics provide nearly constant observation geometry and low illumination variation through the scene. However, the illumination variation throughout the year impacts the surface reflectance by producing higher values for low solar zenith angles and lower reflectance for large zenith angles. In this work, we present a model to derive the bidirectional reflectance distribution function (BRDF) normalization and apply it to the HLS product at 30 m spatial resolution. It is based on the BRDF parameters estimated from the MODerate Resolution Imaging Spectroradiometer (MODIS) surface reflectance product (M{O,Y}D09) at 1 km spatial resolution using the VJB method (Vermote et al., 2009). Unsupervised classification (segmentation) of HLS images is used to disaggregate the BRDF parameters to the HLS spatial resolution and to build a BRDF parameters database at HLS scale. We first test the proposed BRDF normalization for different solar zenith angles over two homogeneous sites, in particular one desert and one Peruvian Amazon forest. The proposed method reduces both the correlation with the solar zenith angle and the coefficient of variation (CV) of the reflectance time series in the red and near infrared bands to 4% in forest and keeps a low CV of 3% to 4% for the deserts. Additionally, we assess the impact of the view zenith angle (VZA) in an area of the Brazilian Amazon forest close to the equator, where impact of the angular variation is stronger because it occurs in the principal plane. The directional reflectance shows a strong dependency with the VZA. The current HLS BRDF correction reduces this dependency but still shows an under-correction, especially in the near infrared, while the proposed method shows no dependency with the view angles. We also evaluate the BRDF parameters using field surface albedo measurements as a reference over seven different sites of the US surface radiation budget observing network (SURFRAD) and five sites of the Australian OzFlux network.

1. Introduction

The Harmonized Landsat/Sentinel-2 (HLS) project [1] provides a surface reflectance product that combines observations from USGS/NASA’s Landsat 8 and ESA’s Sentinel-2 satellites at moderate spatial resolution (30 m). The main goal is to provide a unique dataset based on both satellites’ data to improve the revisit time to 3–5 days depending on the latitude [1]. Along with a common atmospheric correction algorithm [2], geometric resampling to 30 m spatial resolution and geographic registration [1], the product is also corrected for BRDF effects and band pass adjustment.

The purpose of the band pass adjustment is to correct discrepancies between data coming from different sensors due differences their Relative Spectral Response (RSR) functions. In the case of the near infrared (NIR) spectral band, the spectral difference between Sentinel-2 and Landsat-8 is ~0.03%, due to the almost identical spectral functions of the analogous bands. In the red band, however, these differences are of ~3% [3] and need to be corrected. To correct the differences in the red band, Villaescusa-Nadal [3] employed a Spectral Band Adjustment Factor (SBAF) exponential model dependent on the NDVI which reduced the discrepancy between the sensors by ~55%.

Regarding the BRDF correction, Landsat and Sentinel-2 sampling characteristics provide nearly constant observation geometry and low illumination variation within each scene. However, when a surface reflectance time series combining measurements from both sensors is created, variations due to differences in view geometry between them arise. In extreme cases, the differences in the view zenith angle for a ground target can reach 20° from adjacent orbits only a few days apart. Additionally, variations in the seasonal illumination also impact the surface reflectance value. Gao [4] concluded that for Landsat-like narrow swath sensors, the major BRDF effect arises from the day of year effect, and can cause variations of 0.04–0.06 reflectance compared to mid-summer observations. Such angular effects can be corrected using a Bidirectional Reflectance Distribution Function (BRDF) model. However, the narrow angular sampling of moderate resolution sensors such as Landsat and Sentinel 2 complicates the BRDF parameters retrieval.

The official HLS product (version 1.4) provides a BRDF normalized product based on the BRDF normalization proposed by Roy [5]. In this product, the view angle is set to nadir and the solar zenith angle is fixed through time but varies for each tile based on the latitude. 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 towards BRDF estimation at moderate resolution for correcting the BRDF effects or for generating an albedo product of medium spatial resolution sensors. Shuai [6] used the MODIS Bidirectional Reflectance Distribution Function (BRDF) MDC43 product at 500 m [7] to generate a Landsat surface albedo product that later was implemented back to the 80s to the “pre-MODIS era” [8]. Their method achieves a Root Mean Square Error (RMSE) generally lower than 0.03 (12%) when compared with in situ data. However, this method does not provide a BRDF product at Landsat scale. Flood [9] used constant BRDF parameters to adjust surface reflectance from Landsat (TM and ETM+) and SPOT 5 images to a standard angular configuration. Gao [4] proposed a method to correct BRDF effects on Landsat and the advanced wide field sensor (AWiFS) based on detecting homogeneous areas with USDA cropland data layer (CDL) map to build a BRDF Look-Up Map (LUM) based on MODIS BRDF data. Van Doninck [10] assessed the directional effects on Landsat TM/ETM+ imagery over Amazonian forests and evaluated five different normalization techniques, including MODIS-based and Landsat-based BRDF. They concluded that the two MODIS-based methods analyzed only removed half of the angular effect in the infrared bands while the best results were obtained parametrizing the BRDF model using pairs of Landsat images. Franch [11] presented an alternative algorithm to calculate the Landsat (TM and ETM+) 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 [12,13]. The validation of the results against field measurements showed an RMSE of 0.015 (7%). Franch [14] tested this method to correct the BRDF effects of the HLS data over the Amazon forest in Peru. Their results show a good performance of the model when applied to HLS data with an error of 0.01 in the surface albedo retrieval.

In this paper, we describe the adaptation of the method of Franch [11,14] for HLS BRDF normalization, apply the proposed method to the HLS product from 2013 to 2017 and validate its applicability to Landsat 8 and Sentinel 2 data, test it on three homogeneous sites to assess its capability to remove BRDF artifacts, and finally evaluate it through the comparison to ground albedo measurements over SURFRAD and OzFlux networks.

2. Materials and Methods

2.1. Materials

2.1.1. 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. The HLS product is also corrected for bandpass difference and BRDF effects. However, in this work, we considered uncorrected data to explore the feasibility of the proposed algorithm for BRDF normalization. We used all HLS data (version 1.4) available from January 2013 to December 2017.

2.1.2. MODIS Data

We processed MODIS surface reflectance Collection 6 data at 1 km spatial resolution (M{OY}D09GA, Vermote, 2015) over the same areas and period, to derive the BRDF parameters using the VJB method [12,13].

2.1.3. Homogeneous Sites

Three homogeneous sites were selected for this study. First, we studied the impact of the solar zenith angle (SZA) variation through the year over two homogeneous sites. The first site is located on a homogeneous area of dense tropical forest in the Tambopata region of the Peruvian Amazon (12.818S, 69.281W, HLS tile 19LDI) where an albedometer is retrieving continuous measurements [15]. In this work, we have not included the albedo validation of this site, where Franch [14] reported an error of 0.01. This area is located in the southwest Amazon moist forest ecoregion of dense forests with a canopy height between 30 to 50 m high (https://www.worldwildlife.org/ecoregions/nt0173). The second site is located on a homogeneous desert area in Yuma, Arizona (32.499N, 114.498W, HLS tile 11SQS). This area is a lower-elevation section of the Sonoran desert in the southwestern United States and northwest of Mexico. We also analyzed the impact of the view zenith angle (VZA) variation. Because in the tropics the Landsat scan direction is closely aligned with the solar principal plane and the angular effects can be as strong as 20% in reflectance [16], we selected a site located close to the equator on the Brazilian Amazon basin (2.144S, 58.999W, HLS tile 21MTT). This area is located in the Uatuma-Trombetas moist forests ecoregion (Figure 1). The forests are characterized by a dense vegetation with a high number of small-diameter (less than 200 mm) to medium-diameter (200 to 600 mm) stems and a canopy height varying from 20 to 30 m, with emergent trees to 40 m (https://www.worldwildlife.org/ecoregions/nt0173).

2.1.4. SURFRAD Data

The global downwelling and upwelling solar radiations are continuously recorded by skyward facing and downward facing Eppley model 8–48 pyranometers mounted on a meteorological tower.

Table 1. Geolocation and brief description of the US surface radiation budget observing network (SURFRAD) sites considered in this study.

Station Name	Network	Location Latitude, Longitude	Land cover Type	HLS Tile	Tower Height above Target
Desert Rock Station	SURFRAD	36.6232N, 116.0196W	Sparse vegetation	11SNA	10 m
Table Mountain	SURFRAD	40.1256N, 105.2378W	Sandy with exposed rocks, sparse grasses and shrubs	13TDE	10 m
Bondville	SURFRAD	40.0516N, 88.3733W	Agriculture	16TCK	10 m
Goodwin creek	SURFRAD	34.2547N, 89.8729W	Pasture grass and sparsely distributed deciduous trees	16SBD	10 m
Penn state university	SURFRAD	40.7203N, 77.9310W	Agriculture Research field	18TTL	10 m
Fort Peck	SURFRAD	48.3078N, 105.1017W	Sparse vegetation	13UDP	10 m
Sioux Falls	SURFRAD	43.7340N, 96.62328W	Prairie grasses	14TPP	10 m

specifies the height of each tower above the target.

2.1.5. OzFlux Data

TERN OzFlux is a national ecosystem research network set up to provide the Australian and global ecosystem modelling communities with nationally consistent observations of energy, carbon, and water exchange between the atmosphere and key Australian ecosystems. TERN OzFlux is part of an international network (FluxNet) of over 500 flux stations that is designed to provide continuous, long-term micrometeorological measurements to monitor the state of ecosystems globally [17].

From the 40 sites in this network, we selected 5 sites for validation purposes according to data availability during the period analyzed (2013 to 2017). In these sites, a pair of pyranometers from Kipp & Zonen CNR4 measure energy balance between incoming and reflected short-wave radiation at different heights (

Table 2. Geolocation and brief description of the OzFlux sites considered in this study.

Station Name	Network	Location Latitude, Longitude	Land cover Type	HLS Tile	Tower Height above Target
Calperum	OzFlux	34.0027S, 140.5877E	Sand dunes with trees and shrubs	54HVH	20 m
Cumberland Plain	OzFlux	33.6152S, 150.7236E	Dry sclerophyll forest	56HKH	29 m (~6 m above canopy)
Whroo	OzFlux	36.6732S, 145.0294E	Eucalyptus forest	55HCV	36 m (~10 m above canopy)
Wombat	OzFlux	37.4222S, 144.0944E	Eucalyptus forest	55HBU	30 m (~5 m above canopy)
Yanco	OzFlux	34.9893S, 146.2907E	Grassland	55HDB	2 m

). The surface albedo is estimated in the spectral range 300–2800 nm.

2.2. Methods

2.2.1. Current HLS BRDF Normalization

The HLS product (version 1.4) is normalized for BRDF effects using the method proposed by Roy [5]. This method uses fixed BRDF coefficients for each spectral band (i.e., a constant BRDF shape), derived from a large number of pixels in the MODIS 500 m BRDF product (MCD43) that are globally and temporally distributed.

$${\rho }^N\left({\theta }_s^{out},0,0;\lambda \right)=\rho \left({\theta }_s,{\theta }_v,\varphi ;\lambda \right)\ast c\left(\lambda \right)$$

(1)

$$c\left(\lambda \right)=\frac{k_0\left(\lambda \right)+k_1\left(\lambda \right)F_1\left({\theta }_s^{out},0,0\right)+k_2\left(\lambda \right)F_2\left({\theta }_s^{out},0,0\right)}{k_0\left(\lambda \right)+k_1\left(\lambda \right)F_1\left({\theta }_s,{\theta }_v,\varphi \right)+k_2\left(\lambda \right)F_2\left({\theta }_s,{\theta }_v,\varphi \right)}$$

(2)

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 [18], and F₂ is the geometric kernel, based on the Li-sparse model [19] but considering the reciprocal form given by Lucht [20] and corrected for the Hot-Spot process proposed by Maignan [21]. ${\theta }_s^{out}$ is the sun zenith angle of the normalized data, which is set to constant value depending on the tile central latitude [1]. The BRDF coefficients of the model (k₀, k₁, and k₂) are provided in [1,5].

2.2.2. Proposed BRDF Normalization Method

For our proposed method, we first subset a 9 × 9 km area centered over the albedometer tower locations. All images are then masked for clouds 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, G.H. and Hall, D.J. (1965). This classification is used to unmix the BRDF parameters retrieved from MODIS (at 1 km spatial resolution) to HLS spatial resolution (30 m). We run the ISODATA algorithm for maximum three iterations with the maximum number of clusters setup at 15.

Following the Vermote et al. (2009) notations, the surface reflectance (ρ) is written as:

$$\rho \left({\theta }_s,{\theta }_v,\varphi \right)=k_0\left[1+\frac{k_1}{k_0}{F}_1\left({\theta }_s,{\theta }_v,\varphi \right)+\frac{k_2}{k_0}{F}_2\left({\theta }_s,{\theta }_v,\varphi \right)\right]$$

(3)

where F₁ and F₂ are fixed functions of the observation geometry, and k₀, k₁, and k₂ are free parameters. From 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 [12,13]. This method assumes that the difference between the successive observations in time is mainly attributed to directional effects while the variation of k₀ is supposed small. Additionally, it assumes that R and V are represented by a linear function of the NDVI.

$$V=V_0+V_1 \mathrm{NDVI}$$

(4)

$$R=R_0+R_1 \mathrm{NDVI}$$

(5)

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:

$$\begin{array}{ll}{k}_0^{1k{m}_x}+k_1^{1k{m}_x}{F}_1& \left({\theta }_s,{\theta }_v,\varnothing \right)+k_2^{1k{m}_x}{F}_2\left({\theta }_s,{\theta }_v,\varnothing \right)\\ & =A_x\left(k_0^{C_1}+k_1^{C_1}{F}_1\left({\theta }_s,{\theta }_v,\varnothing \right)+k_2^{C_1}{F}_2\left({\theta }_s,{\theta }_v,\varnothing \right)\right)\\ & =B_x\left(k_0^{C_2}+k_1^{C_2}{F}_1\left({\theta }_s,{\theta }_v,\varnothing \right)+k_2^{C_2}{F}_2\left({\theta }_s,{\theta }_v,\varnothing \right)\right)+\cdots \\ & =N_x\left(k_0^{C_n}+k_1^{C_n}{F}_1\left({\theta }_s,{\theta }_v,\varnothing \right)+k_2^{C_n}{F}_2\left({\theta }_s,{\theta }_v,\varnothing \right)\right)\end{array}$$

(6)

where A_x, B_x …, N_x represent proportions of each the n classes within the x MODIS pixel. BRDF parameters for each class, which are unknowns in Equation (6), can be derived through matrix inversion. We only invert the 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^{HLS}$ as shown in Equation (7).

$$k_0^{HLS}={\rho }^{HLS}\left({\theta }_s,{\theta }_v,\varphi \right)-k_1^{HLS}{F}_1\left({\theta }_s,{\theta }_v,\varphi \right)-k_2^{HLS}{F}_2\left({\theta }_s,{\theta }_v,\varphi \right)$$

(7)

Figure 2 shows the red band BRDF parameters inverted using this method on an HLS image on June 20th of 2017 over the Yuma desert in Arizona (US).

The BRDF inversion is applied to each image. However, this may generate noise in the BRDF inversion caused mainly by clouds. The presence of clouds reduces the number of cloud-free pixels considered in the BRDF inversion, limiting the number of observations available for inversion (note that just images with a minimum of 30% cloud free pixels can be processed to have enough information to invert the model). Additionally, the inaccuracy of the cloud mask [22] may introduce noise in the BRDF unmixing of a given class.

Based on the individual BRDF unmixing of each image in the HLS database (starting from 2013 for Landsat 8, 2015 for Sentinel 2A and 2017 for Sentinel 2B), our goal is to reproduce the BRDF parameters V₀, V₁, R₀, and R₁ in Equations (4) and (5) at HLS spatial resolution. Next, we apply a linear regression for the V and R parameters versus the NDVI to derive the V₀, V₁, R₀, and R₁ parameters at HLS pixel level.

Finally, with the BRDF parameters at HLS resolution, we apply Equation (8) to derive the normalized surface reflectance (ρ^N) at 45 degrees of solar zenith angle and nadir observation (Vermote et al., 2009):

$${\rho }^N\left(45,0,0\right)=\rho \left({\theta }_s,{\theta }_v,\varphi \right)\ast \frac{1+V{F}_1\left(45,0,0\right)+R{F}_2\left(45,0,0\right)}{1+V{F}_1\left({\theta }_s,{\theta }_v,\varphi \right)+R{F}_2\left({\theta }_s,{\theta }_v,\varphi \right)}$$

(8)

2.2.3. Temporal Evaluation of Homogeneous Sites

We test the performance of the BRDF normalization by applying the proposed method to two homogeneous sites: forest and desert. To evaluate the correction and assuming that these surfaces remain stationary, we estimate the coefficient of variation (CV), for the reflectance time series. The CV is defined as the ratio of the standard deviation to the mean (Equation (9)).

$$\mathrm{CV}=\frac{stddev\left({\rho }_0,{\rho }_1,\dots {\rho }_n\right)}{average\left({\rho }_0,{\rho }_1,\dots {\rho }_n\right)}$$

(9)

where ρ_i is the surface reflectance of day i.

2.2.4. Spatial Evaluation of an Equatorial Region

We evaluate the spatial impact of the view zenith angle variation in a dense Amazonian forest close to the equator. To do this, we extract a transect consisting of seven circular homogeneous regions of 800 pixels next to the Urubu river, which cover different values of the VZA. The condition for homogeneity is that the standard deviation is lower than 5% in any band. From these we can estimate the average and standard deviation for each band. The results from using the directional reflectance, the HLS BRDF normalized, and the proposed BRDF normalization are then compared. A p-value for each case was calculated which provides a test of the null hypothesis that a coefficient for each regression is equal to zero, i.e., the regressor has no effect on the dependent variable (VZA in our case). A small p-value indicates that the null hypothesis can be rejected, meaning that there is a statistically significant relationship between the regressor and dependent variable. We used this statistic index to assess how significant the reflectance dependency with the VZA was.

2.2.5. Albedo Validation

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 [7] by integrals of the BRDF model through the black-sky albedo and the white-sky albedo [7]. Surface albedo is typically estimated in the spectral range 280–2800 nm and is therefore comparable with the broadband MODIS albedo (300–5000 nm). Therefore, we estimate the broadband albedo using the narrow to broadband equation proposed by Liang [23].

For validation purposes, assuming that every pyranometer has an effective field of view of 81°, the downward pyranometer covers an area equal to tan(81°)*height. According to this equation and considering that the tower height in every SURFRAD is 10 m, the area covered by the pyranometer is 63 m in radius (3 by 3 HLS pixels). Then, the weighted average of the nine HLS pixels is estimated considering the angle of the instrument and the center of each pixel. The same approach is applied to OzFlux sites to account for the field of view as function of the height of the instrument.

SURFRAD instruments provide radiation measurements every minute, and we just consider the observations that match the time of Landsat 8 and Sentinel 2 overpass time. On the other hand, since OzFlux instruments provide radiation measurements every 30 min, a linear interpolation is done to estimate the value at the Landsat 8 and Sentinel 2 overpass time. A careful analysis of cloud-free conditions around the tower was performed to discard observations that may include residual effects.

3. Results

3.1. Temporal Evaluation of Homogeneous Sites

First, focusing on the Peruvian Amazon forest, Figure 3 shows the directional (red), BRDF-normalized using Roy [5] (green), and BRDF-normalized using the proposed BRDF normalization (black) surface reflectance versus solar zenith angle 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. Results are shown for (a) the red band, (b) the near infrared (NIR) band, (c) the Normalized Difference Vegetation Index (NDVI), and (d) the SZA for each observation. Given the high frequency of cloud cover in this area, a total of 23 scenes were used in this analysis through the entire period 2013–2017.

All parameters show good consistency for both satellites’ products. The directional surface reflectance (red) shows a clear dependency with the SZA in the RED (r² = 0.56) but mostly in the NIR band (r² = 0.82), showing higher values for low solar angles and lower values for higher angles. This effect is not fully corrected with the current BRDF-normalization of the HLS product (green), which after the correction still shows a dependency on the SZA in both bands. However, it is corrected on the BRDF-normalized RED reflectance and minimized in the NIR reflectance by decreasing the correlation coefficient and providing a slope closer to zero when using the proposed algorithm (black), which normalize all observations to SZA = 45° and nadir observation. The NDVI barely shows any dependency on the SZA for any product. Figure 4 shows a subset of the NIR band directional (left) and BRDF-normalized using the proposed algorithm (right) surface reflectance of the HLS image that is mostly affected by the SZA effect according to Figure 3. The directional reflectance image shows higher values than the BRDF-normalized one over the scene.

Similar to Figure 3, Figure 5 shows the directional (red), BRDF-normalized using Roy [5] (green), and BRDF-normalized using the proposed BRDF normalization (black) surface reflectance versus the SZA in the desert site in Yuma (Arizona, US). The SZA of this site shows a greater variation during the year (from 20 to 60 degrees) compared to the previous site (from 20 to 45 degrees). Despite this, the directional reflectance barely shows any dependency with the SZA. However, the HLS BRDF-normalized data shows a high dependency (green) with higher values for high SZA and lower values for low SZA. This is not the case when applying the proposed algorithm (black), that shows more stable surface reflectance values.

Table 3. Coefficient of variation (CV) of the Peruvian Amazon pixel.

CV (%)	RED	NIR	NDVI
Directional reflectance	11.4	8.3	1.6
Current HLS BRDF normalization	9.3	6.0	1.6
Proposed BRDF normalization	7.6	4.5	1.6

and

Table 4. Coefficient of variation (CV) of a desert site in Arizona (USA).

CV (%)	RED	NIR	NDVI
Directional reflectance	2.8	2.3	6.8
Current HLS BRDF normalization	5.4	4.2	10.1
Proposed BRDF normalization	3.3	2.5	8.2

show the coefficient of variation of the two homogeneous sites considered. The forest shows a high CV in the red (11%) and the NIR (8%) directional surface reflectance. The CV reduction in both red and NIR bands is around 2% when applying the current HLS correction and 4% when applying the proposed algorithm. However, the NDVI CV remains low and constant for all surface reflectances. The desert’s CV, on the other hand, is much lower and does not show any significant change for any band or the NDVI when applying the proposed algorithm, but it does show an increase when applying the current BRDF correction of about 2% to 3%.

3.2. Spatial Evaluation of an Equatorial Region

Figure 6 shows the HLS NIR surface reflectance image. The directional reflectance (Figure 6a) shows higher values along the western part (compared to the eastern) where the VZA is larger (Figure 6c). Note that the solar azimuth angle image shows nearly constant values of 63 degrees. This means that the western part is observed in the backscattering direction, which explains the higher values observed. The current HLS correction (Figure 6b) reduces this angular effect. However, the west area still shows higher values. Finally, when applying the proposed normalization (Figure 6c) this illumination effect is minimized.

Figure 7 shows the transect surface reflectance versus the view zenith angle. Though there is variability through the different samples, we have added a trendline to show the angular dependency of each dataset. The directional reflectance shows the highest values and the largest slope with a p-value of 0.02, and it is followed by the current HLS BRDF normalization with a p-value of 0.08. The proposed method shows the lowest and more stable values with a slope closer to zero and a p-value of 0.31, meaning that there is no significant relationship between the corrected reflectance and the VZA. Additionally, the similar intercepts of the linear regressions show that the three datasets converge to a similar value at zero degrees observation. The analysis of the red band and NDVI (Appendix A) show similar conclusions but the directional effect is lower.

3.3. Surface Albedo Validation

We also evaluate the BRDF normalization by using the BRDF parameters to derive the surface albedo, which can be validated with field measurements. Figure 8 shows the validation of the surface albedo over (a) the SURFRAD network, (b) the OzFlux network, and (c) both networks. Each marker in the plot represents an individual Landsat 8 or Sentinel 2 overpass. OzFlux validation shows lower errors (0.015–0.016) than SURFRAD (0.020). Sentinel 2 albedo validation in SURFRAD shows an overestimation for albedos higher than 0.2. However, this is observed in a total of 5 dates in Desert Rock and may be caused by an inaccurate atmospheric correction during those dates. Nevertheless, the overestimation is included within the RMSE of the sites. Combining both networks, the surface albedo validation shows a RMSE of 0.019 for Landsat 8 and 0.018 for Sentinel 2. Figure 8d shows the result of applying the narrow to broadband equation to the directional reflectance. In this case, the errors increase to 0.028 for Landsat 8 and 0.030 for Sentinel 2. Note that an equivalent evaluation cannot be included for the current HLS BRDF normalization method because in Roy [5] it is explicitly stated that “this approach is not applicable for generation of Landsat surface albedo, which requires a full characterization of the surface BRDF”.

4. Discussion and Conclusions

Quality of atmospheric and BRDF correction is essential for any application of the HLS dataset. In this work, we present a new method to improve the BRDF normalization of the HLS data. It is based on the Franch [11] method, which was originally designed and evaluated to derive the surface albedo of Landsat TM and ETM+ data. In this manuscript we improve it, apply it to HLS data, explore its feasibility to correct the BRDF effects of the dataset, and validate it over two homogeneous sites and several land cover types of the SURFRAD and OzFlux networks.

The proposed improvement is based on reproducing the BRDF parameters’ NDVI dependency [12,13] at moderate spatial resolution (HLS). In the operational context, this means that we create an HLS pixel-based database of the BRDF parameters (V₀, V₁, R₀, and R₁) that is regularly updated to account for any land cover change. In this way, the effect of outliers and poor-quality pixels on the BRDF inversion is minimized, resulting in a more stable and robust BRDF model. We acknowledge that the assumption of a BRDF model being a function of the NDVI has limitations. For example, on sparse forests where the NDVI is not a good descriptor of canopy structure [24]. However, this simple model might be used for a rough correction of BRDF effects in reflectance time series. Although a full inversion of the BRDF model will give better results, some applications, such as real time processing, may want to trade accuracy for simplicity [25].

When applying the proposed method, the results show a decrease in the surface reflectance timeseries CV of 4% and a decrease of the correlation coefficient with the SZA for the forest site, and little to no dependency on the desert site (which is barely affected by angular effects). In contrast, the current HLS BRDF normalization algorithm under-corrects the BRDF effects on the forests site, and increases the CV on the desert site. It is known that vegetation in the Amazon has a phenological cycle as a result of small variations in temperature and precipitation throughout the year [26]. However, directional effects on surface reflectance can be larger than the differences in reflectance between floristically different vegetation types [27,28].

The evaluation of the spatial variability of the view zenith angle showed a clear dependency in the backscatter direction, leading to higher directional reflectance for larger view angles. The current BRDF method reduces this dependency, but still shows an under-correction of the signal in the dense forest area analyzed. This under correction was already reported by Roy [5]. Finally, the proposed method shows no dependency with the view zenith angle.

The evaluation against surface albedo field measurements show an RMSE error of 0.019 for Landsat 8 and 0.018 for Sentinel-2. The OzFlux network error is similar to Franch [11] (RMSE of 0.016), while SURFRAD sites show slightly larger errors. Note that Franch [11] evaluated the method over five SURFRAD sites from 2003 to 2006. The higher error in the SURFRAD network compared to OzFlux might be caused by higher errors in the atmospheric correction algorithm since the atmosphere over Australia exhibit lower aerosol content than over the US. Compared to previous studies that provide albedo validation, the proposed method shows lower errors than Shuai [6] (RMSE 0.024). These results are further evidence of the good performance of the surface reflectance product of Landsat 8 and Sentinel 2.