Facilitating Advanced Sentinel-2 Analysis Through a Simplified Computation of Nadir BRDF Adjusted Reflectance
David Montero Miguel D. Mahecha César Aybar Clemens Mosig Sebastian Wieneke


Framework for Polarized Magnetic Neutron Scattering from Nanoparticle Assemblies with Vortex-Type Spin Textures
    Andreas Michels
    April 29, 2024
===============================================================================================================




The Sentinel-2 (S2) mission from the European Space Agency's Copernicus program provides essential data for Earth surface analysis. Its Level-2A products deliver high-to-medium resolution (10-60 m) surface reflectance (SR) data through the MultiSpectral Instrument (MSI).
To enhance the accuracy and comparability of SR data, adjustments simulating a nadir viewing perspective are essential. These corrections address the anisotropic nature of SR and the variability in sun and observation angles, ensuring consistent image comparisons over time and under different conditions. The c-factor method, a simple yet effective algorithm, adjusts observed S2 SR by using the MODIS BRDF model to achieve Nadir BRDF Adjusted Reflectance (NBAR).
Despite the straightforward application of the c-factor to individual images, a cohesive Python framework for its application across multiple S2 images and Earth System Data Cubes (ESDCs) from cloud-stored data has been lacking.
Here we introduce , a Python package crafted to convert S2 SR data to NBAR, supporting both individual images and ESDCs derived from cloud-stored data.
This package simplifies the conversion of S2 SR data to NBAR via a single function, organized into modules for efficient process management.
By facilitating NBAR conversion for both SAFE files and ESDCs from SpatioTemporal Asset Catalogs (STAC),  is developed as a flexible tool that can handle diverse data format requirements.
We anticipate that  will considerably contribute to the standardization and harmonization of S2 data, offering a robust solution for a diverse range of users across various applications. 
 is an open-source tool available at <https://github.com/ESDS-Leipzig/sen2nbar>.




§ INTRODUCTION


The Sentinel-2 (S2) mission from the Copernicus program of the European Space Agency (ESA), comprising the Sentinel-2A (S2A) and Sentinel-2B (S2B) satellites, is equipped with the MultiSpectral Instrument (MSI) sensor <cit.>. This sensor is designed to capture data across various spectral bands with high-to-medium spatial resolutions (Table <ref>). The Level-2A product from S2 provides surface reflectance (SR) data, which are indispensable for detailed Earth surface analysis. This product has been employed in diverse applications, examples include the estimation of carbon fluxes <cit.>, the investigation of vegetation dynamics using spectral indices <cit.>, and Land Use and Land Cover (LULC) products generation via Artificial Intelligence (AI) models <cit.>. The utility of S2 SR data is significantly enhanced through the incorporation into Earth System Data Cubes <cit.>, which offer an organized spatiotemporal framework, allowing for simplified multitemporal analyses.



While S2 SR data is widely used, its viewing angle (±10.3°) and field of view (20.6°) amplify Bidirectional Reflectance Distribution Function (BRDF) effects due to SR anisotropy <cit.>. To minimize BRDF effects, adjustments simulating a nadir viewing perspective are needed <cit.>. These corrections address the directional effects arising from the anisotropic nature of SR and the variability of sunlight and satellite viewing angles. Such adjustments are crucial for the consistent comparison of images taken at different times and under various sensor acquisition conditions. This is particularly important for the processing and analysis of analysis-ready ESDCs, which are increasingly utilized due to their organized spatiotemporal structure and the simplicity of generating them from cloud-stored data <cit.>.

The spectral parameters derived from the MODIS BRDF model facilitate the computation of directional reflectance across any specified sensor viewing and solar angles <cit.>. Using this framework, <cit.> introduced the c-factor method, a straightforward approach for adjusting Landsat SR data by applying the MODIS spectral BRDF model parameters. This adjustment yields Nadir BRDF Adjusted Reflectance (NBAR) by multiplying the observed Landsat SR with the ratio of reflectances predicted by the MODIS BRDF model for both the observed Landsat SR and a standard nadir view under fixed solar zenith conditions. <cit.> and <cit.> extended this methodology for multiple S2 spectral bands. Despite the simplicity of the c-factor method for individual S2 images, a unified Python framework for applying this conversion uniformly across multiple images, especially for ESDCs derived from SpatioTemporal Asset Catalogs (STAC), is missing so far.

