create a color palette for the gdaldem color-relief command

gistfile1.txt

import math
import colorsys

MIN = -500 # m
MAX = 9000 # m
CYCLE = 500 # m

"""
make an elevation color palette suitable for global use, while still preserving
small variations.

the eye is most sensitive to luminosity, so the main component is zig-zagging hsv "value"
between max and min once every 'CYCLE' meters

then we vary saturation over the same cycle to distinguish the ascending and descending
phases of luminosity.

lastly, hue creeps forward to distinguish different cycles from each other. initial hue
was changed from red to orange because the initial result was a bit too... bloody

also add a phase discontinuity at sea-level to easily distinguish landforms
"""

def mix(k, min, max):
    k = .5 * k + .5
    return (1-k)*min + k*max

def hsv(elev):
    k = float(elev) / (MAX - MIN)
    theta = (float(elev) / CYCLE + (.2 if elev <= 0 else 0)) * 2*math.pi

    kval = 1 - 2*abs(1 - (theta / math.pi) % 2.)
    ksat = math.sin(theta)

    hue = (k + .08) % 1.
    val = mix(kval, .02, 1.)
    sat = mix(ksat, .2, 1.)

    return (hue, sat, val)
    
def print_step(elev):
    color = colorsys.hsv_to_rgb(*hsv(elev))
    print elev, ' '.join(str(int(256*c)) for c in color)

STEP = float(CYCLE) / 200
for i in range(int(math.floor(MIN / STEP)), 1):
    print_step(i * STEP)
print_step(.001)
for i in range(1, int(math.ceil(MAX / STEP)) + 1):
    print_step(i * STEP)