Fast Reversible Encoding/Decoding of DEM (16bit) into RGB

Copyright (C) 2011, Klokan Technologies GmbH. Author: Petr Pridal, Vaclav Klusak, Petr Sloup

The problem we are facing in WebGL Earth implementation of the terrain:

WebGL seems to not be able to read the 16-bit textures, which are normally used in OpenGL ES 2.0 for rendering of the terrain with elevation (in the vertex shader). To go over this limit we need to encode standard elevation model of Earth (-10,701m to 8,248m) into RGB space and prepare tile pyramid from the available data (such as SRTM3 Plus). The encoding algorithm should be:

  1. very fast (only bit shifting and logical operations or basic algebra if possible, simpler means better)
  2. symmetric (used colors must be decodable 1:1 to the elevation model and vice versa - without any loss)
  3. optional: friendly to lossy compression (if the DEM data are saved in an image format which is using lossless compression then the differences between input elevation and decoded elevation should be minimal) - if this is not possible then lossless compression must be used (such as PNG)
  4. optional: friendly to humans (people naturally expect sea with blue tones, land with green/yellow/brown, bright to dark for actual elevation values) - this can be practical for visual validation of correctness of the dataset too
Here we present possible solutions. Please contact us if you know about a better one.

Alternative A: Simple bit encoding of a float into RGB 8-bit

A naive encoding of a float with elevation with known range into R and G channels of the output.
Decoding algorithm:

r =[0] / 256.0
g =[1] / 256.0
elev = Math.ceil((TERRAIN_MAX-TERRAIN_MIN)*(r + g/256.0) - TERRAIN_MAX)

Encoding is possible with the utility gdaldem_web:

The alternative A is what we are now using in WebGL Earth.

Alternative B: Z-order curve

Terrestrial Canvas
Bathymetry Canvas
The elevation:

The suggested algorithm: The Z-order curve (known as Morton order) which goes trough the RGB cube in 2 dimensions - with the blue color set to minimal value for land and maximal for water. We have available value range (65536) which equals to 25 cm steps for elevation model. The tiles are rendered in popular mercator system, same as the map data such as OpenStreetMap.

Encoder/Decoder Source code:

function part1by1(n) {
    n &= 0x0000ffff
    n = (n | (n << 8)) & 0x00FF00FF
    n = (n | (n << 4)) & 0x0F0F0F0F
    n = (n | (n << 2)) & 0x33333333
    n = (n | (n << 1)) & 0x55555555
    return n

function unpart1by1(n) {
    n &= 0x55555555
    n = (n ^ (n >> 1)) & 0x33333333
    n = (n ^ (n >> 2)) & 0x0f0f0f0f
    n = (n ^ (n >> 4)) & 0x00ff00ff
    n = (n ^ (n >> 8)) & 0x0000ffff
    return n

function interleave2(x, y) {
    return part1by1(x) | (part1by1(y) << 1)

function deinterleave2(n) {
    return [ unpart1by1(n), unpart1by1(n >> 1) ]

Note: this code is port of an adaptation of the recipes from Sean Eron Anderson and Fabian "ryg" Giesen prepared by Romain Dartigues; all credits goes to the respective authors.