Rather than venting about a project I've recently decoupled myself from, I'll try to do something constructive... I also posted this to the gispython mailing list, but I've had to figure it out a couple times, so I'll put it here for the record. Given an N * N numpy array, and a bounding box, it's actually fairly simple to make a georeferenced tiff:
from osgeo import gdal, gdal_array
from osgeo.gdalconst import GDT_Float64
xsize, ysize = 10, 10
a = numpy.random.random((xsize, ysize)).astype(numpy.float64)
xmin, xmax = -121., -119.
ymin, ymax = 41., 43.
driver = gdal.GetDriverByName('GTiff')
# bad: out = driver.Create('a.tiff', a.shape, a.shape, 1, GDT_Float64)
# the args to Create are 'name', xsize, ysize. and .shape is rows, which is y.
driver.Create('a.tiff', a.shape, a.shape, 1, GDT_Float64)
, (xmax - xmin)/a.shape
, (ymax - ymin)/a.shape])
where the SetGeoTransform bit is (I believe) the same stuff you'd stick in a world file.
and plottable in pylab:
tif = gdal.Open('a.tiff')
a = tif.ReadAsArray()
Anyone checked out matplotlib's basemap recently? I just wish they didn't rely on geos < 3...