eARTh Engine: Turn cold pixels to a colorful Terrain

I am relatively new to Google Earth Engine and I have so much more to learn from this wonderful platform. During the process of reading its guides and examples and taking the first steps to writing small snippets on the Code Editor I came up with producing a colorful global terrain.

I was so excited with the final outcome, that I decided to share my process, as well as the entire code to all Earth Engine users. I expect some of you to benefit from this, while some others to identify mistakes and suggest improvements.

A brief summary of the process is that I produce a custom multidirectional hillshade from eight traditional hillshades, upon which I add a slope raster to produce a crisp shaded relief. On top of the shaded relief I add a hypsometric tint, water bodies, the seas and the oceans. The whole process is mostly artistic rather than scientific, hence the capital letters at the word eARTh at the title of this blog post!

At the following paragraphs I analyse step-by-step the overall process and explain what each part of the code does. You may skip all this chatter, if you wish, and scroll to the very end of the post, where I include a direct link to the source code.

The Data Source

First of all, I call the necessary datasets from the Earth Engine Data Catalog, as shown at the following list.

• ALOS DSM: Global 30m, which will be used for the production of the multi-directional shaded relief, as well as for the hypsometric tinting of the final terrain. This dataset will be also used for the production of the shoreline.
• ETOPO1: Global 1 Arc-Minute Elevation, which will be used for the bathymetry layer at the seas and the oceans.
• JRC Global Surface Water Mapping Layers, v1.2, which will be used for the inland water bodies, such as lakes and rivers.

I declare all three datasets as variables, which will be used later at the code.

/* DATA SOURCE */
/* ETOPO1: Global 1 Arc-Minute Elevation */
var NOAA = ee.Image('NOAA/NGDC/ETOPO1');

/* ALOS DSM: Global 30m */
var ALOS = ee.Image('JAXA/ALOS/AW3D30/V2_2').multiply(4);

/* JRC Global Surface Water Mapping Layers, v1.2 */
var GSWM = ee.Image('JRC/GSW1_0/GlobalSurfaceWater');

The Traditional Hillshade(s)

Among the numerous libraries available on Google Earth Engine, one may find the ee.Terrain library, which includes all the basic algorithms for terrain creation and analysis.

For shaded relief or hillshade, there is the classic ee.Terrain.hillshade, which computes a simple, traditional hillshade out of a digital elevation model (DEM) or a digital surface model (DSM).

As very clearly explained at the Hillshade function specification of ArcGIS Pro: […] “The traditional method calculates the hillshade using an illumination source from one direction using the altitude and azimuth properties to specify the sun’s position” […] And it further specifies:

• Azimuth is the sun’s relative position along the horizon (in degrees). This position is indicated by the angle of the sun measured clockwise from due north. An azimuth of 0 degrees indicates north, east is 90 degrees, south is 180 degrees, and west is 270 degrees.
• Altitude is the sun’s angle of elevation above the horizon and ranges from 0 to 90 degrees. A value of 0 degrees indicates that the sun is on the horizon, that is, on the same horizontal plane as the frame of reference. A value of 90 degrees indicates that the sun is directly overhead.

Since I want to produce a multi-directional hillshade, which will be calculated from more than one illumination source, I must somehow combine a number of traditional hillshades.

To do this, in Earth Engine I create eight different traditional hillshades and I declare them as variables. For input I assign the ALOS variable which calls the Global 30m ALOS DSM and for each one of them I assign an appropriate pair of azimuth/altitude values, as shown in Table 1.

Variable	Cardinal Direction	Azimuth	Altitude	Weight
N	North	0	36	0
NE	NorthEast	45	44	0
E	East	90	56	0
SE	SouthEast	135	68	0
S	South	180	80	0.1
SW	SouthWest	225	68	0.2
W	West	270	56	0.2
NW	NorthWest	315	44	0.5

Note that one may adjust the values of azimuth/altitude to meet local illumination conditions. For instance, in my blog post Shadowplay I use very different values from those illustrated on Table 1, adjusted to the local conditions of my area of interest.

Note also the Weight column. Later on at the scope I will combine all eight traditional hillshades to a single multi-directional hillshade. This column indicates the weights with which I multiply each traditional hillshade. I actually decrease their intensity, while I diminish some of them completely (where weight equals zero).

The following lines of code produce all eight traditinal hillshades and at the same time they multiply each one of them with their desired weight. I am doing this with the ee.Image.multiply function.

/* TERRAIN */
/* Traditional Hillshade (input,azimuth,altitude).multiply(weight) */
var N = ee.Terrain.hillshade(ALOS,0,36).multiply(0);
var NE = ee.Terrain.hillshade(ALOS,45,44).multiply(0);
var E = ee.Terrain.hillshade(ALOS,90,56).multiply(0);
var SE = ee.Terrain.hillshade(ALOS,135,68).multiply(0);
var S = ee.Terrain.hillshade(ALOS,180,80).multiply(0.1);
var SW = ee.Terrain.hillshade(ALOS,225,68).multiply(0.2);
var W = ee.Terrain.hillshade(ALOS,270,56).multiply(0.2);
var NW = ee.Terrain.hillshade(ALOS,315,44).multiply(0.5);

