I recently had to convert some geographical data given as points on a map projection (x,y) to geographical coordinates (lon, lat). As the map projection was unknown, I resorted to fitting a polynomial to some selected points, and transforming the remainder of the data set using the fit. This can be achieved easily using linear least squares, but writing out terms in the bivariate expansion can quickly get tedious.
The code below implements the LLS algorithm for bivariate polynomials, along with a helper function to calculate the results of the fit.
