Add a bidirectional reflectance distribution function for sentinel 2

The Bidirectional Reflectance Distribution Function (BRDF) describes the directional dependence of the reflected energy of a target as a function of illumination and viewing geometries. The BRDF depends on wavelength.

The BRDF is needed in remote sensing for the correction of view and illumination angle effects.

1. Create a new script called brdf_module with the code below and save in your repository

var PI = ee.Number(3.14159265359);
var MAX_SATELLITE_ZENITH = 7.5;
var MAX_DISTANCE = 1000000;
var UPPER_LEFT = 0;
var LOWER_LEFT = 1;
var LOWER_RIGHT = 2;
var UPPER_RIGHT = 3;

/* Creates a collection of mosaics with a given temporal interval.
 *
 * collection - the collection from which to make composites.
 * start - the date of the first composite (either a string or an ee.Date)
 * count - the number of composites to make
 * interval - The time between composites, in units of "units".
 * units - The units of step (day, week, month, year; see ee ee.Date.advance)
 */
exports.brdfS2 = function(collection) {

  collection = collection.map(applyBRDF);
  return collection;

  function applyBRDF(image){
    var date = image.date();
    var footprint = ee.List(image.geometry().bounds().bounds().coordinates().get(0));
    var angles =  getsunAngles(date, footprint);
    var sunAz = angles[0];
    var sunZen = angles[1];

    var viewAz = azimuth(footprint);
    var viewZen = zenith(footprint);

    var kval = _kvol(sunAz, sunZen, viewAz, viewZen);
    var kvol = kval[0];
    var kvol0 = kval[1];
    var result = _apply(image, kvol.multiply(PI), kvol0.multiply(PI));

    return result;}

  /* Get sunAngles from the map given the data.
  *
  * date:  ee.date object
  * footprint: geometry of the image
  */
  function getsunAngles(date, footprint){
    var jdp = date.getFraction('year');
    var seconds_in_hour = 3600;
    var  hourGMT = ee.Number(date.getRelative('second', 'day')).divide(seconds_in_hour);

    var latRad = ee.Image.pixelLonLat().select('latitude').multiply(PI.divide(180));
    var longDeg = ee.Image.pixelLonLat().select('longitude');

    // Julian day proportion in radians
    var jdpr = jdp.multiply(PI).multiply(2);

    var a = ee.List([0.000075, 0.001868, 0.032077, 0.014615, 0.040849]);
    var meanSolarTime = longDeg.divide(15.0).add(ee.Number(hourGMT));
    var localSolarDiff1 = value(a, 0)
            .add(value(a, 1).multiply(jdpr.cos()))
            .subtract(value(a, 2).multiply(jdpr.sin()))
            .subtract(value(a, 3).multiply(jdpr.multiply(2).cos()))
            .subtract(value(a, 4).multiply(jdpr.multiply(2).sin()));

    var localSolarDiff2 = localSolarDiff1.multiply(12 * 60);

    var localSolarDiff = localSolarDiff2.divide(PI);
    var trueSolarTime = meanSolarTime
            .add(localSolarDiff.divide(60))
            .subtract(12.0);

    // Hour as an angle;
    var ah = trueSolarTime.multiply(ee.Number(MAX_SATELLITE_ZENITH * 2).multiply(PI.divide(180))) ;
    var b = ee.List([0.006918, 0.399912, 0.070257, 0.006758, 0.000907, 0.002697, 0.001480]);
    var delta = value(b, 0)
          .subtract(value(b, 1).multiply(jdpr.cos()))
          .add(value(b, 2).multiply(jdpr.sin()))
          .subtract(value(b, 3).multiply(jdpr.multiply(2).cos()))
          .add(value(b, 4).multiply(jdpr.multiply(2).sin()))
          .subtract(value(b, 5).multiply(jdpr.multiply(3).cos()))
          .add(value(b, 6).multiply(jdpr.multiply(3).sin()));

    var cosSunZen = latRad.sin().multiply(delta.sin())
          .add(latRad.cos().multiply(ah.cos()).multiply(delta.cos()));
    var sunZen = cosSunZen.acos();

    // sun azimuth from south, turning west
    var sinSunAzSW = ah.sin().multiply(delta.cos()).divide(sunZen.sin());
    sinSunAzSW = sinSunAzSW.clamp(-1.0, 1.0);

    var cosSunAzSW = (latRad.cos().multiply(-1).multiply(delta.sin())
                    .add(latRad.sin().multiply(delta.cos()).multiply(ah.cos())))
                    .divide(sunZen.sin());
    var sunAzSW = sinSunAzSW.asin();

    sunAzSW = where(cosSunAzSW.lte(0), sunAzSW.multiply(-1).add(PI), sunAzSW);
    sunAzSW = where(cosSunAzSW.gt(0).and(sinSunAzSW.lte(0)), sunAzSW.add(PI.multiply(2)), sunAzSW);

    var sunAz = sunAzSW.add(PI);
     // # Keep within [0, 2pi] range
    sunAz = where(sunAz.gt(PI.multiply(2)), sunAz.subtract(PI.multiply(2)), sunAz);

    var footprint_polygon = ee.Geometry.Polygon(footprint);
    sunAz = sunAz.clip(footprint_polygon);
    sunAz = sunAz.rename(['sunAz']);
    sunZen = sunZen.clip(footprint_polygon).rename(['sunZen']);

    return [sunAz, sunZen];
  }

  /* Get azimuth.
  *
  *
  * footprint: geometry of the image
  */
  function azimuth(footprint){

    function x(point){return ee.Number(ee.List(point).get(0))}
    function  y(point){return ee.Number(ee.List(point).get(1))}

      var upperCenter = line_from_coords(footprint, UPPER_LEFT, UPPER_RIGHT).centroid().coordinates();
      var lowerCenter = line_from_coords(footprint, LOWER_LEFT, LOWER_RIGHT).centroid().coordinates();
      var slope = ((y(lowerCenter)).subtract(y(upperCenter))).divide((x(lowerCenter)).subtract(x(upperCenter)));
      var slopePerp = ee.Number(-1).divide(slope);
      var azimuthLeft = ee.Image(PI.divide(2).subtract((slopePerp).atan()));
      return azimuthLeft.rename(['viewAz']);
  }

  /* Get zenith.
  *
  *
  * footprint: geometry of the image
  */
  function zenith(footprint){
      var leftLine = line_from_coords(footprint, UPPER_LEFT, LOWER_LEFT);
      var rightLine = line_from_coords(footprint, UPPER_RIGHT, LOWER_RIGHT);
      var leftDistance = ee.FeatureCollection(leftLine).distance(MAX_DISTANCE);
      var rightDistance = ee.FeatureCollection(rightLine).distance(MAX_DISTANCE);
      var viewZenith = rightDistance.multiply(ee.Number(MAX_SATELLITE_ZENITH * 2))
          .divide(rightDistance.add(leftDistance))
          .subtract(ee.Number(MAX_SATELLITE_ZENITH))
          .clip(ee.Geometry.Polygon(footprint))
          .rename(['viewZen']);
    return viewZenith.multiply(PI.divide(180));
  }

  /* apply function to all bands
  *
  * http://www.mdpi.com/2072-4292/9/12/1325/htm#sec3dot2-remotesensing-09-01325
  * https://www.sciencedirect.com/science/article/pii/S0034425717302791
  *
  * image : the image to apply the function to
  * kvol:
  * kvol0
  *
  */
function _apply(image, kvol, kvol0){
      var f_iso = 0;
      var f_geo = 0;
      var f_vol = 0;
			var blue = _correct_band(image, 'blue', kvol, kvol0, f_iso=0.0774, f_geo=0.0079, f_vol=0.0372);
			var green = _correct_band(image, 'green', kvol, kvol0, f_iso=0.1306, f_geo=0.0178, f_vol=0.0580);
			var red = _correct_band(image, 'red', kvol, kvol0, f_iso=0.1690, f_geo=0.0227, f_vol=0.0574);
			var re1 = _correct_band(image, 're1', kvol, kvol0, f_iso=0.2085, f_geo=0.0256, f_vol=0.0845);
			var re2 = _correct_band(image, 're2', kvol, kvol0, f_iso=0.2316, f_geo=0.0273, f_vol=0.1003);
			var re3 = _correct_band(image, 're3', kvol, kvol0, f_iso=0.2599, f_geo=0.0294, f_vol=0.1197);
      var nir = _correct_band(image, 'nir', kvol, kvol0, f_iso=0.3093, f_geo=0.0330, f_vol=0.1535);
      var re4 = _correct_band(image, 're4', kvol, kvol0, f_iso=0.2907, f_geo=0.0410, f_vol=0.1611);
      var swir1 = _correct_band(image, 'swir1', kvol, kvol0, f_iso=0.3430, f_geo=0.0453, f_vol=0.1154);
      var swir2 = _correct_band(image, 'swir2', kvol, kvol0, f_iso=0.2658, f_geo=0.0387, f_vol=0.0639);
			return image.select([]).addBands([blue, green, red, nir,re1,re2,re3,nir,re4,swir1, swir2]);
}

  /* correct band function
  *
  *
  * image : the image to apply the function to
  * band_name
  * kvol
  * kvol0
  * f_iso
  * f_geo
  * f_vol
  *
  */
  function _correct_band(image, band_name, kvol, kvol0, f_iso, f_geo, f_vol){
			//"""fiso + fvol * kvol + fgeo * kgeo"""
			var iso = ee.Image(f_iso);
			var geo = ee.Image(f_geo);
			var vol = ee.Image(f_vol);
			var pred = vol.multiply(kvol).add(geo.multiply(kvol)).add(iso).rename(['pred']);
			var pred0 = vol.multiply(kvol0).add(geo.multiply(kvol0)).add(iso).rename(['pred0']);
			var cfac = pred0.divide(pred).rename(['cfac']);
			var corr = image.select(band_name).multiply(cfac).rename([band_name]);
			return corr;
  }

  /* calculate kvol and kvol0
  *
  * sunAZ
  * sunZen
  * viewAz
  * viewZen
  *
  */
  function _kvol(sunAz, sunZen, viewAz, viewZen){
			//"""Calculate kvol kernel.
			//From Lucht et al. 2000
			//Phase angle = cos(solar zenith) cos(view zenith) + sin(solar zenith) sin(view zenith) cos(relative azimuth)"""

			var relative_azimuth = sunAz.subtract(viewAz).rename(['relAz']);
			var pa1 = viewZen.cos().multiply(sunZen.cos());
			var pa2 = viewZen.sin().multiply(sunZen.sin()).multiply(relative_azimuth.cos());
			var phase_angle1 = pa1.add(pa2);
			var phase_angle = phase_angle1.acos();
			var p1 = ee.Image(PI.divide(2)).subtract(phase_angle);
			var p2 = p1.multiply(phase_angle1);
			var p3 = p2.add(phase_angle.sin());
			var p4 = sunZen.cos().add(viewZen.cos());
			var p5 = ee.Image(PI.divide(4));

			var kvol = p3.divide(p4).subtract(p5).rename(['kvol']);

			var viewZen0 = ee.Image(0);
			var pa10 = viewZen0.cos().multiply(sunZen.cos());
			var pa20 = viewZen0.sin().multiply(sunZen.sin()).multiply(relative_azimuth.cos());
			var phase_angle10 = pa10.add(pa20);
			var phase_angle0 = phase_angle10.acos();
			var p10 = ee.Image(PI.divide(2)).subtract(phase_angle0);
			var p20 = p10.multiply(phase_angle10);
			var p30 = p20.add(phase_angle0.sin());
			var p40 = sunZen.cos().add(viewZen0.cos());
			var p50 = ee.Image(PI.divide(4));

			var kvol0 = p30.divide(p40).subtract(p50).rename(['kvol0']);

			return [kvol, kvol0]}

  /* helper function
  *
  *
  *
  */  

  function line_from_coords(coordinates, fromIndex, toIndex){
      return ee.Geometry.LineString(ee.List([
        coordinates.get(fromIndex),
        coordinates.get(toIndex)]));
  }

  function where(condition, trueValue, falseValue){
        var trueMasked = trueValue.mask(condition);
        var falseMasked = falseValue.mask(invertMask(condition));
        return trueMasked.unmask(falseMasked);
    }

    function invertMask(mask){
      return mask.multiply(-1).add(1);
    }

    function value(list,index){
      return ee.Number(list.get(index));
    }

};

