1

I'm using the Google Earth Engine JavaScript API to do a terrain correction on Sentinel-1 imagery, as in Andreas Vollrath's work. However, I want to use a higher-resolution Canadian digital elevation model, CDEM (see documentation here) instead of the SRTM DEM.

CDEM comes as an Image Collection rather than an image, so needs to be mosaicked into an image before being put into a larger terrain correction function, where I use function terrain.slope() to get the slope from the DEM. However, when I mosaic this dataset, terrain.slope() does not produce the slope output at the same high resolution of the CDEM dataset (though the values and data types look correct). However, if I do not use mosaic but call terrain.slope() on the mapped CDEM ImageCollection, it turns out correctly (but with no-data borders around individual images because of edge effects).

My questions:

  • can I make a mosaic of the CDEM that I can use in the terrain.slope() function, and have an appropriately scaled slope layer as output that can be used for terrain correction? - What part of the mosaicking process results in the lower resolution slope layer?

Note that I've tried mosaicking by using mosaic() and a median reducer, and both lead to this problem.

You can view code that illustrates my problem here: https://code.earthengine.google.com/331ddb0c3324256ebd4ec1ea0db4402c

The section that is problematic is:

var dataset2 = ee.ImageCollection('NRCan/CDEM');
var elevation2 = dataset2.mosaic()   //somehow this line causes problems for the slope layer***
Map.addLayer(elevation2, elevationVis, "CDEM")  //as a DEM it shows up fine on the map

//calculate slope in radians
var alpha_sRad3 = ee.Terrain.slope(elevation2).select('slope')
            .multiply(Math.PI/180).setDefaultProjection(proj);
Map.addLayer(alpha_sRad3, {}, 'slope CDEM') //lower resolution than the elevation2 layer - why??***

Whereas, if you do the slope function on each individual image in the CDEM dataset, it works:

var alpha_sRad4 = dataset2.map(function(image){
  return ee.Terrain.slope(image).select('slope')
            .multiply(Math.PI/180).setDefaultProjection(proj);
  })
Map.addLayer(alpha_sRad4, {}, 'slope CDEM collection')

Note that 1) this problem also exists for terrain.aspect(), 2) I moved this question here from StackOverflow, and 3) the terrain correction is also discussed here: Terrain correction (flattening) of Sentinel 1 images

0

2 Answers 2

2

On the Earth Engine documentation, it says that

"A few cases that require a fixed projection include:

  • Computing gradients (e.g. ee.Terrain.gradient or ee.Terrain.slope)
  • reduceResolution, for when you want to aggregate higher resolution pixels into lower resolution."

The problem is that ee.Terrain.slope() doesn't work on computed projections, and the solution is to use reproject()

2

I came across this same problem today and it took me awhile to realize I could just .map() the slope and aspect calculations (since they're being applied to an collection and not a single image) then mosaic (or mean/median) the resulting collection to get a single image.

// small test region
var ROI = ee.Geometry.Polygon( 
        [[[-138.2524779323139, 63.245633977207305],
          [-138.2524779323139, 62.229182841861466],
          [-135.2202513698139, 62.229182841861466],
          [-135.2202513698139, 63.245633977207305]]], null, false);

var dem = ee.ImageCollection('NRCan/CDEM') // note this is an IC not an Image

// now .map() the slope and aspect functions across the IC
var slope = dem.map(function(image){ return ee.Terrain.slope(image) }).mosaic().clip(ROI)
var aspect = dem.map(function(image){ return ee.Terrain.aspect(image) }).mosaic().clip(ROI)


Map.addLayer(dem, {min:0, max:2000}, "dem")
Map.addLayer(slope, {min:0, max:50}, "slope")
Map.addLayer(aspect, {min:0, max:360}, "aspect")

UPDATE: the above solution isn't very good, as it results in missing data gaps (a single pixel wide) along the borders of each image in the original collection. This seems to work better:

//get projection info relevant to my study
var L8proj = ee.ImageCollection("LANDSAT/LC08/C01/T1_SR")
        .filterBounds(ROI).first().projection()
    
// now mosaic IC, clip to ROI, and define projection (from ~23m to 30)
// .reproject() also works but is way slower for large areas
var dem = ee.ImageCollection('NRCan/CDEM').mosaic().clip(ROI).setDefaultProjection(L8proj)


var slope = ee.Terrain.slope(dem)

var aspect = ee.Terrain.aspect(dem)

Map.addLayer(dem, {min:0, max:2000}, "dem")
Map.addLayer(slope, {min:0, max:50}, "slope")
Map.addLayer(aspect, {min:0, max:360}, "aspect")

Not the answer you're looking for? Browse other questions tagged or ask your own question.