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.map(illuminationCondition);
  collection = collection.map(illuminationCorrection);

  return(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,
                                           '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 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 = 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));

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

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

    function apply_SCSccorr(band){
      var method = 'SCSc';
      var out = img_plus_ic_mask2.select('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 ){
       return img_plus_ic_mask2.select(band);
       }

  else{
      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))", {
        'image': img_plus_ic_mask2.select(band),
        '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 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();
    img_SCSccorr = img_SCSccorr.unmask(img_plus_ic.select(bandList_IC)).select(bandList);

    return img_SCSccorr.addBands(nonCorrectBands)
      .setMulti(props)
      .set('system:time_start',st);
  }

}  

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)
                   .filterDate("2017-01-01","2017-12-31")
                   .sort("CLOUD_COVER")
                   .select(inBands,outBands);

var img = ee.Image(collection.first());
print(img)
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');
Advertisements

8 comments

  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!

    Like

  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).

    Like

Leave a Reply

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

WordPress.com Logo

You are commenting using your WordPress.com 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