Modular landcover system step-5: topographic corrections

Add a bidirectional reflectance distribution function for sentinel 2

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

var scale = 300;

// get terrain layers
var dem = ee.Image("USGS/SRTMGL1_003");
var degree2radian = 0.01745;

exports.topoCorrection = function(collection) {

  collection =;
  collection =;

  //collection = correction.merge(notcorrection).sort("system:time_start");


  // Function to calculate illumination condition (IC). Function by Patrick Burns
  // ( and Matt Macander
  // (
  function illuminationCondition(img){

  // Extract image metadata about solar position
  var SZ_rad = ee.Image.constant(ee.Number(img.get('MEAN_SOLAR_ZENITH_ANGLE'))).multiply(3.14159265359).divide(180).clip(img.geometry().buffer(10000));
  var SA_rad = ee.Image.constant(ee.Number(img.get('MEAN_SOLAR_AZIMUTH_ANGLE')).multiply(3.14159265359).divide(180)).clip(img.geometry().buffer(10000));
  // Creat terrain layers
  var slp = ee.Terrain.slope(dem).clip(img.geometry().buffer(10000));
  var slp_rad = ee.Terrain.slope(dem).multiply(3.14159265359).divide(180).clip(img.geometry().buffer(10000));
  var asp_rad = ee.Terrain.aspect(dem).multiply(3.14159265359).divide(180).clip(img.geometry().buffer(10000));

  // Calculate the Illumination Condition (IC)
  // slope part of the illumination condition
  var cosZ = SZ_rad.cos();
  var cosS = slp_rad.cos();
  var slope_illumination = cosS.expression("cosZ * cosS",
                                          {'cosZ': cosZ,
  // aspect part of the illumination condition
  var sinZ = SZ_rad.sin();
  var sinS = slp_rad.sin();
  var cosAziDiff = (SA_rad.subtract(asp_rad)).cos();
  var aspect_illumination = sinZ.expression("sinZ * sinS * cosAziDiff",
                                           {'sinZ': sinZ,
                                            'sinS': sinS,
                                            'cosAziDiff': cosAziDiff});
  // full illumination condition (IC)
  var ic = slope_illumination.add(aspect_illumination);

  // Add IC to original image
  var img_plus_ic = ee.Image(img.addBands(ic.rename('IC')).addBands(cosZ.rename('cosZ')).addBands(cosS.rename('cosS')).addBands(slp.rename('slope')));
  return img_plus_ic;
  // Function to apply the Sun-Canopy-Sensor + C (SCSc) correction method to each
  // image. Function by Patrick Burns ( and Matt Macander
  // (
  function illuminationCorrection(img){
    var props = img.toDictionary();
    var st = img.get('system:time_start');

    var img_plus_ic = img;
    var mask1 ='nir').gt(-0.1);
    var mask2 ='slope').gte(5)
    var img_plus_ic_mask2 = ee.Image(img_plus_ic.updateMask(mask2));

    // Specify Bands to topographically correct
    var bandList = ['blue','green','red','nir','swir1','swir2'];
    var compositeBands = img.bandNames();
    var nonCorrectBands =;

    var geom = ee.Geometry(img.get('system:footprint')).bounds().buffer(10000);

    function apply_SCSccorr(band){
      var method = 'SCSc';
      var out =  ee.Image(1).addBands('IC', band))
                            .reduceRegion({reducer: ee.Reducer.linearRegression(2,1),
                                           geometry: ee.Geometry(img.geometry()),
                                           scale: scale,
                                           bestEffort :true,

        var fit = out.combine({"coefficients": ee.Array([[1],[1]])}, false);

				//Get the coefficients as a nested list,
				// ast it to an array, and get just the selected column
				var out_a = (ee.Array(fit.get('coefficients')).get([0,0]));
				var out_b = (ee.Array(fit.get('coefficients')).get([1,0]));
				var out_c = out_a.divide(out_b);

        // Apply the SCSc correction
        var SCSc_output = img_plus_ic_mask2.expression(
        "((image * (cosB * cosZ + cvalue)) / (ic + cvalue))", {
        'cvalue': out_c

      return SCSc_output;


    var img_SCSccorr = ee.Image('IC'));
    var bandList_IC = ee.List([bandList, 'IC']).flatten();
    img_SCSccorr = img_SCSccorr.unmask(;

    return img_SCSccorr.addBands(nonCorrectBands)



2. add the code below to your main script:

var topography = require("users/username/landcoverS2:topographic_module");

print("applying illumination correction");
s2 =  topography.topoCorrection(s2);
Map.addLayer(ee.Image(s2.first()),{min:0,max:0.3,bands:"swir2,nir,red"},"after terrain correction");


Leave a Reply

Fill in your details below or click an icon to log in: Logo

You are commenting using your account. Log Out /  Change )

Google photo

You are commenting using your Google account. Log Out /  Change )

Twitter picture

You are commenting using your Twitter account. Log Out /  Change )

Facebook photo

You are commenting using your Facebook account. Log Out /  Change )

Connecting to %s