# Calculating inverse polynomial transforms for pixel sampling when map georeferencing

In map georeferencing, one typically chooses a set of ground control points that represent pixel locations in the map image along with their real world geographic coordinates, which are then used to estimate a polynomial forward transformation of 1st/affine, 2nd, or 3rd order, for instance using least squares. The forward transform is then used to calculate the bounds of the geographic coordinate system.

Next, one needs the equivalent backwards transformation for determining which pixel coordinate in the original map image to "sample" for each of the pixels in the output warped/resampled map image.

My question is what is the correct or common approach used to calculate the backwards/inverse of polynomial transforms when map georeferencing? How is it done in common GIS software/libraries such as gdalwarp? The forward transform is simple enough using a least squares approach, but I can't figure out how to get the inverse transformation, especially for 2nd and 3d order polynomials.

My solution so far has been to simply estimate a separate backwards transformation using least squares to predict the image coordinates based on the geographic coordinates of the GCPs (instead of the other way around). The problem with this however, is that reestimating a new transformation does not create an exact inverse of the original forward transform, i.e. forward and then backward transforming an image pixel coordinate does not return the same image pixel that I started with, it ends up slightly offset. For the image transformation to be accurate the backwards sampling needs to be the exact inverse of the forward transform, so you can go back and forth without losing information.

Is there another mathematical or programmatic approach to get the exact inverse of such polynomial map transforms, involving for instance matrix algebra or numpy? I'm not looking to outsource the problem to external libraries, I want to implement the nuts and bolts, preferably in numpy.

I've dived under the hood of GDAL and gdalwarp but so far I've been unable to find the source code for calculating the inverse transform, so any pointers here would be useful as well.

No need for estimations, we can invert the matrix exactly using GDAL itself:

``````def get_index_for_GPS(lat, lon):
transform = ds.GetGeoTransform()
reverse = gdal.InvGeoTransform(transform)
x, y = [int(v) for v in gdal.ApplyGeoTransform(reverse, lon, lat)]
return x, y
``````

However, not all ports of GDAL expose this function, so we can also exactly compute the reverse operation by applying a common trick. We first observe that the forward transform is a list of six numbers, defining the following operation:

``````xgeo = a + b * xpix + c * yline
ygeo = d + e * xpix + f * yline
``````

That is:

``````                     ⎛   1   ⎞
⎛ xgeo ⎞ = ⎛ a  b  c ⎞ ⎜ xpix  ⎟
⎝ ygeo ⎠   ⎝ d  e  f ⎠ ⎝ yline ⎠
``````

So, if we want the reverse operation, we just invert the forward transform matrix. Except we can't, because only square matrices can be inverted. Thankfully, we can make our transform square:

``````⎛  1   ⎞   ⎛ 1  0  0 ⎞ ⎛   1   ⎞
⎜ xgeo ⎟ = ⎜ a  b  c ⎟ ⎜ xpix  ⎟
⎝ ygeo ⎠   ⎝ d  e  f ⎠ ⎝ yline ⎠
``````

The top row doesn't actually "do" anything here other than pad the matrix to square form, which is perfect: we now have a matrix that we can feed into any matrix inversion utility to get the inverse. For example, if we ask WolframAlpha, we get the following inverse:

But depending on whether there are any zeroes in our forward transform - and there pretty much always are - we can get a much more easily computed inverse. For instance, the GDAL forward transform almost always has zeroes for c and e:

``````⎛  1   ⎞   ⎛ 1  0  0 ⎞ ⎛   1   ⎞
⎜ xgeo ⎟ = ⎜ a  b  0 ⎟ ⎜ xpix  ⎟
⎝ ygeo ⎠   ⎝ d  0  f ⎠ ⎝ yline ⎠
``````

And that, we can trivially invert. By hand even: first we left-pad the matrix with a 3x3 identity matrix:

``````⎛ 1  0  0  ⎜  1  0  0 ⎞
⎜ 0  1  0  ⎜  a  b  0 ⎟
⎝ 0  0  1  ⎜  d  0  f ⎠
``````

and then we row-reduce. Subtract (a * first row) from the second row, and (d * first row) from the third row:

``````⎛  1  0  0  ⎜  1  0  0 ⎞
⎜ -a  1  0  ⎜  0  b  0 ⎟
⎝ -d  0  1  ⎜  0  0  f ⎠
``````

then divide the second row by b, and divide the third row by f:

``````⎛    1    0    0  ⎜  1  0  0 ⎞
⎜ -a/b  1/b    0  ⎜  0  1  0 ⎟
⎝ -d/f    0  1/f  ⎜  0  0  1 ⎠
``````

And done: we've now "moved" the identity matrix from the left to the right, so the values on the left "reverse" the original matrix we had on the right. Let's remove that 3x3 identity on the right, and remove the top row, and presto: our inverse operation is:

``````                  ⎛  1   ⎞
⎛ -a/b  1/b   0  ⎞ ⎜ xgeo ⎟ = ⎛ xpix  ⎞
⎝ -d/f   0   1/f ⎠ ⎝ ygeo ⎠   ⎝ yline ⎠
``````

and we can feed GDAL (etc) that reversed matrix directly:

``````def get_index_for_GPS(lat, lon):
forward = ds.GetGeoTransform()

# find the reverse transform for this forward transform...
[a,b,c,d,e,f,] = forward
if c == 0 and e == 0:
# Use the nice and easy reverse transform.
reverse = [-a/b, 1/b, 0, -d/f, 0, 1/f]
else:
# Use the big guns to compute it, instead.
m = np.array([[1,0,0],[a,b,c],[d,e,f]])
inverted = np.linalg.inv(m)
reverse = list(inverted[1]) + list(inverted[2])

# And now we can get our pixel coordinate given our geo coordinate.
x, y = [int(v) for v in gdal.ApplyGeoTransform(reverse, lon, lat)]
return x, y
``````