This paper presents , an open-source Python package designed to facilitate converting S2 SR data to NBAR using the c-factor method.  is engineered to support both automatic NBAR conversions through a single function, accommodating SAFE files and ESDCs derived from STAC. The document is organized into the following sections: Section <ref> presents the  framework, elaborating on the modules that simplify the processing steps; Section <ref> illustrates the application of , providing practical examples; Section <ref> explores the limitations of  as well as its potential, with a particular emphasis on AI; and Section <ref> summarizes our conclusions.



§ THE  PYTHON PACKAGE



The  package, developed in Python, facilitates the conversion of S2 L2A Surface Reflectance values into NBAR values through a single function. To achieve this higher level, the package is structured into several modules, ensuring a methodical conversion process leveraging multidimensional arrays via  <cit.> and  <cit.>:





 §.§ The  module


This module stores the predefined MODIS BRDF spectral model parameters, namely f_iso(λ), f_vol(λ), and f_geo(λ) (as shown in Table <ref>), in addition to preserving the values corresponding to the spatial resolution of each spectral band. In accordance to the methodology proposed by <cit.>, our implementation encompasses all spectral bands except for the Aerosols (B01) and the narrow NIR (B8A) bands. Consequently, the NBAR values are not computed for B01 and B8A bands. Also note that for the S2 Red Edge bands, the BRDF spectral model parameters were determined through a process of linear interpolation between the red and NIR MODIS BRDF spectral model parameters <cit.>.



 §.§ The  module



This module is tasked with the extraction of two critical sets of data from the scene's metadata files ( and , Figure <ref>a), essential for the computation and harmonization of NBAR values: firstly, it extracts the sun and sensor viewing angles; secondly, it determines the processing baseline of the scene. The acquisition of this information is conducted directly from the metadata file when dealing with SAFE files[<https://sentinels.copernicus.eu/web/sentinel/user-guides/sentinel-2-msi/data-formats>]. Alternatively, when operating with ESDCs constructed via  or  <cit.>, the module retrieves the requisite data from the corresponding metadata STAC asset within a specific item in a STAC catalog using the  library. This dual-path approach ensures flexibility and adaptability in processing various data sources for NBAR computation and harmonization.



  §.§.§ Solar and sensor viewing angles:


The extraction process for solar and sensor viewing angles in degrees involves the aggregation of detector-specific data across each spectral band into a unified 4-dimensional array 𝐀∈ℝ^|Λ|×|Θ|× H × W. This array is structured with dimensions |Λ|=|B|+1, where B are the bands undergoing conversion (with an additional index allocated for solar information) and Λ is the set of bands added to the solar component; |Θ|=2 (the number of extracted angles, i.e. Θ={zenith, azimuth}); and H=W=23, reflecting the spatial resolution of this information (5 km). The assignment of spatial coordinate values within this grid uses the Upper Left (UL) X (ULX) and Y (ULY) coordinates, derived from the tile's geocoding information.



  §.§.§ Processing Baseline:


The Processing Baseline (PB) of the scene is extracted and converted into a floating-point numerical format. This PB version serves as a criterion for initiating a harmonization procedure for scenes possessing a PB version of 4.00 or greater.



 §.§ The  module



