plot_geotiff

#!/usr/bin/env python

from osgeo import gdal
import matplotlib.pyplot as mpl
import sys
import numpy.ma as ma

def plotit (filename):
    ds = gdal.Open(filename, gdal.GA_ReadOnly)
    (X, deltaX, rotation, Y, rotation, deltaY) = ds.GetGeoTransform()
    srs_wkt = ds.GetProjection()
    Nx = ds.RasterXSize
    Ny = ds.RasterYSize
    arys=[]
    for i in xrange(1, ds.RasterCount+1):
        data = ds.GetRasterBand(i).ReadAsArray()
        ndv = ds.GetRasterBand(i).GetNoDataValue()
        if  ndv:
            data = ma.masked_values(data, ndv)
        else:
            ndv=-9e9
        data=ma.masked_invalid(data)
        print data
        arys.append(data)
    cmap = mpl.cm.afmhot
    cmap.set_bad('w',1.)
    print arys[0][2,2]
    for band  in  arys:
        mpl.figure()
        mpl.imshow(band,interpolation="nearest", cmap=cmap)
        mpl.colorbar()
        mpl.show()

for i in sys.argv[1:]:
    plotit(i)