numpy to tiff via gdal

EDIT: 7 months later I came back to this and found an error. update in the code below with old line commented out.

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
import numpy
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[0], a.shape[1], 1, GDT_Float64)
# the args to Create are 'name', xsize, ysize. and .shape[0] is rows, which is y.
driver.Create('a.tiff', a.shape[1], a.shape[0], 1, GDT_Float64)

out.SetGeoTransform([xmin
                  , (xmax - xmin)/a.shape[0]
                  , 0
                  , ymin
                  , 0
                  , (ymax - ymin)/a.shape[1]])

gdal_array.BandWriteArray(out.GetRasterBand(1), a)

where the SetGeoTransform bit is (I believe) the same stuff you'd stick in a world file.

and plottable in pylab:

import pylab
tif = gdal.Open('a.tiff')
a = tif.ReadAsArray()
pylab.imshow(a)
pylab.show()

Anyone checked out matplotlib's basemap recently? I just wish they didn't rely on geos < 3...