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"},"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,
// 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 ='NIR').gt(-0.1);
var mask2 ='slope').gte(5)
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 ='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))", {
'cvalue': out_c

return ee.Image(SCSc_output);

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



  1. Hi, thanks for all your examples! I’m working with a map in Australia using Landsat level 1 Surface Reflectance data, I’m not sure if I need to perform the terrain correction or not. If yes, how to do preserve the properties of the images (e.g. MEAN_SOLAR_ZENITH_ANGLE) after using a reducer. In my case I’m taking the median of all the images in a 5 year period. However, the resulting image only contains the bands not the properties. Any help would be great.



  2. Hello, thanks for the code! Is there a literature reference for the algorithms for topographic correction that you used?

  3. Hello, Thank you so much. I would like to export the output image (to an asset), but I could not find the way to do it.

Leave a Reply