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

demgdalgeotiff-tiffglobal-mappertiles

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.

Best Answer

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.

Related Question