The Multidirectional Hillshade

To produce a multidirectional hillshade out of the eight traditional hillshades, all I have to do is add them with the ee.Image.add function. I am actually adding the weighted traditional hillshades, since I have passed the ee.Image.multiply function at the same variable that creates them. The code is as follows.

/* Multidirectional Hillshade */
var MULTI = N
    .add(NE)
    .add(E)
    .add(SE)
    .add(S)
    .add(SW)
    .add(W)
    .add(NW)
    .visualize({
    min:0,
    max:255,
    palette:[
        '#000000',
        '#ffffff'
        ],
    })
    .resample('bicubic')
    .updateMask(0.5);

Pay attention to the Image Visualization parameters. I stretch the whole image to its miximum – maximum (0 – 255) and I define a palette. The reason I have to define a palette is for making sure the multidirectional layer will mosaic later on at the process (keep reading).

Note also the .resample(‘bicubic’) argument. This is an instance of ee.Image.resample and I use it to make sure that the derived image will be as smooth as possible (I may also use other methods, like adding a Gaussian Blur for instance, but I want to save as much computing time and effort as possible from the overall process).

Finally, have a look at the .updateMask(0.5) argument. Here I use the ee.Image.updateMask to make the derived image 50% transparent. This is important because later in the scope I am about to mosaic it with another image and I want it to blend somehow properly.

The Slope Raster

Slope represents the rate of change of elevation for each digital elevation model (DEM) cell, as stated at the Slope function specification of ArcGIS Pro. As a cartographer, I very often use the slope raster to enhance the terrain, for example where it happens to have steep slopes, so it is a good idea to also apply it here at my Earth Engine scope.

So, I create a new variable to store the slope and I apply the ee.Terrain.slope function on the ALOS digital elevation model, as shown at the following lines of my code.

/* Slope */
var SLOPE = ee.Terrain.slope(ALOS)
    .multiply(2)
    .visualize({
        min:100,
        max:180,
        palette:[
            '#ffffff',
            '#000000'
            ]
        })
    .resample('bicubic')
    .updateMask(1);

I exaggerate a little bit the image by applying .multiply(2) and at the visualization parameters I narrow down the range of minimum-maximum values (100 – 180). This is an alternative way of filtering the image to not show the flat areas of the slope.

I define a palette with the exact same colors with those defined at the multidirectional hillshade, only this time I change their order. This inverts the slope image and will darken the steeper slopes.

The Shaded Relief

Now I will mosaic the multidirectional hillshade and the slope raster to a single image, which will produce the final Shaded Relief.

To do this, I declare a new ee.ImageCollection, I pass both the multidirectional hillshade and the slope raster as a list of images and I then I composite all the images in the collection, using the ee.ImageCollection.mosaic.

/* Shaded Relief */
var SHADED_RELIEF = ee.ImageCollection([
    SLOPE,
    MULTI
    ])
    .mosaic()
    .reduce(
      ee.Reducer.median()
      )
    .updateMask(1);

I then add the derived Shaded Relief image on the Map. My screen should look similar to Picture 1, depending on which part of the World I have centered my Map.

Map.addLayer(
  SHADED_RELIEF,{
      min:0,
      max:255,
      gamma:1
      },
    'Shaded Relief',
  true
);

The Hypsometric Tinting

The Shaded Relief itself is very impressive and I can spend hours staring at it. Still, a little color may not harm. One common cartographic technique is depicting elevation with color, where usually green hues are applied to lowland regions that gradually change to brown or yellow hues as elevation increases and end to white at mountain tops (even though there exist numerous other palette examples)!

This process is known as Hypsometric Tinting. The word “Hypsometric” is actually Greek! “Hypso” (ύψος) in Greek means altitude and “metric” (μέτρο) is a measurement or the process of making a measurement. So hypsometric (υψόμετρο) means measuring the altitude!

So let’s do this in Earth Engine. I pass once again the ALOS digital elevation model to a new variable (var ELEVATION) and at the visualization parameters I define which band will be used (read the dataset’s bands specification).

/* Elevation */
var ELEVATION = ALOS
    .visualize({
        bands:['AVE_DSM'],
        min:0,
        max:12500,
        palette:[
            '#386641',
            '#6a994e',
            '#a7c957',
            '#fdf7d6',
            '#ffffff'
            ]
        })
    .resample('bicubic')
    .updateMask(0.4);

I then define a color scheme for the palette. I selected a classic color scheme, even though one may select far more creative, crazy combinations. My colors will produce the following linear gradient (Gradient 1), which represents elevation from the sea level to the highest mountain tops.

At the minimum and maximum options, of the visualization parameters, I select the range of values that the colors will be applied. I would expect the maximum value to be the peak of Everest Mountain (8848m) but for some reason the values of that specific dataset exceed this limit (I don’t know why). So I have to adjust the maximum value to the mean altitude of my area of interest, every time I hit Run on my code.

