GDAL geographic coordinate to pixel coordinate

Question

I've got a function that converts x,y of raster to geographic coordinates in given transformation (where gt is GetGeoTransform() from GDAL):

def p(col, row, gt): # p: pixel coords to map coords
    c, a, b, f, d, e = gt
    x_geo = a * col + b * row + a * 0.5 + b * 0.5 + c
    y_geo = d * col + e * row + d * 0.5 + e * 0.5 + f
    return x_geo, y_geo  # map coordinates

But how do I reverse it, so it would be geographic coordinates to pixel coordinates? Becasue mathematically it doesn't seem to work correctly:

def p_reverse(col, row, gt):
    c, a, b, f, d, e = gt
    x_geo = a / col - b / row - a / 0.5 - b / 0.5 - c
    y_geo = d / col - e / row - d / 0.5 - e / 0.5 - f
    return x_raster, y_raster

What am I doing wrong?

Answer 1

from osgeo import gdal
layer = gdal.Open('a.tif')
gt =layer.GetGeoTransform()
# from your function
col = 100
row = 75
x_geo = a * col + b * row + a * 0.5 + b * 0.5 + c
y_geo = d * col + e * row + d * 0.5 + e * 0.5 + f
print(x_geo,y_geo)
159531.45050461384 77980.4986499551
# reverse
col = int((x_geo - c) / a)
row = int((y_geo - f) / e)
print(col, row)
100 75

Answer 2

Using numpy, your "forward" transform could be written like this:

import numpy as np
P = np.array([col + .5, row + .5, 1])
M = np.array([[a, d, 0],
              [b, e, 0],
              [c, f, 1]])
G = P @ M # this will be [x_geo, y_geo, 1]

To go in the opposite direction, you need to multiply a geographic coordinate by the inverse of M, and then subtract the half-pixel offset:

P1 = np.array([x_g1, y_g1, 1]) @ np.linalg.inv(M) - 0.5

where x_g1 and y_g1 are your geographic coordinates.

(Not tested, so I have probably got the matrix the wrong way around. Anyway, the key concepts are "matrix inverse" and "homogenous coordinates". Time to get out your old linear algebra textbook — or get a new one!)