Sentinel 2 topographic correction

Flatten the topography

Add a geometry and use the code below or click here.

var dem = ee.Image("USGS/SRTMGL1_003")

var S2 = ee.ImageCollection("COPERNICUS/S2");
var img = ee.Image(S2.filterBounds(geometry).sort("CLOUD_COVERAGE_ASSESSMENT").first())
.select(['B2','B3','B4','B8','B11','B12'], ['BLUE', 'GREEN', 'RED', 'NIR', 'SWIR1', 'SWIR2']) ;
Map.addLayer(img,{min:0,max:3000,bands:"SWIR1,NIR,RED"},"original");
Map.addLayer(img,{min:0,max:3000,bands:"SWIR1,NIR,RED"},"false color");

img = topoCorr_IC_L4toL8(img);
img = topoCorr_SCSc_L4toL8(img);
function topoCorr_IC_L4toL8(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,
'cosS': cosS.select('slope')});
// 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 topoCorr_SCSc_L4toL8(img){
var img_plus_ic = img;
var mask1 = img_plus_ic.select('NIR').gt(-0.1);
var mask2 = img_plus_ic.select('slope').gte(5)
.and(img_plus_ic.select('IC').gte(0))
.and(img_plus_ic.select('NIR').gt(-0.1));
var img_plus_ic_mask2 = ee.Image(img_plus_ic.updateMask(mask2));

var bandList = ['BLUE', 'GREEN', 'RED', 'NIR', 'SWIR1', 'SWIR2']; // Specify Bands to topographically correct

function apply_SCSccorr(bandList){
var method = 'SCSc';
var out = img_plus_ic_mask2.select('IC', bandList).reduceRegion({
reducer: ee.Reducer.linearFit(), // Compute coefficients: a(slope), b(offset), c(b/a)
geometry: ee.Geometry(img.geometry().buffer(-5000)), // trim off the outer edges of the image for linear relationship
scale: 10,
maxPixels: 1000000000
});
var out_a = ee.Number(out.get('scale'));
var out_b = ee.Number(out.get('offset'));
var out_c = ee.Number(out.get('offset')).divide(ee.Number(out.get('scale')));

//apply the SCSc correction
var SCSc_output = img_plus_ic_mask2.expression("((image * (cosB * cosZ + cvalue)) / (ic + cvalue))", {
'image': img_plus_ic_mask2.select(bandList),
'ic': img_plus_ic_mask2.select('IC'),
'cosB': img_plus_ic_mask2.select('cosS'),
'cosZ': img_plus_ic_mask2.select('cosZ'),
'cvalue': out_c
});

return ee.Image(SCSc_output);
}

var img_SCSccorr = ee.Image(bandList.map(apply_SCSccorr)).addBands(img_plus_ic.select('IC'));
var bandList_IC = ee.List([bandList, 'IC']).flatten();
return img_SCSccorr.unmask(img_plus_ic.select(bandList_IC))
.addBands(mask1.rename('initMask'))
.addBands(mask2.rename('corrMask'));
}