Remote sensing observations are influenced by the Sun-Sensor Geometry Effects, which are described by the BRDF. The surface reflectance measured by an orbital sensor is influenced by changes in the Sun-Sensor geometry due to differences in azimuth and zenith angles. Scattering models can be used to estimate Bidirectional Reflectance Factors (BRFs) as if they were observed in a specific view and Sun zenith angles. Nadir BRDF-Adjusted Reflectance (NBAR) products are surface reflectance estimated through a scattering model at nadir (0° view zenith) and a specified solar zenith angle.

Recently, global coefficients were used to correct BRDF effects through the c-factor approach for each spectral band (Roy et al., 2016). It has been evaluated for Landsat Time Series (4–8) (Roy et al., 2016; Qiu et al., 2019; Guan et al., 2020) and Sentinel-2/MSI (Mandanici and Bitelli, 2016; Roy et al., 2017; Kremezi and Karathanassi, 2019; Nguyen et al., 2020), and it is also being used to generate the Harmonized Landsat Sentinel (HLS) products (Claverie et al., 2018). The method generates adjusted reflectance calculated using the original reflectance and a c-factor (Eq. 1), which is calculated as the ratio between normalized and sensor BRFs (Eq. 2). ρ λ , θ Norm = c λ ⋅ ρ λ , θ sensor (1) c λ = B R F θ Norm B R F θ sensor (2)

The c-factor BRF is based on the Ross-Thick/Li-Sparse-Reciprocal (RTLSR) model (Eq. 3), which defines reflectance in terms of three basic scattering types: isotropic, volumetric, and geometric. The isotropic scattering describes Lambertian reflectance assumed to be an intrinsic property of an observed target, independent of Sun and observation geometry (Wagner et al., 2015).

RTLSR is calculated as a weighted sum of the isotropic, volumetric and geometric scattering terms, f _iso , f _vol , f _geo , respectively, and the volumetric and geometric scattering kernels, K _vol and K _geo , respectively, as can be seen in Eq. 3. f _iso , f _vol , f _geo are spectrally dependent BRDF model parameters (Schaaf et al., 2002; 2011), with fixed values for each spectral band (Roy et al., 2016) as can be seen in Table 1. B R F θ s , θ v , ϕ , λ = f iso λ + f vol λ K vol θ s , θ v , ϕ , λ + f geo λ K geo θ s , θ v , ϕ , λ (3)

The K _vol is calculated using the Ross-Thick from the radiative transfer model (Ross, 1981) through Eqs 4, 5. cos ξ = cos θ s cos θ v + sin θ s sin θ v cos ϕ (4) K vol θ s , θ v , ϕ = π 2 − ξ cos ξ + sin ξ cos θ s + cos θ v − π 4 (5)

The K _geo is calculated using the Li-Sparce-Reciprocal from surface scattering and geometric shadow casting theory (Li and Strahler, 1992) through Eqs 6–12. θ ′ = tan − 1 b r tan θ (6) cos ξ ′ = cos θ s ′ cos θ v ′ + sin θ s ′ sin θ v ′ cos ϕ (7) D = tan 2 θ s ′ + tan 2 θ v ′ − 2 t a n 2 θ s ′ tan θ v ′ cos ϕ (8) cos t = h b D 2 + tan θ s ′ tan θ v ′ sin ϕ 2 sec θ s ′ + sec θ v ′ (9) t = sec min max cos t , − 1 , 1 (10) O θ s , θ v , ϕ = 1 π t − sin t cos t sec θ s ′ + sec θ v ′ (11) K geo θ s , θ v , ϕ = O θ s , θ v , ϕ − sec θ s ′ − sec θ v ′ + 1 2 1 + cos ξ ′ sec θ v ′ sec θ s ′ (12)

