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).
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/ .
Let's check out the image metadata
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
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
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
Done.
This is not the best DEM for these colors but feel free to check out color brewer to make your own!