BRDF Normalisation Script - Sentinel 2

//VERSION=3

// TRUE COLOR
var bands = ["B04", "B03", "B02"]
var lum = 2.5
var sample_type = "AUTO"

// FALSE COLOR
// var bands = ["B08", "B04", "B03"]
// var lum = 2.5
// var sample_type = "AUTO"

// RAW VALUES
// var bands = ["B02", "B03", "B04", "B08", "B11", "B12"]
// var lum = 1
// var sample_type = "FLOAT32"


function setup() {
    return {
        input: bands.concat(["sunAzimuthAngles", "sunZenithAngles", "viewAzimuthMean", "viewZenithMean"]),
        output: {
            bands: bands.length,
            sampleType: sample_type
        },
        mosaicking: "TILE"
    };
}

function evaluatePixel(sample) {
    for (let i = 0; i < sample.length; i++) {
        // Only calculate the BRDF for pixels where the 
        // view geometry is available (on the edges of tiles
        // reflectance is available but not view geometry)
        if (sample[i].viewAzimuthMean > 1) {
            var available = sample[i]
            var saa = deg2rad(available.sunAzimuthAngles);
            var sza = deg2rad(available.sunZenithAngles);
            var vaa = deg2rad(available.viewAzimuthMean);
            var vza = deg2rad(available.viewZenithMean);

            let constant = build_constants(sza, vza, saa, vaa);
            let c_vza_zero = build_constants(sza, 0, saa, vaa);
            let kernels = {
                kgeo: calc_kgeo(constant),
                kvol: calc_kvol(constant),
                kgeo_vza_zero: calc_kgeo(c_vza_zero),
                kvol_vza_zero: calc_kvol(c_vza_zero)
            }

            return bands.map(band => lum * calc_nbar(
                available[band],
                f_values[band],
                kernels));
        }
    }
    return [0, 0, 0]
}

// Kernel Parameters (Roy et al. 2017, Table 1)
// [f_iso, f_geo, f_vol]
var f_values = {
    "B02": [0.0774, 0.0079, 0.0372],
    "B03": [0.1306, 0.0178, 0.0580],
    "B04": [0.1690, 0.0227, 0.0574],
    "B08": [0.3093, 0.0330, 0.1535],
    "B11": [0.3430, 0.0453, 0.1154],
    "B12": [0.2658, 0.0387, 0.0639]
};

function build_constants(sza, vza, saa, vaa) {
    // calculates constants from viewing geometry that are often needed in the
    // calculations and are expensive to calculate (i.e. tan)
    const phi = relative_azimuth(saa, vaa);
    let c = {
        cos: {
            sza: Math.cos(sza),
            vza: Math.cos(vza),
            phi: Math.cos(phi),
        },
        sin: {
            sza: Math.sin(sza),
            vza: Math.sin(vza),
            phi: Math.sin(phi),
        },
        sec: {
            sza: sec(sza),
            vza: sec(vza),
        },
        tan: {}
    };
    c.tan.sza = c.sin.sza / c.cos.sza;
    c.tan.vza = c.sin.vza / c.cos.vza;
    //c 43 Lucht
    c.cos_xi = c.cos.sza * c.cos.vza + c.sin.sza * c.sin.vza * c.cos.phi;
    return c
}

function deg2rad(x) {
    // Convert degrees to radians
    return x * Math.PI / 180;
}

function sec(x) {
    // Calculate the secant of a value
    return 1 / Math.cos(x);
}

function relative_azimuth(saa, vaa) {
    // Calculate relative azimuth angle
    // Angles in RAD !
    let phi = Math.abs(saa - vaa)
    let diff = 0
    if (phi > Math.PI) {
        diff = 2 * Math.PI - phi;
    } else {
        diff = phi;
    }
    return diff;
}

function calc_kgeo(c) {
    // Calculate the LiSparse kernel from Lucht et al. 2000
    // Angles in RAD !

    //c 42 Lucht
    let Dsq = Math.pow(c.tan.sza, 2) + Math.pow(c.tan.vza, 2) - 2 * c.tan.sza * c.tan.vza * c.cos.phi;
    let tantansin = c.tan.sza * c.tan.vza * c.sin.phi;
    let costtop = Math.sqrt(Dsq + Math.pow(tantansin, 2))

    //c 41 Lucht
    let cost = 2 * costtop / (c.sec.sza + c.sec.vza)
    let t = Math.acos(Math.min(1, cost));

    // c 40 Lucht
    let big_o = (1 / Math.PI) * (t - Math.sin(t) * Math.cos(t)) * (c.sec.sza + c.sec.vza);

    // Kgeo
    let kgeo = big_o - c.sec.sza - c.sec.vza + 1 / 2 * (1 + c.cos_xi) * c.sec.sza * c.sec.vza;

    return kgeo;
}