Landsat-8 can revisit the same location in at least 16 days, with an overlap between adjacent Landsat-8 orbits that acquire images with 7–9 days of difference. This characteristic also happens for Sentinel-2/MSI images, with a 5-day revisit rate (considering both Sentinel-2A and Sentinel-2B), also presenting more frequent revisits for adjacent orbits overlap areas. This effect is beneficial to this study since the border between the images is more affected by view angle variation. However, there are a few regions where the Landsat overlapping orbit areas overlap with the sentinel-2 overlap areas. One of these areas occurs at Rio Grande do Sul in the south of Brazil, as can be seen in Figure 1, in the WRS2 path rows 223081 and 223082 and at MGRS Granule 22JBM.

For this study, we have selected images ranging from 2017 to 2020 from both Sentinel-2 and Landsat-8. Figure 2 illustrates the experiment processing diagram. It shows that Landsat-8/OLI Collection 2 images were acquired already processed as SR, Level 2, while Sentinel-2 images were processed to SR by Sen2cor or LaSRC through docker containers prepared for each processor. Once with the SR products, the angle bands were generated. Another docker containing the USGS Landsat Angles Creation Tool was used to generate the angle bands for the Landsat-8 data, while the s2angs Python package was used to generate the angle bands for the Sentinel-2 data. After that, the c-factor could be calculated and used to generate NBAR products. The sensor-harm Python package is capable of calculating and applying the c-factor on images from the sensors on board the Landsat satellites, while for Sentinel-2 data it also applies a bandpass coefficient, as performed by (Claverie et al., 2018). Then, it considers only pixels that were not masked as cloud, cloud shadows, or snow pixels, by each sensor’s cloud mask. Table 2 demonstrates the comparison performed to evaluate the differences between the used surface reflectance products.

These comparisons were performed using the Mean Absolute Difference ( Δ ρ ̄ λ ) described in Eq. 13 and Mean Relative Absolute Difference ( Δ ρ ̄ λ * ) described in Eq. (14). Δ ρ ̄ λ = ∑ i = 1 n ρ f o r w a r d , λ − ρ b a c k w a r d , λ n (13) Δ ρ ̄ λ * = ∑ i = 1 n 2 ρ f o r w a r d , λ − ρ b a c k w a r d , λ ρ f o r w a r d , λ + ρ b a c k w a r d , λ n × 100 (14)

Note that several works have already compared the Surface Reflectance products from Landsat and Sentinel products (Claverie et al., 2018) and obtained feasible results using images acquired on the same day. Since our objective consists of analyzing these images as a series, we adopted a close-in-time period to perform the comparison. Considering that adjacent Landsat images can be acquired with 7–9 days of difference, Landsat-8 x Landsat-8 pairs were acquired considering images with fewer than 10 days of difference. For Sentinel-2, considering that it has a higher revisit rate, this difference was adopted considering a maximum distance of 5 days, which was also adopted when comparing Landsat-8 images to Sentinel-2 images. This approach was used to obtain close-in-time images. The assumption here is that each pair of values being compared should be almost equal, ignoring differences from atmospheric conditions and changes on the imaging point, e.g., vegetation growth.

After the SR products comparison, all these products were processed to NBAR through the c-factor approach and all comparisons were repeated using this level of processing, as demonstrated in Table 3.

Since a few terminology differences can be found in the literature regarding the terms Replicability and Reproducibility, here we followed the Claerbout/Donoho/Peng terminology convention as presented by (Barba, 2018).

• Reproducible research: authors provide all the necessary data and the computer code to run the analysis again, re-creating the results.

To organize all the input data, code, processed data, and extracted information, the repository was organized as a Research Compendium (Gentleman and Temple Lang, 2007). A Research Compendium is a centralized unit that contains all objects of the research, which includes code, examples, and other research components that are controlled by a software system. This Research Compendium organizes the files in a structure that facilitates file sharing between different tools, maintaining a clear separation of which of these files are methods, input data (read-only), or output data (disposable) and specifies which software, considering its versioning, was originally used (Marwick et al., 2018). This Research Compedium is provided on the repository https://github.com/brazil-data-cube/compendium-harmonization.

