Using GDAL and/or NASA ASP to create a colorize hillshade

There are many good posts out there to help create colorshaded reliefs. Here I just want a quick method to create a "data" ready colorize hillshade. This means not a figure but a geospatial ready colorshade (e.g. GeoTiff) to use in a GIS or WMS server. And I need the capability to process very large GeoTiffs ( > 4GB).

Get the Tools

I recommend GDAL or NASA Ames Stereo Pipeline (ASP). Now ASP v2.6.0 actually contains the GDAL commands for this tutorial or you can use their own versions.

For GDAL binaries I recommend the Anaconda Python Environment (for all OSs). I am still digging out of Python 2.7 into Python 3.x but GDAL should work in either. First install the Anaconda environment and then run "conda install gdal" to get the gdal binaries. https://store.continuum.io/cshop/anaconda/ .

• For a more plain GDAL environment Windows binaries see: https://trac.osgeo.org/osgeo4w/
• For Mac: http://www.kyngchaos.com/software:frameworks
• For Linux: It's not too hard to download and build yourself.

If you just want ASP there are binaries for Mac and Linux: https://ti.arc.nasa.gov/tech/asr/intelligent-robotics/ngt/stereo/

Get some Data

Raster digital elevation models (DEMs / DTMs) are getting pretty easy to find. And when needed you can even create your own in ASP for some instruments. For a simple example, lets grab a small 8MB HiRISE 32bit DTM from UofA: http://www.uahirise.org/dtm/dtm.php?ID=PSP_001612_1780

wget http://www.uahirise.org/PDS/DTM/PSP/ORB_001400_001499/PSP_001414_1780_PSP_001612_1780/DTEEC_001414_1780_001612_1780_U01.IMG

Let's check out the image metadata

ast{91}> gdalinfo -mm DTEEC_001414_1780_001612_1780_U01.IMG
Driver: PDS/NASA Planetary Data System
Files: DTEEC_001414_1780_001612_1780_U01.IMG
Size is 1279, 1694
Coordinate System is:
PROJCS["EQUIRECTANGULAR MARS",
    GEOGCS["GCS_MARS",
        DATUM["D_MARS",
            SPHEROID["MARS_localRadius",3396190,0]],
        PRIMEM["Reference_Meridian",0],
        UNIT["degree",0.0174532925199433]],
    PROJECTION["Equirectangular"],
    PARAMETER["latitude_of_origin",0],
    PARAMETER["central_meridian",354.5],
    PARAMETER["standard_parallel_1",0],
    PARAMETER["false_easting",0],
    PARAMETER["false_northing",0]]
Origin = (-401.706464281687772,-120698.120618077693507)
Pixel Size = (1.011855073757400,-1.011855073757400)
Metadata:
  BANDWIDTH=
  CENTER_FILTER_WAVELENGTH=
  DATA_SET_ID="MRO-M-HIRISE-5-DTM-V1.0"
  FILTER_NAME=
  INSTRUMENT_ID="HIRISE"
...
  TARGET_NAME="MARS"
Corner Coordinates:
Upper Left  (    -401.706, -120698.121) (  5d30'24.40"W,  2d 2'10.50"S)
Lower Left  (    -401.706, -122412.203) (  5d30'24.40"W,  2d 3'54.60"S)
Upper Right (     892.456, -120698.121) (  5d29' 5.80"W,  2d 2'10.50"S)
Lower Right (     892.456, -122412.203) (  5d29' 5.80"W,  2d 3'54.60"S)
Center      (     245.375, -121555.162) (  5d29'45.10"W,  2d 3' 2.55"S)
Band 1 Block=1279x1 Type=Float32, ColorInterp=Undefined
    Computed Min/Max=-1452.895,-1373.310
  NoData Value=-3.40282265508890445e+38

Time to process using GDAL colorshade

1.) Hillshade

gdaldem hillshade DTEEC_001414_1780_001612_1780_U01.IMG Victoria_shade.tif -z 2

--note -z is the exaggeration amount

2.) Color ramp -- here we are using  a inverted diverging Blue to Brown from Color Brewer:
http://colorbrewer2.org/#type=diverging&scheme=BrBG&n=6
this look-up uses percent (as calculated min / max) and then Red Green Blue (8bit values)
create a new file and copy and paste this into it - filename: BlueBrown.lut

nv   0 0 0       //NoData
0%   1   102  94 //brown
20%  90  180 172 //
40%  199 234 229 //light brwon
60%  246 232 195 //light blue
80%  216 179 101 //
100% 140  81  10 //blue

here we are using  an inverted Yellow to Green to Blue from Color Brewer:
http://colorbrewer2.org/#type=sequential&scheme=YlGnBu&n=6
create a new file and copy and paste this into it - filename: YlGnBlue.lut

nv   0 0 0       //NoData
0%   37  52 148  //dark blue
20%  44  127 184 //med blue
40%  65  182 196 //light blue
60%  127 205 187 //light blue/green
80%  199 233 180 //light yellow/green
100% 255 255 204 //light yellow

Below you can also use discrete values in the lookup - not just percentages

nv 0 0 0               //NoData
-1452.8945  37  52 148 //dark blue
-1436.9777  44 127 184 //med blue
-1421.0608  65 182 196 //light blue
-1405.1440 127 205 187 //light blue/green
-1389.2272 199 233 180 //light yellow/green
-1373.3104 255 255 204 //light yellow

--note nv is for the NoData value - here set to black. Use 255 255 255 for white.

3.) Apply color ramp to DEM - this is a temporary file which can be deleted.

gdaldem color-relief DTEEC_001414_1780_001612_1780_U01.IMG BlueBrown.lut Victoria_ClrTemp.tif

4.1) merge the colorized file and hillshade (ASP routine)

hsv_merge Victoria_ClrTemp.tif Victoria_shade.tif -o Victoria_colorshade.tif

4.2) merge the colorized file and hillshade (GDAL Python routine in Anaconda - slower)

wget https://svn.osgeo.org/gdal/trunk/gdal/swig/python/samples/hsv_merge.py
python hsv_merge.py Victoria_ClrTemp.tif Victoria_shade.tif Victoria_colorshade.tif

OR process using single ASP routine - color blend might be a little darker

colormap -s Victoria_shade.tif --legend -o Victoria_colorshade.tif --colormap-style BlueBrown.lut DTEEC_001414_1780_001612_1780_U01.IMG

--note there is a legend ramp but you will need to add min and max text values in a photo editor app. See the min and max values from the "gdalinfo" run above.

Done.

This is not the best DEM for these colors but feel free to check out color brewer to make your own!