Note also that I apply an image mask on the derived image. This will produce an image with 40% transparency.

Then, I add the Elevation layer on the Map. My screen should look similar to Picture 2, depending on which part of the World I have centered my Map.

Map.addLayer(
  ELEVATION,{
      // visualization
      },
    'Elevation',
  true
);

Surface Water

For the inland water bodies, such as lakes and rivers, I will use the JRC Global Surface Water Mapping Layers, v1.2 dataset from the catalogue. So, I declare a new variable to hold the water dataset (the variable declared at the beginning of the code) and at the visualization parameters I define which specific band should be consumed.

/* WATER */
/* Surface Water */
var SURFACE_WATER = GSWM
    .visualize({
        bands:['occurrence'],
        min:0,
        max:100,
        palette:[
            '#B9E9E7'
            ]
        })
    .resample('bicubic');

I give a yellowish/greenish blue color at the palette parameter and then I add the layer on the Map. My screen should look similar to Picture 3, depending on which part of the World I have centered my Map.

Map.addLayer(
  SURFACE_WATER,{
      // visualization
      },
    'Surface Water',
  true
);

The Sea

The Sea basically represents level zero. Everything above zero is land and everything below zero is undersea. To add the Sea in Earth Engine all I have to do is filter the ALOS digital surface model to only show values less than or equal to zero. This can be achieved with a mask .updateMask(ALOS.lte(0)).

For a more harmonious result, I give to the Sea the same blue color as the one I gave to the water bodies earlier in the scope. I just want a solid color (not a gradient) so at the palette option of the visualization parameters I assign only that color.

/* Sea */
var SEA = ALOS
    .updateMask(ALOS.lte(0))
    .visualize({
        bands:['AVE_DSM'],
        min:0,
        max:0,
        palette:[
            'B9E9E7'
            ]
        })
    .resample('bicubic');

I then add the Sea layer on the Map. My screen should look similar to Picture 4, depending on which part of the World I have centered my Map.

Map.addLayer(
  SEA,{
      // visualization
      },
    'Sea',
  true
);

Bathymetry

The ALOS DSM does not include information underseas. It does not show the depth of the ocean bottom. So I need somehow to fill its no data areas and at the same time to add on the Map the missing information of the bathymetry.

The word “Bathymetry” is also Greek! “Bathy” (βάθος) means depth in Greek language and “metry” (μέτρο) is a measurement or the process of making a measurement. So Bathymetry (Βαθυμετρία) means measuring the depth!

For Bathymetry I will use the ETOPO1 GDEM which includes ocean bathymetry. Of course, I need to mask it to only show values less than or equal to zero. Because the resolution of that DEM is lower than that of the ALOS, I apply the mask to actually show values less than or equal to -10 to ensure that the large pixels of ETOPO1 will not cover parts of the shoreline.

/* Bathymetry */
var BATHYMETRY = NOAA
    .updateMask(NOAA.lte(-10))
    .visualize({
        bands:['bedrock'],
        min:-5000,
        max:0,
        palette:[
            '#8ECCCB',
            '#ABE0DF',
            'B9E9E7'
            ]
        })
    .resample('bicubic');

Since this layer is intended to show bathymetry, I have to assign a palette that darkens as the depth increases. I also want this palette to include darker tones of the same blue hue that I assigned earlier to the Sea and the Water Bodies layers.

To do this, I generate a new palette at the awesome Coolors web app, which starts with that hue and gradually darkens. I then assign the colors at the palette option of the visualization parameters of the Bathymetry layer. The produced color scheme is illustrated at the following gradient (Gradient 2).

I then add the Bathymetry layer on the Map. My screen should look similar to Picture 5, depending on which part of the World I have centered my Map.

Map.addLayer(
  BATHYMETRY,{
      // visualization
      },
    'Bathymetry',
  false
);

Note that for the Bathymetry layer I have declared to be false when added to the Map. This means that the layer is being loaded when I hit Run, but it is hidden and I have to manually make it visible from the Layers manager. I do this for speed purposes and because that layer is suitable only for smaller scales.

Center the Map

The last step is to center the Map. I need to define the coordinates (longitude and latitude) of the center of the Map, as well as the zoom level of the initial extent.

Map.setCenter(
/* VIKOS - AOOS GORGE */
      20.7619,
      39.9186,
      12
    );

And that’s it!

This is how I mixed and matched datasets and functions from the Earth Engine library to produce a global terrain, which may be completely useless or enormously useful. It’s up to you to decide! As for me, next time I present an indexed or a classified image, I may blend a little of that terrain to enchance the final outcome.

I have shared the entire code, which is directly accessible from your Earth Engine Code Editor, provided your Google account has been approved for Earth Engine access. Feel free to modify it or further develop it, provided you follow the requirements of its associated licence CC BY-NC-SA 4.0.

Kindest regards from Crete, Greece

Spiros