gistfile1.txt

#!/usr/bin/env python

import gdal
import gdalnumeric
import sys
import numpy

# =============================================================================
def Usage():
    print('Usage: elev2rgb.py elevation_file blue_file output_file')
    print('')
    sys.exit(1)

# =============================================================================

argv = gdal.GeneralCmdLineProcessor( sys.argv )

elev_filename = None
blue_filename = None
out_filename = None

# Parse command line arguments.
i = 1
while i < len(argv):
    arg = argv[i]

    if elev_filename is None:
      elev_filename = arg

    elif blue_filename is None:
      blue_filename = arg

    elif out_filename is None:
      out_filename = arg

    else:
        Usage()

    i = i + 1

if out_filename is None:
    Usage()

elev_db = gdal.Open( elev_filename )
blue_db = gdal.Open( blue_filename )

drv = gdal.GetDriverByName('GTiff')
outdb = drv.Create( out_filename, elev_db.RasterXSize, elev_db.RasterYSize, 3,
                    gdal.GDT_Byte )
outdb.SetProjection( elev_db.GetProjectionRef() )
outdb.SetGeoTransform( elev_db.GetGeoTransform() )

outdb.GetRasterBand(3).WriteArray(blue_db.GetRasterBand(1).ReadAsArray())

elev_data = elev_db.ReadAsArray()
elev_data = numpy.multiply(elev_data, 3)
#elev_data = elev_data + 10000

#outdb.GetRasterBand(1).WriteArray(elev_data / 256)
outdb.GetRasterBand(1).WriteArray(numpy.floor(elev_data / 256))
outdb.GetRasterBand(2).WriteArray(elev_data % 256)