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');
Can this be applied to sentinel-2.
???…great work
LikeLike
I think the same method should work with sentinel 2
LikeLike
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!
LikeLike
That is correct, it is a paired band operation. I just updated the blog post with our updated code.
LikeLike
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).
LikeLike
That’s an interesting question and I wouldn’t know the answer. I always thought they are like synonyms.. Like your explanation though.
LikeLike
Thanks. I found the following paper useful for understanding the different types of corrections that can be applied to satellite data. The authors are precise/explicit in the definitions they use. I hope you find it useful!
https://esajournals.onlinelibrary.wiley.com/doi/full/10.1002/ecy.1730
LikeLike
That makes sense! Highly appreciated!
LikeLike
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? 🙂
LikeLike
you can cite the original papers. We just published this paper where we use the correction for landsat https://www.mdpi.com/2072-4292/11/7/831
LikeLike
Thank you for providing this great work – I will also be citing your work in our research!
LikeLiked by 1 person
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
LikeLike
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).
🙂
LikeLike
He used the SCS + C, it´s a improved method of SCS (2005).
In some place of the script they specify it:
// Function to apply the Sun-Canopy-Sensor + C (SCSc) correction method to each
LikeLike
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?
LikeLike
Thanks for this helpful post. Any idea why I get strange results when applied to Landsat 7? https://code.earthengine.google.com/b95004461fe1b7af661aaad126749d2d
It seems to work OK on Landsat 5.
LikeLike
Sorry, by ‘strange results’ I mean the spattering of white or light-colored pixels.
LikeLike
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.
LikeLike
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.
LikeLike
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!
LikeLike
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
LikeLike
Hi, I have the same error, how could you solve it?
LikeLike
I would guess there is no l8 image for the area you are doing the calculation for..
LikeLike
In my case, I added .copyProperties(img) in cloud mask function return output which required to run this terrainCorrection function
LikeLike
Hi, I have the same error, and I checked there is the l8 image for the area. How could you solve it?
LikeLike
Greetings everyone!!!
I am a newcomer to GEE.
How could I add a cloud cover and a median reducer to a terrain corrected collection?
Thank you!
LikeLike
Thanks for sharing this tutorial (topographic correcttion, right?)
Does someone know why the SCALE is 300?
LikeLike
to speed up the calculation. you could use 30.. but might hit the computational limit..
LikeLike
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?
LikeLike
Sorry to bother you again. Theoretically, your code is available for Landsat Collection 2.
But it doesn’t work. Could you help me to find the question?
https://code.earthengine.google.com/435ce9aaea1fa79fd8823d74cf3faf57
Thanks for your help!
LikeLike
The SOLAR_ZENITH_ANGLE and SOLAR_AZIMUTH_ANGLE changed for level 2.. also make sure you apply the new scaling
https://code.earthengine.google.com/13e4a529397f5bfc203052fc3c2c73b5
LikeLike
Thanks for your help! If I use your code to Sentinel 2, should some scales be modified?
LikeLike
I dont know if all assumption in the model apply to sentinel-2 as well. There is also a scale difference between the DEM (30m) and the sentinel-2 bands (10, 20 and 60m). You can find an implementation here https://mygeoblog.com/2018/07/27/sentinel-2-terrain-correction/
LikeLike
Thanks for your sharing. There are some problems when I use code to correct Sentinel 2 images, could you give me some advice? Thank you very much. https://code.earthengine.google.com/34ac55bfd623c0894e439516892d86b0
LikeLike
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.’
https://code.earthengine.google.co.in/de0b1916693599579f2d15103ef43da5
If anyone would be able to help in the slightest it would be greatly appreciated.
Thanks,
Ellie
LikeLike
Hi,
I’m also new in GEE and I wanted to export an image from Landsat 8 corected using this code but I keep getting an error related to incompatible data types: Float64 and Float32. What shoul I do?
Here is the code: https://code.earthengine.google.com/687582668c2bcceaf459e5c6c98dff5f
Thanks, Joaquín.
LikeLike
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));
I wonder what these two mask do
LikeLike