Fusión de datos hillshade / DEM en relieve de color geotiff único con QGIS y GDAL

Question

Tengo un DEM geotiff de una sola banda y he conseguido crear una versión en color del mismo (con pesudocolor de una sola banda) y un hillshade (con Terrain Analysis --> Hillshade). Según la lección en línea ( https://docs.qgis.org/2.2/en/docs/training_manual/rasters/terrain_analysis.html ), puedo superponer la sombra de la colina sobre el mapa de colores y queda genial. Mi único problema es que quiero exportarlo fuera de QGIS como un único GeoTIFF unificado para su uso en otros softwares.

Este ( http://dirkraffel.com/2011/07/05/best-way-to-merge-color-relief-with-shaded-relief-map/ ) sugiere usar el conversor de ColorMagick, pero esto resulta en archivos muy grandes y a menudo se bloquea, por no mencionar que los datos de coordenadas se pierden.

Este ( https://anitagraser.com/2012/01/19/a-guide-to-beautiful-reliefs-in-qgis/ ) es otra opción, pero me resulta muy difícil poner escalas de color adecuadas y no se puede controlar el sombreado.

Este ( http://blog.mastermaps.com/2012/07/terrain-mapping-with-mapnik.html ) depende de Mapnik, que actualmente no es una opción para mí.

¿Puede alguien sugerir una solución eficaz que permita mezclar dos geotiffs (grandes) superpuestos con el superior (hillshade) ligeramente transparente?

Para quien esté interesado, he improvisado este script GDAL a partir de la primera fuente, que hace el trabajo a medias:

   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 es una buena herramienta para realizar operaciones y funciones matemáticas en imágenes rasterizadas. Las funciones de gamma y superposición son bastante sencillas de escribir y puedes encadenar algunos archivos temporales para segmentar la complejidad. Lo mejor de todo es que puedes conservar la geolocalización del archivo.

la función gamma parece :

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

el modo de fusión de la superposición parece :

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

con a y b las imágenes inferior y superior (cada una normalizada entre 0 y 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

Puede utilizar --outfile=colored_hillshade.png si desea que el archivo sea más transportable.