Topographic correction in GEE

Handle mountains and their shadows

Remote sensing image analysis can be negatively affected by the influence of terrain slope and aspect. Differences in terrain orientation create variation in signal values between pixels with similar land cover and biophysical-structural properties. The code below can be used to correct for terrain effects.

Note that I changed the title based on James Ball his comments below. There is a difference between terrain and topographic correction, as described in this paper:

“An important distinction should be made between topographic and terrain correction. Topographic correction is a radiometric process while terrain correction is geometric in nature. Although Landsat Level‐1 products are terrain corrected, this does not account for the same effects as a topographic correction. Terrain correction ensures each pixel is displayed as viewed from directly above regardless of topography or view angle, and, while important, does not account for the same effects as topographic correction.”

See an example.

var scale = 300;

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

var terrainCorrection = function(collection) {

  collection =;
  collection =;


  // 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('SOLAR_ZENITH_ANGLE'))).multiply(3.14159265359).divide(180).clip(img.geometry().buffer(10000));
  var SA_rad = ee.Image.constant(ee.Number(img.get('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 ='IC', band).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: 300,
      maxPixels: 1000000000

   if (out === null || out === undefined ){

      var out_a = ee.Number(out.get('scale'));
      var out_b = ee.Number(out.get('offset'));
      var out_c = out_b.divide(out_a);
      // 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)


var l8 = ee.ImageCollection("LANDSAT/LC08/C01/T1_SR");

var inBands = ee.List(['B2','B3','B4','B5','B6','B7'])
var outBands = ee.List(['blue','green','red','nir','swir1','swir2']);
var collection = l8.filterBounds(geometry)

var img = ee.Image(collection.first());
var collection = terrainCorrection(collection);
var newimg = ee.Image(collection.first())

Map.addLayer(ee.Image(newimg),{ bands: 'red,green,blue',min: 0, max: 3000},'corrected');
Map.addLayer(ee.Image(img),{ bands: 'red,green,blue',min: 0, max: 3000},'original');


  1. Hi,

    First that all, thank you very much for publishing this!

    I think there is a problem in the code when retrieving the parameters from the linear regressions. As it is now, the ratio of the intercept and the slope of the regression for the first band is used as the c parameter for all bands. However, as I understand the method, we should be using the ratio of the coefficients from band k to construct the c parameter for band k. That is, for each band k, we derive a parameter c_k that is equal to the ratio of the intercept and the slope of the regression of cosi on the reflectance of band k. I modified the last part of the code to correct this problem (this comes after fitting the linear regression):

    var coeff_array = ee.Array(fit.get(‘coefficients’));
    var int = ee.Array(coeff_array.toList().get(0))
    var slo = ee.Array(coeff_array.toList().get(1))
    var C = int.divide(slo)
    var Cimg = ee.Image.constant(C.toList())

    var newimg = img.expression(
    ‘((img * ((cosz) + C))/(cosi + C))’,
    ‘img’: img,
    ‘cosp’: p.cos(),
    ‘cosz’: z.cos(),
    ‘cosi’: cosi,
    ‘C’: Cimg

    Hope this is helpful!


  2. Is this not topographic correction rather than terrain correction? My understanding is that terrain correction ensures all pixels are viewed from above whereas topographic correction deals with influence of slope and aspect on illumination (as described here).


  3. Loving your blog! I ended up deciding to write a master thesis using Earth Engine to process Sentinel-2 time series. I am using this script with slight modifications to topographically correct my data. How do you want me to cite your work? 🙂


  4. Hello, I would like to modify the code to execute it in a polygon, but I try it and only gives me the correction of a point
    Any recommendation?
    thank you


  5. First of all, thank you so much for publishing this. I would like to ask that exactly kind of method did you use to create your scripts for the topographic correction in Landsat image. I just heard that there are currently eight methods available: “cosine”, “improvedcosine”, “minnaert”, “ccorrection” (first four from Riano et al. 2003), “minslope” (Minnaert with slope correction, also from Riano et al. 2003), “gamma” (from Richter et al. 2009), “SCS” (Gu and Gillespie 1998, Gao and Zhang 2009), “illumination” (uncorrected illumination).


      1. Hi, Cuevas. Will the moderator C, which is calculated from linear regression between radiance and illumination conditions, be biased if the image is containimated by cloud coverage?


  6. Hi, thanks for the script. I have one question. Why in the “mask” ‘nir’ should be greater than -0.1? What is the actual range of the values. Are they in TOC reflectance (i.e. 0 to 1). How the mask can prevent cloudy pixels to be sampled? Thanks.


    1. They are SR product, whose data range are 0 to 1. Question still unsovled. This criteria (nir > -0.1) is meaningless if SR product is used.


  7. Hi everyone!
    Great work!
    How do apply amedian reducer over the corrected collection?
    I am a newcomer to javaScript and GEE, and I woul love if you could give any advise!
    Thank you!


  8. Thank you for the code. I tried to combine the script with the cloud masking script by using .map() on the landsat image collection but everytime I did this, the topographic correction also failed with the report of “ImageCollection (Error)
    reduce.median: Error in map(ID=LC08_119064_20190109): Number.divide: Parameter ‘left’ is required.” do you have any explanation for this? Thank you


    1. In my case, I added .copyProperties(img) in cloud mask function return output which required to run this terrainCorrection function


  9. Thanks for sharing this tutorial.
    And I read your paper. It seems that the topographic correction was only used for Landsat 8 not used for Sentinel 2.
    Do not Sentinel 2 data need to be topographic corrected?


  10. Hello,
    I am verynnew to using GEE, so was wondering if anyone could help. I’ve tried to adapt this code to work on my cloud masked composite made from Landsat 8 collection 2 data, however I keep getting this error:

    ‘corrected: Layer error: reduce.median: Error in map(ID=LC08_009047_20160910):
    Image.constant: Parameter ‘value’ is required.’

    If anyone would be able to help in the slightest it would be greatly appreciated.



  11. var mask1 =‘nir’).gt(-0.1);
    var mask2 =‘slope’).gte(5)

    I wonder what these two mask do


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 )

Facebook photo

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

Connecting to %s