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

See an example.

Method is adopted from the paper below.

Soenen, S. A., Peddle, D. R., & Coburn, C. A. (2005).SCS+ C: A modified sun-canopy-sensor topographic correction in forested terrain. IEEE Transactions on geoscience and remote sensing, 43(9), 2148-2159.


var l8 = ee.ImageCollection("LANDSAT/LC08/C01/T1_SR");
var img = ee.Image(l8.filterBounds(geometry).sort("CLOUD_COVER").first());
var region = img.geometry();

//////////////////////////////////////////////////////////
///////////Correcting topographic effect in///////////////
/////////Landsat 8 surface reflectance collection/////////
//////////////////////////////////////////////////////////

//Same equations and expressions applied for Landsat 5, but changing the bands
//to their equivalents.

var degree2radian = 0.01745;

//Extract solar zenith and calculate incidence angle (i)
//Load USGS/SRTMGL1_003 DEM
var terrain = ee.call('Terrain', ee.Image('USGS/SRTMGL1_003'));
//Extract slope in radians for each pixel in the image
var p = terrain.select(['slope']).multiply(degree2radian);
//Extract solar zenith angle from the image
var z = ee.Number(img.get('SOLAR_ZENITH_ANGLE')).multiply(degree2radian);
//Extract solar azimuth from the image
var az = ee.Image(ee.Number(img.get('SOLAR_AZIMUTH_ANGLE')).multiply(degree2radian));
//Extract aspect in radians for each pixel in the image
var o = terrain.select(['aspect']).multiply(degree2radian);
var cosao = (az.subtract(o)).cos();//cos(ϕa−ϕo)
//Calculate the cosine of the local solar incidence for every pixel in the image in radians (cosi=cosp*cosz+sinp*sinz*cos(ϕa−ϕo)
var cosi = img.expression(
'((cosp*cosz) + ((sinp*sinz)*(cosao)))',
{
'cosp': p.cos(),
'cosz': z.cos(),
'sinp': p.sin(),
'sinz': z.sin(),
'az' : az,
'o' : o,
'cosao': cosao
});

// Create the image to apply the linear regression.The first band
//is the cosi and the second band is the response variable, the reflectance (the bands).
//L (y) = a + b*cosi(x); a = intercept, b = slope

// Dependent: Reflectance
var y = img.select(['B4', 'B5', 'B2']).multiply(0.0001);
// Independent: (cosi)
var x = cosi;
// Intercept: a
var a = ee.Image(1).rename('a');
// create an image collection with the three variables by concatenating them
var reg_img = ee.Image.cat(a,x,y);
// specify the linear regression reducer
var lr_reducer = ee.Reducer.linearRegression({
numX: 2,
numY: 3
});

// fit the model
var fit = reg_img.reduceRegion({
reducer: lr_reducer,
geometry: region,
scale: 30,
maxPixels: 1e9
});
fit = fit.combine({"coefficients": ee.Array([[1],[1]])}, false);
print(fit)

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

//Calculate C parameter C= a/b

var C = int.divide(slo);
print(C)

var newimg = img.expression(
'((img * ((cosp*cosz) + C))/(cosi + C))',
{
'img': img,
'cosp': p.cos(),
'cosz': z.cos(),
'cosi': cosi,
'C': C
});

Map.addLayer(ee.Image(newimg),{ bands: 'B4,B3,B2',min: 0, max: 3000},'corrected');
Map.addLayer(ee.Image(img),{ bands: 'B4,B3,B2',min: 0, max: 3000},'original');

 

2 thoughts on “Terrain correction in GEE”

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