/* Creates a collection of mosaics with a given temporal interval.
 *
 * collection - the collection from which to make composites.
 * start - the date of the first composite (either a string or an ee.Date)
 * count - the number of composites to make
 * interval - The time between composites, in units of "units".
 * units - The units of step (day, week, month, year; see ee ee.Date.advance)
 */
exports.brdfL8 = function(image) {

  var date = image.date();
  var footprint = ee.List(image.geometry().bounds().bounds().coordinates().get(0));
  var angles =  getsunAngles(date, footprint);
  var sunAz = angles[0];
  var sunZen = angles[1];

  var viewAz = azimuth(footprint);
  var viewZen = zenith(footprint);

  var kval = _kvol(sunAz, sunZen, viewAz, viewZen);
  var kvol = kval[0];
  var kvol0 = kval[1];
  var result = _apply(image, kvol.multiply(PI), kvol0.multiply(PI));

  return result;

  /* Get sunAngles from the map given the data.
  *
  * date:  ee.date object
  * footprint: geometry of the image
  */
  function getsunAngles(date, footprint){
    var jdp = date.getFraction('year');
    var seconds_in_hour = 3600;
    var  hourGMT = ee.Number(date.getRelative('second', 'day')).divide(seconds_in_hour);

    var latRad = ee.Image.pixelLonLat().select('latitude').multiply(PI.divide(180));
    var longDeg = ee.Image.pixelLonLat().select('longitude');

    // Julian day proportion in radians
    var jdpr = jdp.multiply(PI).multiply(2);

    var a = ee.List([0.000075, 0.001868, 0.032077, 0.014615, 0.040849]);
    var meanSolarTime = longDeg.divide(15.0).add(ee.Number(hourGMT));
    var localSolarDiff1 = value(a, 0)
            .add(value(a, 1).multiply(jdpr.cos()))
            .subtract(value(a, 2).multiply(jdpr.sin()))
            .subtract(value(a, 3).multiply(jdpr.multiply(2).cos()))
            .subtract(value(a, 4).multiply(jdpr.multiply(2).sin()));

    var localSolarDiff2 = localSolarDiff1.multiply(12 * 60);

    var localSolarDiff = localSolarDiff2.divide(PI);
    var trueSolarTime = meanSolarTime
            .add(localSolarDiff.divide(60))
            .subtract(12.0);

    // Hour as an angle;
    var ah = trueSolarTime.multiply(ee.Number(MAX_SATELLITE_ZENITH * 2).multiply(PI.divide(180))) ;
    var b = ee.List([0.006918, 0.399912, 0.070257, 0.006758, 0.000907, 0.002697, 0.001480]);
    var delta = value(b, 0)
          .subtract(value(b, 1).multiply(jdpr.cos()))
          .add(value(b, 2).multiply(jdpr.sin()))
          .subtract(value(b, 3).multiply(jdpr.multiply(2).cos()))
          .add(value(b, 4).multiply(jdpr.multiply(2).sin()))
          .subtract(value(b, 5).multiply(jdpr.multiply(3).cos()))
          .add(value(b, 6).multiply(jdpr.multiply(3).sin()));

    var cosSunZen = latRad.sin().multiply(delta.sin())
          .add(latRad.cos().multiply(ah.cos()).multiply(delta.cos()));
    var sunZen = cosSunZen.acos();

    // sun azimuth from south, turning west
    var sinSunAzSW = ah.sin().multiply(delta.cos()).divide(sunZen.sin());
    sinSunAzSW = sinSunAzSW.clamp(-1.0, 1.0);

    var cosSunAzSW = (latRad.cos().multiply(-1).multiply(delta.sin())
                    .add(latRad.sin().multiply(delta.cos()).multiply(ah.cos())))
                    .divide(sunZen.sin());
    var sunAzSW = sinSunAzSW.asin();

    sunAzSW = where(cosSunAzSW.lte(0), sunAzSW.multiply(-1).add(PI), sunAzSW);
    sunAzSW = where(cosSunAzSW.gt(0).and(sinSunAzSW.lte(0)), sunAzSW.add(PI.multiply(2)), sunAzSW);

    var sunAz = sunAzSW.add(PI);
     // # Keep within [0, 2pi] range
    sunAz = where(sunAz.gt(PI.multiply(2)), sunAz.subtract(PI.multiply(2)), sunAz);

    var footprint_polygon = ee.Geometry.Polygon(footprint);
    sunAz = sunAz.clip(footprint_polygon);
    sunAz = sunAz.rename(['sunAz']);
    sunZen = sunZen.clip(footprint_polygon).rename(['sunZen']);

    return [sunAz, sunZen];
  }

  /* Get azimuth.
  *
  *
  * footprint: geometry of the image
  */
  function azimuth(footprint){

    function x(point){return ee.Number(ee.List(point).get(0))}
    function  y(point){return ee.Number(ee.List(point).get(1))}

      var upperCenter = line_from_coords(footprint, UPPER_LEFT, UPPER_RIGHT).centroid().coordinates();
      var lowerCenter = line_from_coords(footprint, LOWER_LEFT, LOWER_RIGHT).centroid().coordinates();
      var slope = ((y(lowerCenter)).subtract(y(upperCenter))).divide((x(lowerCenter)).subtract(x(upperCenter)));
      var slopePerp = ee.Number(-1).divide(slope);
      var azimuthLeft = ee.Image(PI.divide(2).subtract((slopePerp).atan()));
      return azimuthLeft.rename(['viewAz']);
  }

  /* Get zenith.
  *
  *
  * footprint: geometry of the image
  */
  function zenith(footprint){
      var leftLine = line_from_coords(footprint, UPPER_LEFT, LOWER_LEFT);
      var rightLine = line_from_coords(footprint, UPPER_RIGHT, LOWER_RIGHT);
      var leftDistance = ee.FeatureCollection(leftLine).distance(MAX_DISTANCE);
      var rightDistance = ee.FeatureCollection(rightLine).distance(MAX_DISTANCE);
      var viewZenith = rightDistance.multiply(ee.Number(MAX_SATELLITE_ZENITH * 2))
          .divide(rightDistance.add(leftDistance))
          .subtract(ee.Number(MAX_SATELLITE_ZENITH))
          .clip(ee.Geometry.Polygon(footprint))
          .rename(['viewZen']);
    return viewZenith.multiply(PI.divide(180));
  }

  /* apply function to all bands
  *
  * http://www.mdpi.com/2072-4292/9/12/1325/htm#sec3dot2-remotesensing-09-01325
  * https://www.sciencedirect.com/science/article/pii/S0034425717302791
  *
  * image : the image to apply the function to
  * kvol:
  * kvol0
  *
  */
function _apply(image, kvol, kvol0){
      var f_iso = 0;
      var f_geo = 0;
      var f_vol = 0;
			var blue = _correct_band(image, 'blue', kvol, kvol0, f_iso=0.0774, f_geo=0.0079, f_vol=0.0372);
			var green = _correct_band(image, 'green', kvol, kvol0, f_iso=0.1306, f_geo=0.0178, f_vol=0.0580);
			var red = _correct_band(image, 'red', kvol, kvol0, f_iso=0.1690, f_geo=0.0227, f_vol=0.0574);
      var nir = _correct_band(image, 'nir', kvol, kvol0, f_iso=0.3093, f_geo=0.0330, f_vol=0.1535);
      var swir1 = _correct_band(image, 'swir1', kvol, kvol0, f_iso=0.3430, f_geo=0.0453, f_vol=0.1154);
      var swir2 = _correct_band(image, 'swir2', kvol, kvol0, f_iso=0.2658, f_geo=0.0387, f_vol=0.0639);
			return image.select([]).addBands([blue, green, red, nir, swir1, swir2]);
}

  /* correct band function
  *
  *
  * image : the image to apply the function to
  * band_name
  * kvol
  * kvol0
  * f_iso
  * f_geo
  * f_vol
  *
  */
  function _correct_band(image, band_name, kvol, kvol0, f_iso, f_geo, f_vol){
			//"""fiso + fvol * kvol + fgeo * kgeo"""
			var iso = ee.Image(f_iso);
			var geo = ee.Image(f_geo);
			var vol = ee.Image(f_vol);
			var pred = vol.multiply(kvol).add(geo.multiply(kvol)).add(iso).rename(['pred']);
			var pred0 = vol.multiply(kvol0).add(geo.multiply(kvol0)).add(iso).rename(['pred0']);
			var cfac = pred0.divide(pred).rename(['cfac']);
			var corr = image.select(band_name).multiply(cfac).rename([band_name]);
			return corr;
  }

  /* calculate kvol and kvol0
  *
  * sunAZ
  * sunZen
  * viewAz
  * viewZen
  *
  */
  function _kvol(sunAz, sunZen, viewAz, viewZen){
			//"""Calculate kvol kernel.
			//From Lucht et al. 2000
			//Phase angle = cos(solar zenith) cos(view zenith) + sin(solar zenith) sin(view zenith) cos(relative azimuth)"""

			var relative_azimuth = sunAz.subtract(viewAz).rename(['relAz']);
			var pa1 = viewZen.cos().multiply(sunZen.cos());
			var pa2 = viewZen.sin().multiply(sunZen.sin()).multiply(relative_azimuth.cos());
			var phase_angle1 = pa1.add(pa2);
			var phase_angle = phase_angle1.acos();
			var p1 = ee.Image(PI.divide(2)).subtract(phase_angle);
			var p2 = p1.multiply(phase_angle1);
			var p3 = p2.add(phase_angle.sin());
			var p4 = sunZen.cos().add(viewZen.cos());
			var p5 = ee.Image(PI.divide(4));

			var kvol = p3.divide(p4).subtract(p5).rename(['kvol']);

			var viewZen0 = ee.Image(0);
			var pa10 = viewZen0.cos().multiply(sunZen.cos());
			var pa20 = viewZen0.sin().multiply(sunZen.sin()).multiply(relative_azimuth.cos());
			var phase_angle10 = pa10.add(pa20);
			var phase_angle0 = phase_angle10.acos();
			var p10 = ee.Image(PI.divide(2)).subtract(phase_angle0);
			var p20 = p10.multiply(phase_angle10);
			var p30 = p20.add(phase_angle0.sin());
			var p40 = sunZen.cos().add(viewZen0.cos());
			var p50 = ee.Image(PI.divide(4));

			var kvol0 = p30.divide(p40).subtract(p50).rename(['kvol0']);

			return [kvol, kvol0]}

  /* helper function
  *
  *
  *
  */  

  function line_from_coords(coordinates, fromIndex, toIndex){
      return ee.Geometry.LineString(ee.List([
        coordinates.get(fromIndex),
        coordinates.get(toIndex)]));
  }

  function where(condition, trueValue, falseValue){
        var trueMasked = trueValue.mask(condition);
        var falseMasked = falseValue.mask(invertMask(condition));
        return trueMasked.unmask(falseMasked);
    }

    function invertMask(mask){
      return mask.multiply(-1).add(1);
    }

    function value(list,index){
      return ee.Number(list.get(index));
    }

};

link

2. add the code below to your main script:

var brdf = require("users/apoortinga/landcoverS2:brdf_module");

print("applying sentinel brdf");
s2 =  brdf.brdfS2(s2);
Map.addLayer(ee.Image(s2.first()),{min:0,max:0.3,bands:"swir2,nir,red"},"after brdf");
print(ee.Image(s2.first()));

link

Advertisements