Even though a software compilation can be in many times very straight forward, users may not be familiar with this process and can encounter difficulties performing it. In this context, Docker containers are a valid option to encapsulate software by providing virtualization, allowing code to be executed on environments with similar conditions to when it was initially created, facilitating reproducibility. Considering that this paper implemented the harmonization of Landsat and Sentinel-2 images through the c-factor approach, which has been described in Section 2.1 and requires other auxiliary methods to be executed, e.g., angle bands creation tool, containerization is a valid option to execute the same processing on different images and to allow the research workflow to be repeated by other users on the same or different data. Based on that, several docker images were created to process the data.

• Sen2cor Atmosphere Correction

To easily control the execution of all the processes, trying to avoid the necessity of human intervention to run the workflow, the Research Compendium provides two workflow tools to orchestrate the execution of the research, a Jupyter Notebook and a dagster application. Both tools can be executed using docker containers in order to achieve reproduction and replicability. These workflow environments handle all details from the processing of the input data to the final results.

The Jupyter Notebook is commonly used in Python and R communities. It allows users to run blocks of code and text notations, providing a controlled flux of the processes. The main advantage of the Jupyter Notebooks is that it allows the combination of code executions with text explanations, equations, and figures that can help in the comprehension of the performing task. The dagster application presents the workflow in a visual interface of the processes and can be used to rerun the processing, execute the tasks concurrently, and even distribute them among multiple processing machines.

Considering that image processing can generate big data volumes and that this study alone generated over a Terabyte for a single study area, within this research repository, a minimal example is provided to allow reproduction. This example consists of four images: a pair from Landsat-8 and a pair from Sentinel-2. The results of the workflow using the minimal examples are presented in Appendix 3.2, while the results using the complete set of data are described in Section 3. Also, the workflow input data can be changed, allowing replication to different study areas or dates.

Sen2cor is the default atmosphere processor used to correct Sentinel-2 Level-1C (L1C) images (Main-Knorn et al., 2017). For processing the images of this research, we created a docker image, which is accessible through the repository at https://github.com/brazil-data-cube/sen2cor-docker. This image facilitates the installation and execution of Sen2cor for atmospheric correction of Sentinel-2 L1C images.

Regarding the LaSRC processor (Vermote et al., 2018), we created a docker image, which can be accessed through the repository at https://github.com/brazil-data-cube/lasrc-docker. This image facilitates the installation and execution of LaSRC. Notably, this image is capable of correcting Landsat-8 Collection-1 Level-1 images and Sentinel-2 L1C images. However, since this study acquired Landsat-8 Collection-2 images, already Level-2, we exclusively employed the LaSRC image to process the Sentinel-2 images.

USGS provides a tool to generate angle bands using auxiliary metadata. The installation of this software requires it to be compiled. To generate View Azimuth, View Zenith, Solar Azimuth, and Solar Zenith angle bands for each Landsat scene, the software uses a ANG.txt file provided with each Landsat scene. We created a docker image that installed the Landsat Angle Creation Tool and given a ANG.txt generates the angle bands. The Docker image is available on the repository https://github.com/brazil-data-cube/landsat-angles-docker.

Concerning Sentinel-2 data, it is important to note that angle bands are not included in the provided data. Instead, several 23 × 23 matrices can be found within image metadata and can be used to estimate the sensor angle bands in a 5000-m spatial resolution (Pahlevan et al., 2017). have developed a code to estimate per-pixel angles. It is worth mentioning that their code was made using Python 2, a deprecated version of Python. To harness the existing code effectively, we undertook the task of adapting it, incorporating novel features as required. This development resulted in a new Python Package. It allows users to extract per-pixel angle bands from Sentinel-2 data, and this package is maintained in https://github.com/brazil-data-cube/s2-angs.

The harmonization procedure was implemented as a Python package in the GitHub repository https://github.com/brazil-data-cube/sensor-harm. This library supports harmonizing images from Landsat four to eight and Sentinel-2 images, even though in this study, Landsat four to seven were not used. It implements the c-factor method (Roy et al., 2016) and bandpass, as used in (Claverie et al., 2018), and, for doing so, it requires per-pixel angle bands, as described in 2.1.

