best way(s) to fit a polynomial to points?

[Kyler Laird]

I'm trying to do some georeferencing - using points of known
location (ground control points, GCPs) on an image to develop
a polynomial that can be used to approximate the locations of
other points. 

I usually use the Python bindings to GDAL
	http://www.remotesensing.org/gdal/
to manipulate geospatial data but there are no Python bindings
for GDALCreateGCPTransformer().
	http://remotesensing.org/cgi-bin/cvsweb.cgi/~checkout~/osrs/gdal/alg/gdal_crs.c?rev=1.4.2.1&content-type=text/plain
I tried using it from C (both through GDAL and directly to
the GRASS code) and failed.

I'd rather have a Python-based solution anyway, so I'm
searching for appropriate tools.  The ones I'm finding (like
matplotlib.mlab.polyfit() and PyBLD's polfit()) only deal with
single variable equations (f(x) = X).  I need something that
operates on two variables (f(x, y) = X).  (I'd have two
polynomials; one would calculate the northing and the other
would do the easting).  I need to be able to fit at least
third degree polynomials.

The solution does not need to be especially fast.  A pure
Python solution is preferred but I'll take about anything
that works right now.  (I'm trying to show someone how to do
this.  I've only used proprietary software to do it but I've
been meaning to write my own.)

Any suggestions?

Thank you.

--kyler