Starting on gdal for SRTM data

Awhile back I made a slippy map of parts of the US over 8000′. I did it using SRTM elevation data turned into images with mapserver, a dynamic program. Easy to get going but awfully slow to serve. I’d like to make the same map but statically. Inspired by Matt’s success with gdal2tiles.py I thought I’d try again.

My source data is the .hgt files for North America. They’re bundled together into a single index.shp file using gdaltindex, as described in the mapserver docs. A quick and dirty render is accomplished with

gdal2tiles.py N37W119.hgt ~/public_html/tmp/gdal2tiles

That does all the expected tiling and rendering at multiple scales, but doesn’t have a meaningful colourmap. But in my work with Mapserver I created a meaningful colormap:

  LAYER
    NAME         srtm
    TILEINDEX    "index.shp"
    TYPE         RASTER
    STATUS       DEFAULT
    PROCESSING   "SCALE=2438,4572"    # 914m = 3000', 2438m = 8000', 4572m = 15000'
    PROCESSING   "BANDS=1,1,1,1"
    PROCESSING   "LUT_1=0:0,1:0,73:80, 74:130, 255:180"
    PROCESSING   "LUT_2=0:0,1:0,73:80, 74:0, 255:0"
    PROCESSING   "LUT_3=0:0,1:0,73:80, 74:0, 255:0"
    PROCESSING   "LUT_4=0:0,1:255,255:255"

What that says, in mapserver-eze, is:

  1. For each input point from SRTM, map the range 2438:4572 to 0:255.
    This maps everything < 8000′ to 0, everything over 15000′ to 255, and a linear range between.
  2. Map the input band “1” (the altitude) to the output channels RGBA. Ie, a full 32 bit image with 8 bits of alpha.
  3. LUT_4 maps everything that is 0 (< 8000′) to transparent, everything > 0 to fully visible (alpha = 255)
  4. LUT_1-3 apply a colourmap. 1 (8000′) is mapped to 0,0,0 (black), 73 (10000′?) maps to 80,80,80 (grey), and 255 (15000′) maps to red. This gives a color ramp: black to grey to red.

Now I need to figure out how to generate the same colourmap using solely GDAL. And then convince gdal2tiles to do that colourmapping, probably by first turning my HGT files into GeoTIFF and then tiling the GeoTIFF.