function calc_kvol(c) {
    //Calculate the RossThick kernel (k_vol) from Lucht et al. 2000 equation 38
    // Angles in RAD !

    // eq 44
    let xi = Math.acos(c.cos_xi);

    let kvol = ((Math.PI / 2 - xi) * c.cos_xi + Math.sin(xi)) / (c.cos.sza + c.cos.vza) - Math.PI / 4;

    return kvol;
}

function calc_rho_modis(kgeo, kvol, f) {
    // Eq. 6 in Roy et al 2017, Eq. 37 in Lucht et al 2000
    return f[0] + f[1] * kgeo + f[2] * kvol;
}

function calc_c_lambda(kernels, f) {
    // Part 2 of Eq. 5 in Roy et al 2017
    return calc_rho_modis(kernels.kgeo_vza_zero, kernels.kvol_vza_zero, f) / calc_rho_modis(kernels.kgeo, kernels.kvol, f);
}

function calc_nbar(r_s2, f, kernels) {
    //Part 1 of Eq. 5 in Roy et al 2017
    // r_s2: reflectance in band 
    // f: f values for band
    let c_lambda = calc_c_lambda(kernels, f);
    return c_lambda * r_s2;
}

Evaluate and Visualize

• Copernicus Browser
• EO Browser

General description of the script

Implementation of a Bidirectional Reflectance Distribution Function (BRDF) normalisation for Sentinel-2 scenes in a Sentinel Hub Evalscript.

The Evalscript proposed here is the “semi-empirical BRDF normalisation c-factor correction approach proposed by Luch et al. 2000 [1]. The formula is detailed in eq. 37 (using eq. 38 & 39). To avoid lookup tables and other complex modelling approaches, we use the fixed BRDF spectral model parameters derived from “the global year of highest quality snow-free MODIS BRDF product” conveniently defined for Sentinel-2 in Roy et al. 2017 [2] that were derived from the Landsat ones in Roy et al. 2016 [3].

The approach seems to perform generally quite well for most surfaces (except snow and ice) and high sun angles. Nevertherless, it relies on a certain number of approximations and may not always be suitable.

Usage

BRDF corrections are available for Sentinel 2 bands 2, 3, 4, 8, 11 and 12. In the beginning of the script there is the option to define the bands which the correction should be applied to. It is also possible to set a luminance value in case RGB images are supposed to be generated. The luminance value is a multiplication factor designed to brighten the image for a better visualisation. By default, a factor of 2.5 is applied which works for most surface types. For darker surfaces, such as forests, a higher value may be required.

There are 3 output options available, a True Color and False Color visualization and the option to return Raw Values. Changing what the script does is just a matter of commenting and uncommenting the output options that are available. The output can also be customized following the same schema of the three available output options. This can for example be done to define custom composites using other bands.

Description of representative images

Top image: True Colour RGB acquired with Sentinel-2 (L2A processing level) over the Australian desert, on the border between Queensland and South Australia. The image is composed of two acquisitions from different orbits from 2nd January 2022 and 4th January 2022.

Bottom image: Identical to the top image, with the BRDF normalisation applied.

References

[1] W. Lucht, C. B. Schaaf and A. H. Strahler, “An algorithm for the retrieval of albedo from space using semiempirical BRDF models,” in IEEE Transactions on Geoscience and Remote Sensing, vol. 38, no. 2, pp. 977-998, March 2000, doi: https://doi.org/10.1109/36.841980.

[2] David P. Roy, Jian Li, Hankui K. Zhang, Lin Yan, Haiyan Huang, Zhongbin Li, Examination of Sentinel-2A multi-spectral instrument (MSI) reflectance anisotropy and the suitability of a general method to normalize MSI reflectance to nadir BRDF adjusted reflectance, Remote Sensing of Environment, Volume 199, 2017, Pages 25-38, ISSN 0034-4257, https://doi.org/10.1016/j.rse.2017.06.019.

[3] D.P. Roy, H.K. Zhang, J. Ju, J.L. Gomez-Dans, P.E. Lewis, C.B. Schaaf, Q. Sun, J. Li, H. Huang, V. Kovalskyy, A general method to normalize Landsat reflectance data to nadir BRDF adjusted reflectance, Remote Sensing of Environment, Volume 176, 2016, Pages 255-271, ISSN 0034-4257, https://doi.org/10.1016/j.rse.2016.01.023.