Merging hillshade / DEM data into color relief single geotiff with QGIS and GDAL

Question

I have a single-band geotiff DEM and I have managed to create a color version of it (with single-band pesudocolor) and a hillshade (with Terrain Analysis --> Hillshade). As per the online lesson (https://docs.qgis.org/2.2/en/docs/training_manual/rasters/terrain_analysis.html), I can overlap the hillshade over the color map and it looks great. My only problem is I want to export it out of QGIS as a single unified GeoTIFF for use in other softwares.

This (http://dirkraffel.com/2011/07/05/best-way-to-merge-color-relief-with-shaded-relief-map/) suggests using ColorMagick's convert, but this results in very large files and often crashes for me, not to mention the coordinate data gets lost.

This (https://anitagraser.com/2012/01/19/a-guide-to-beautiful-reliefs-in-qgis/) is another option, but I find it very hard to put fitting color scales and the hillshading cannot be controlled.

This (http://blog.mastermaps.com/2012/07/terrain-mapping-with-mapnik.html) relies on Mapnik, which is currently not an option for me.

Can anyone suggest an efficient solution which allows to blend two (large) overlapping geotiffs with the top one (hillshade) being slightly transparent?

For anyone interested, I cobbled together this GDAL script from the first source, which gets the job half-done:

   set input_dem=test-mbes2.tif
   set input_hill=test-illumination.tif

   convert -gamma .5 %input_hill% hills_gamma.tif
   convert %input_dem% hills_gamma.tif -compose Overlay -composite output.tif
   listgeo %input_dem% > meta.txt
   geotifcp -g meta.txt output.tif final.tif 

Answer 1

gdal_calc.py is a good tool to perform mathematical operations and functions on raster images. The gamma and overlay functions are fairly simple to write and you can string together some temporary files to segment the complexity. Best of all you can retain the geolocation of the file!

the gamma function looks like:

encoded = ((original / 255) ^ (1 / gamma)) * 255

the overlay blend mode looks like:

2*a*b if a<0.5 else (1-2*(1-a)*(1-b))

with a and b the lower and upper images (each normalized between 0 and 1)

# gamma hillshade
gdal_calc.py -A hillshade.tif --outfile=gamma_hillshade.tif \
  --calc="uint8(((A / 255.)**(1/0.5)) * 255)"



# overlay
 gdal_calc.py -A gamma_hillshade.tif -B color.tif --allBands=B \
--calc="uint8( ( \
                 2 * (A/255.)*(B/255.)*(A<128) + \
                 ( 1 - 2 * (1-(A/255.))*(1-(B/255.)) ) * (A>=128) \
               ) * 255 )" --outfile=colored_hillshade.tif

You may wish to use --outfile=colored_hillshade.png if you want to make the file more transportable.