A kernel-based BRDF model integrates various scattering mechanisms as a linear combination of distinct scattering modes, termed kernels <cit.>. The Ross-Li BRDF model (Section <ref>), also known as the Ross-Thick/Li-Sparse Reciprocal BRDF model, which is the standard for MODIS BRDF <cit.>, is composed of three components: the isotropic scattering parameter, the radiative transfer-type volumetric scattering kernel (Volumetric Kernel K_vol), and geometric-optical surface scattering kernel (Geometric Kernel K_geo). The objective of this module is to calculate the volumetric and geometric kernels. The computation of these kernels adheres to the mathematical formulations presented by <cit.>. Here we use the solar zenith angle θ=𝐀_{sun, zenith}, the view zenith angle ϑ=𝐀_{B, zenith}, and the view-sun relative azimuth angle ϕ=𝐀_{sun, azimuth}-𝐀_{B, azimuth} (all angles expressed in radians). Note that θ∈ℝ^H× W (since it is a 2-dimensional array representing the solar zenith angle) while ϑ,ϕ∈ℝ^|B|× H× W (since they are 3-dimensional arrays representing the zenith and relative azimuth angles for all bands). This implies that K_vol(θ,ϑ,ϕ), K_geo(θ,ϑ,ϕ)∈ℝ^|B|× H× W:



  §.§.§ Volumetric Kernel:


The Volumetric Kernel K_vol in the Ross-Li BRDF model is represented by the RossThick kernel <cit.>, which is based on the radiative transfer theory presented by <cit.>:


    K_vol(θ,ϑ,ϕ) = (π/2 - ξ) cosξ + sinξ/cosθ + cosϑ - π/4


Where cosξ = cosθcosϑ + sinθsinϑcosϕ.



  §.§.§ Geometric Kernel:


