relieff.py

#!/usr/bin/env python

import math
import sys

import gdal
import gdalconst
import numpy
from PIL import Image

CONTRAST = 10

def main(demfile, imagefile):
    # http://geoinformaticstutorial.blogspot.no/2012/09/reading-raster-data-with-python-and-gdal.html
    dataset = gdal.Open(demfile, gdalconst.GA_ReadOnly)
    width = dataset.RasterXSize
    height = dataset.RasterYSize
    band = dataset.GetRasterBand(1)
    data = band.ReadAsArray(0, 0, width, height).astype(numpy.float)

    im = Image.new('RGB', (width-1, height-1))
    pix = im.load()

    for x in range(1, width):
        if x % 50 == 0:
            sys.stdout.write('\r[%s]' % ('-' * (x*60/width)).ljust(60))
        for y in range(1, height):
            d = data[y, x]
            h = int((d - data[y-1, x-1]) * CONTRAST + 128) if d != 0 else 0
            pix[x-1, y-1] = (h, h, h)
    
    im.save(imagefile)
    sys.stdout.write('\r[%s]\n' % ('-' * 60))

if __name__ == '__main__':
    if len(sys.argv) == 3:
        main(sys.argv[1], sys.argv[2])
    else:
        print 'usage: relieff.py input output'