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

[Cameron Laird]

In article <iuu4g1-kl1.ln1 at jowls.lairds.org>,
Kyler Laird  <Kyler at news.Lairds.org> wrote:
>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.)
			.
			.
			.
"I tried using it ... and failed."  I naturally wonder *how*
you (it, more likely) failed--the specific details of what
you actually observed.  I can imagine, though, that if you 
had the time to understand and describe what happened with
precision, it'd be easier fixing it than talking about it.
We can live this part aside for now.

Curve-fitting is a domain where a lot of subtlety accompanies
the transition from one to two dimensions.  Let me repeat in
my own words some of what you've said, to be certain I under-
stand:  you have a table of observations, in columns

  x      y     X

with a large number of rows.  You have reason to believe that
the table can be adequately modeled by a polynomial of low 
degree in x and y, jointly, plus small non-systematic errors.
The use of this polynomial will be to estimate new X-s from
new values of x and y.

Will the new (x,y)-s be interior to the convex hull formed by
the tabulated ones, or might they be outside?

I don't know PyBLD and matplotlib.  It's *very* likely, though,
that there's a simple way to restate the problem so that the
least-squares' solvers they must embed can be applied to 
simply transformed data to yield polynomial coefficients.  How
soon does this thing need to go live?
-- 

Cameron Laird <claird at phaseit.net>
Business:  http://www.Phaseit.net