Tables 4, 5, 6 present, using the same metrics as in (Roy et al., 2017), the mean absolute reflectance ( Δ ρ ̄ λ ) differences and mean absolute relative percentage differences ( Δ ρ ̄ λ * ) of SR and NBAR products comparison composed by only Landsat-8 images, only Sentinel-2 images, and series with both Landsat-8 and Sentinel-2 images, respectively. Figures 3, 4 summarize the mean absolute reflectance ( Δ ρ ̄ λ ) differences and mean absolute relative percentage differences ( Δ ρ ̄ λ * ) from these tables.

Although NBAR products also represent surface reflectance measurements (adjusted), here we will address by SR the products without the BRDF corrections while NBAR the products containing the c-factor correction. The results show that SR products are more heterogeneous than NBAR products, since for all sensor combinations, all bands of the NBAR products presented less difference than the Surface Reflectance products when considering mean absolute relative percentage. This corroborates the feasibility of the c-factor approach (Roy et al., 2016) to standardize Landsat-8 and Sentinel-2 measurements and the results obtained by (Roy et al., 2016; 2017; Claverie et al., 2018; Poortinga et al., 2019; Guan et al., 2020; Nguyen et al., 2020).

Analyzing the results on Tables 4, Table 5 we can observe that Landsat-8 is less affected by the BRDF effects than Sentinel-2 data, as, when considering the single sensor approach, Landsat-8 NBAR is more similar to Landsat-8 SR than Sentinel-2 NBAR is to Sentinel-2 SR. This was expected since the Sentinel-2 swath is wider than Landsat-8’s.

Table 6 illustrates a noteworthy observation: Landsat-8/OLI images, subjected to surface reflectance correction via LaSRC, exhibit a closer resemblance to Sentinel-2/MSI images that have undergone LaSRC correction, as opposed to those that have been processed with Sen2cor. This alignment in results was anticipated, a pattern also identified by (Marujo et al., 2021). Such consistency is to be expected given that LaSRC and Sen2cor employ distinct models for atmospheric correction. Consequently, these differences manifest in the surface reflectance (SR) products.

It is crucial to note, however, that this observation does not imply a superiority of one atmospheric correction method over the other. Rather, it indicates to use the same atmosphere correction algorithm to correct all images on the analyzed series to avoid differences in data due to differences in the atmosphere correction modeling.

It is expected that this effect occurs in other areas, as have also been observed by (Roy et al., 2016; 2017; Claverie et al., 2018; Poortinga et al., 2019; Guan et al., 2020; Nguyen et al., 2020). Based on that, to be able to quantify the SR and NBAR differences, it is recommended to replicate this study on different areas, a task that can be performed using the free and open-source code developed and provided alongside this research.

For this study area, we could observe that NBAR images were more similar to themselves than the SR products, reducing the mean absolute relative percentage differences. This can be noted in Figures 3, 4 as the dashed lines (SR) are mainly above the continuous lines (NBAR) for each experiment pair comparison (same color).

From Figure 3 we can also note that NIR and SWIR bands presented for Sentinel-2, higher differences when considering absolute values, which we can see from Figure 4 only occurred due to the high surface reflectance values measurements, probably caused by vegetated samples.

Still from Figures 3, 4 one can note that the pattern observed for each SR comparison is the same obtained for the NBAR comparison, though the NBAR pattern is smaller. This reinforces that the NBAR product is reducing the differences between the compared products even for the most affected pixels in a scene, the ones near the image border that are most affected from greater angle variation.

Considering that the reproduction of our experiment would be a time-consuming task, we prepared a minimum example that uses few images. More specifically two Landsat-8 and two Sentinel-2 images, as bellow.

• LC08_L2SP_222081_20171120_20200902_02_T1

Our entire methodology can be reproduced using the compendium repository. Even the input images are being stored in it. As a result, we obtained similar results in comparison to the experiment using the entire image dataset. The mean absolute reflectance ( Δ ρ ̄ λ ) differences and mean absolute relative percentage differences ( Δ ρ ̄ λ * ) of SR and NBAR products comparison were composed of only Landsat-8 images, only Sentinel-2 images, and series with both Landsat-8 and Sentinel-2 images, respectively. These results can be seen in Tables 7–9.