How to transform a GeoTiff DEM to png tiles encoding elevation as RGB colors

Question

I want to convert a GeoTiff DEM file, to PNG tiles where the elevation is encoded as RGB colors with a formula like that :

-10000 + ((color.r * 256 * 256 + color.g * 256 + color.b) * 0.1)

I know how to do the raster Geotiff to PNG tiles with Global Mapper or GDAL, but I don't know how to make to colourize the elevation with tri-linear interpolation like the formula.

I tried with Global Mapper and a grayscale shader between min an max elevation, but the elevation is as a consequence only encoded in 8bit, and it is not precise enough.

The same with gdaldem color-relief, the color configuration file only do linear interpolation between 2 elevations and colors so it is the same as a grayscale and result in 8bit precision. Or i would have to give a 255*255 long color configuration file to express the tri-linear interpolation as a linear interpolation.

I am not sure if gdal-translate can expand the elevation band as 3 band in RGB.

I am thinking of doing a python program to do it now, but before losing time I want to know if there is other options.

Answer 1

To convert a GeoTIFF into png, I use a 2-step process. In the first step, I use a small Python script (depends on rasterio) to translate the elevation band of the original GeoTIFF into r, g, and b bands. The result is another GeoTIFF, which I use as input for gdal2tiles.py to create a PNG tile set.

The Python code to encode the elevation into r, g, and b looks like this:

with rasterio.open('infile.tif') as src:
    dem = src.read(1)

r = np.zeros(dem.shape)
g = np.zeros(dem.shape)
b = np.zeros(dem.shape)

r += np.floor_divide((100000 + dem * 10), 65536)
g += np.floor_divide((100000 + dem * 10), 256) - r * 256
b += np.floor(100000 + dem * 10) - r * 65536 - g * 256

meta = src.meta
meta(
    dtype=rasterio.uint8,
    nodata=0,
    count=3)

with rasterio.open('outfile.tif', 'w', **meta) as dst:
    dst.write_band(1, r.astype(rasterio.uint8))
    dst.write_band(2, g.astype(rasterio.uint8))
    dst.write_band(3, b.astype(rasterio.uint8))

The png tiles generated from outfile.tif will have the elevation encoded. To decode, apply the formula

-10000 + (pixel.r * 256 * 256 + pixel.g * 256 + pixel.b) * 0.1

to the pixel you want to get the elevation for.

Answer 2

You can accomplish that in Global Mapper creating a custom shader. Go to the "Shader Options" tab on the configuration windows and in the section "Custom Shaders" select "New". There you can create a multi-linear sahder selecting different colors for different elevations. In the example below I did create one going red-green between 0 and 3000 m and green-blue between 3000 and 6000 m Then, after selecting this shader you just have to export your DEM as raster and specify the tiling on the "Tiling" tab of the PNG export options.

Answer 3

You can do this in GRASS GIS. If you are not familiar with GRASS GIS, you can read more about how to get started here. The explanation below assumes you run the commands from the command line. However, you can find the same functions in the menu if you prefer to work with the GUI.

Let's assume you have a GRASS GIS database and your geotiff file is called 'dem'.

Step 1: First step is to start GRASS in the location and mapset with a projection that matches that of your layer.

Step 2: Import your geotiff in a GRASS GIS database. Run the g.region function to ensure that the working resolution and extent matches that of the imported layer.

r.in.gdal input=dem.tif output=dem
g.region raster=dem

Step 3: You use the r.mapcalc function to create a new layer according to your equation.

r.mapcalc "dem2 = -10000 + ((r#dem * 256 * 256 + g#dem * 256 + b#dem) * 0.1)"

Step 4: You can now export the newly created layer as geotif (or any other format supported by GRASS):

r.out.gdal input=dem2 output=dem2.tif format=GTiff

For a detailed explanation of the syntax used in the r.mapcalc function, see the r.mapcalc manual page