The Geometric Kernel K_geo in the Ross-Li BRDF model is represented by the LiSparse kernel <cit.>, which is based on the geometric-optical mutual shadowing BRDF model by <cit.>:

    K_geo(θ,ϑ,ϕ) = O - θ' - ϑ' + 1/2 (1 + cosξ') θ' ϑ'

where

    O    = 1/π (t - sin t cos t) (θ' + ϑ') 
    cos t    = h/b√(D^2 + (tanθ' tanϑ' sinϕ)^2)/θ' + ϑ'
    
      D    = √(tan^2 θ' + tan^2 ϑ' - 2 tanθ' tanϑ' cosϕ)
    cosξ'    = cosθ' cosϑ' + sinθ' sinϑ' cosϕ
    θ'    = tan^-1( b/rtanθ) 
    ϑ'    = tan^-1( b/rtanϑ)

Here, the cos t term is constrained to [-1,1], h/b= 2, and b/r= 1 <cit.>.



 §.§ The  module



This module is dedicated to the computation of the Ross-Li BRDF model, leveraging the scattering kernels sourced from the  module and the BRDF spectral parameters obtained from the  module, in accordance with the mathematical framework outlined by <cit.>. Prior to computation, the BRDF spectral model parameters are transformed into  objects such that f_iso, f_vol, f_geo ∈ℝ^|B|. This enhances the efficiency and scalability of the BRDF model computation. Consequently, this computation of the Ross-Li BRDF model outputs values as a multidimensional array BRDF(θ,ϑ,ϕ)∈ℝ^|B|× H× W:


    BRDF(θ,ϑ,ϕ) = f_iso + f_volK_vol(θ,ϑ,ϕ) + f_geoK_geo(θ,ϑ,ϕ)




 §.§ The  module



This module is tasked with calculating the  c-factor, adhering to the methodology delineated by <cit.>. The c-factor is instrumental in modifying the reflectance values to align with any given viewing or solar geometry, as outlined by <cit.>. To standardize the reflectance measurements to a nadir viewing zenith angle, the value of the viewing zenith angle, ϑ, is predetermined to be 0. Note that c(θ,ϑ,ϕ)∈ℝ^|B|× H× W:


    c(θ,ϑ,ϕ) = BRDF(θ,0,ϕ)/BRDF(θ,ϑ,ϕ)


Furthermore, this module integrates information and functionalities from preceding modules to directly derive the c-factor from a metadata file, establishing a comprehensive function specifically tailored for SAFE files. Additionally, recognizing the accessibility of metadata files through STAC items, a separate function has been developed to compute the c-factor directly from STAC items. This latter function is notably equipped with the capability to reproject the c-factor to a designated CRS when required. This reprojection feature is essential for ensuring the alignment and consistency of ESDCs created from scenes that originate in varying CRS, thereby facilitating their integration and analysis within a unified spatial framework.





 §.§ The  module


This module is designed to compute NBAR values, tailored to the format of the input data, whether they are SAFE files or ESDCs constructed via  or :



  §.§.§ SAFE files:


For SAFE files (Figure <ref>a), the spatial resolution is band-specific, denoted as λ∈ B, leading to the necessity of computing NBAR values for each λ independently. Let ρ_λ∈ℝ^H_λ× W_λ represent the SR of band λ, with H_λ and W_λ indicating its spatial dimensions. Initially, metadata files are used to extract the PB ∈ℝ_+ and to compute the c-factor c(θ,ϑ,ϕ)∈ℝ^|B|× H× W. The PB is then applied to harmonize ρ_λ if required:


    ρ*_λ =
        ρ_λ - 1000,    if PB≥4, 
    ρ_λ,    otherwise


Here, ρ*_λ signifies the harmonized SR values. Subsequently, the NBAR values are computed. It's important to note that the c-factor is derived using the original spatial resolution of the angle information (5 km), which makes the computation faster compared to performing angle interpolation beforehand. Therefore, c(θ,ϑ_λ,ϕ_λ) undergoes bilinear interpolation, denoted by the function I such that I:ℝ^H× W→ℝ^H_λ× W_λ, to align with the original resolution of band λ. The NBAR values are then generated as follows:


    NBAR_λ = I(c(θ,ϑ_λ,ϕ_λ))×ρ*_λ


Resulting NBAR images are subsequently exported individually into the same SAFE file directory, available in either COG or GeoTIFF formats.



  §.§.§ Earth System Data Cubes:


For ESDCs constructed from STAC (Figure <ref>b), the spatial resolution is predefined, enabling the generation of NBAR values for the entire ESDC via operations on multidimensional arrays. We denoted ρ_T∈ℝ^|T|× |B|× H_r× W_r as the SR of the ESDC, where T represents the temporal dimension items, and H_r and W_r are its spatial dimensions at resolution r, alongside the bounding box (BBox) of the ESDC. Initially, the module accesses the STAC, retrieving metadata files and PB values for each item τ∈ T. For each τ, the c-factor c_τ(θ_τ,ϑ_τ,ϕ_τ) is calculated and subsequently concatenated along the temporal dimension such that c_T=[c_τ(θ_τ,ϑ_τ,ϕ_τ)]_τ∈ T, resulting in c_T∈ℝ^|T|× |B|× H× W. The PB values for each τ are also concatenated as PB_T∈ℝ^|T|. This array facilitates the harmonization of ρ_T:


    ρ*_T = [
        ρ_τ - 1000,    if PB_τ≥ 4, 
    ρ_τ,    otherwise]_τ∈ T


Here, ρ*_T represents the harmonized SR values of the ESDC. Following this, NBAR values are calculated. Given that the c-factor is derived from the original spatial resolution of angle information, it undergoes bilinear interpolation through function I such that I:ℝ^|T|× |B|× H× W→ℝ^|T|× |B|× H_r× W_r, aligning with the predetermined resolution and BBox of ρ*_T. This process enables NBAR computation for ESDCs with uniform spatial dimensions (H_r=W_r, as in cubes generated via ) and for those with non-uniform dimensions (H_r≠ W_r, as in cubes generated via ). Note that I only interpolates the spatial dimensions H and W. The NBAR values are then derived as follows:


    NBAR_T = I(c_T)×ρ*_T




§ SHOWCASE



In a case study designed to showcase the capabilities of the  package in handling multidimensional data, we utilized an ESDC of S2 L2A data generated for the Hainich Natural National Park area, proximate to the DE-Hai Eddy Covariance tower <cit.>, a deciduous broadleaf forest in Germany, using  <cit.>. The ESDC was crafted with a uniform resolution of 10 meters (r=10) across all bands, featuring an edge size of 128 pixels, which translates to H_r=W_r=128. This ESDC covers the period from 2021-01-01 to 2023-12-31, focusing on images with less than 10% cloud coverage for analysis, encompassing a total of 90 images (|T|=90). It includes all S2 bands with the exception of B01 and B8A. Data from both satellites (S2A and S2B) were incorporated.

Utilizing the aforementioned ESDC, NBAR values were computed, followed by the calculation of the difference:


    Δρ=NBAR_T-ρ*_T




This metric shows the variations in reflectance values when transitioning from SR to NBAR data, particularly within the spatiotemporal framework of the ESDC. Figure <ref> illustrates the Δρ_λ for the ESDCs. Spatially, while the reflectance scales according to the c-factor and shows minimal variance owing to the confined area of study, the temporal variations are more pronounced, underscoring the significance of NBAR for accurate multitemporal analysis. This effect is especially marked in the Red Edge-NIR region, where discrepancies between values are more substantial compared to other bands. However, it's important to note that the maximum and minimum Δρ values exhibit similar ranges across all bands, highlighting the impact of NBAR adjustments.

Reflectance changes, as highlighted previously, can greatly influence the derivation of various products. To illustrate this effect, we computed four widely used vegetation indices for both the NBAR_T and ρ*_T: the Normalized Difference Vegetation Index <cit.>, the Near-Infrared Reflectance of Vegetation <cit.>, the Kernel NDVI <cit.>, and the Inverted Red Edge Chlorophyll Index <cit.>. Additionally, the difference between the indices derived from NBAR values and those from SR values was calculated as: 


    Δψ=ψ_NBAR-ψ_ρ*


Here, ψ represents the index in question. Figure <ref> displays the Δψ for the ESDCs. It is noted that the impact on NDVI is minimal, but it considerably increases for the other indices, particularly for kNDVI and IRECI, where the largest absolute changes observed reach values of 0.051 and 0.087, respectively. This variance underlines that the influence of reflectance adjustments on derived indices is index-dependent and can be substantial enough to potentially skew analysis reliant on these derived features.






§ DISCUSSION



Here we present various factors that either constrain or enhance the capabilities of , focusing on data providers and scalability. As we discuss below, this approach enhances the possibilities for big data computations, workflow generation towards analysis-ready data cubes, and therefore eases the application of novel artificial intelligence (AI) methods.



 §.§ The Impact of Metadata Availability


 is fully operational with SAFE files, which are widely accessible as they are publicly available (including metadata files). When it comes to ESDCs,  is compatible with those generated from STAC, particularly from public collections. A prime example of this compatibility is with the S2 L2A collection from the Microsoft Planetary Computer STAC, which possesses the required metadata files for transforming SR to NBAR values. Conversely,  does not support the S2 L2A collection from the Element84 STAC due to the absence of necessary metadata files. Similarly, while tools like  facilitate the creation of ESDCs from Google Earth Engine (GEE), the lack of metadata files in this context precludes the calculation of NBAR values. Future versions of  may overcome these limitations by enabling the interconnection of STAC assets (leveraging data from one STAC provider while sourcing metadata from another). This approach, however, depends on the availability of identical scenes across different collections.



 §.§ Parallel Processing and Lazy Evaluation


 employs multidimensional arrays through  and , enabling the execution of various array operations. In the context of ESDCs, the outcome of NBAR processing is another multidimensional array mirroring the dimensions of the input array. This architecture facilitates the acceleration of computations by exploiting lazy evaluations and parallel processing when necessary. Specifically, when an ESDC is generated via  or , the cube initially remains unevaluated, and its actual creation occurs only upon data retrieval. This setup allows the NBAR processing to be integrated into the lazy evaluation chain, as  accommodates  arrays <cit.>, ensuring that the NBAR computation is conducted concurrently with cube creation. Moreover, if an ESDC is pre-existing, either stored in memory or on disk, it can be divided into chunks for processing.  is then capable of executing the NBAR transformation in a lazy manner on these chunks. This approach is particularly effective with inherently chunked data formats, such as zarr, which can be processed lazily and re-written to disk in chunks, enhancing efficiency.



 §.§ Refining Multidimensional Dataset Production


 can enhance the accuracy of S2 data and derived features, supporting a wide range of applications, from individual snapshot analyses to comprehensive multitemporal studies. This extends to datasets generated from S2 data. This improvement will be largely attributed to the streamlined generation of ESDCs via STAC coupled with the data accuracy refinement process offered by . This conversion process ensures interoperability across the temporal dimension for all kind of applications. As a result,  streamlines the processing of multidimensional datasets like FluxnetEO <cit.>, which originates from MODIS NBAR data designed for vegetation monitoring but is adaptable to S2 data, or BigEarthNet <cit.>, derived from S2 patches aimed at land cover classification benchmarks. Therefore, the  package emerges as a crucial tool in remote sensing, improving the accuracy and utility of data for a broad spectrum of applications.



 §.§ Artificial Intelligence


The ultimate impact of  is its potential to greatly improve AI models in remote sensing. With the introduction of Transformer architectures <cit.>, AI has seen the development of new models that surpass traditional architectures across a variety of tasks, including Convolutional Neural Networks (CNNs) and Long-Short Term Memory (LSTM) networks. In remote sensing, the adoption of Transformer models has notably advanced, especially with pretrained models such as SITS-Former <cit.>. These foundation models, capable of learning from unlabeled datasets, are instrumental for various specific applications. Notable examples include Prithvi-100M <cit.> and SpectraGPT <cit.>, which are trained on Harmonized Landsat Sentinel <cit.> data and BigEarthNet, respectively. As these models form the basis for diverse applications, the role of  in improving the quality of input datasets becomes invaluable. By ensuring that these models have access to clean and accurate data, particularly in a multitemporal context,  will greatly contribute to the improvement of AI models' performance in remote sensing applications.



§ CONCLUSIONS



In this paper, we introduced , an open-source Python package crafted to facilitate the computation of Nadir BRDF Adjusted Reflectance (NBAR) values for Sentinel-2 L2A data. By enabling the computation of NBAR values through a single function,  greatly simplifies the process, demonstrating its adaptability and compatibility with SAFE files as well as Earth System Data Cubes (ESDCs) derived from STAC through  or . We foresee  playing an instrumental role in enhancing the accuracy of remote sensing data analyses within Earth system research. Its utility is especially pronounced in spatio-temporal analysis and in strengthening advanced AI-driven tasks via refined Sentinel-2 datasets.



§ ACKNOWLEDGEMENTS


D.M., M.D.M and S.W. acknowledge support from the “Digital Forest” project, Ministry of Lower-Saxony for Science and Culture (MWK) via the program Niedersächsisches Vorab (ZN 3679); and the “RS4BEF” project via the iDiv's Flexpool program. D.M. and M.D.M acknowledge support from the “DeepESDL” project from the European Space Agency (ESA). C.A. acknowledges support from the National Council of Science, Technology, and Technological Innovation (CONCYTEC, Perú) through the “PROYECTOS DE INVESTIGACIÓN BÁSICA – 2023-01” program with contract number PE501083135-2023-PROCIENCIA. M.D.M. and C.M. acknowledge support by the German Aerospace Center, DLR (“ML4Earth”) and by the Federal Ministry of Education and Research of Germany and by Sächsische Staatsministerium für Wissenschaft, Kultur und Tourismus in the programme Center of Excellence for AI-research “Center for Scalable Data Analytics and Artificial Intelligence Dresden/Leipzig”, project identification number: ScaDS.AI. Arrays within all Figures were visualized utilizing  <cit.>. Figures <ref> and <ref> have been designed using images from Flaticon.com.

abbrvnat


