# GDAL geographic coordinate to pixel coordinate

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?

• – gene
Feb 20, 2022 at 14:45
• I'm sorry but that is not clear enough for me, I still don't know how to correct my function. I really wouldn't like to use Rasterio. Feb 20, 2022 at 14:55

``````from osgeo import gdal
layer = gdal.Open('a.tif')
gt =layer.GetGeoTransform()
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
``````
• What if b or d are non-zero? Feb 20, 2022 at 15:38
• Use the solution of perrygeo
– gene
Feb 20, 2022 at 15:41
• @TurePålsson good question... nothing about it on perrygeo Feb 24, 2022 at 15:56
• Very helpful, just reverse without int() gives 0.5 pixel more. I ended up doing like this `col = (x_geo - upper_left_x - x_size * 0.5 )/ x_size ` `row = (y_geo - upper_left_y - y_size * 0.5 )/ y_size ` Nov 11, 2023 at 1:57
• just to clear things up, adding 0.5 is not always good, that depends on which quarter of globe you are, sometimes you need to put "-" on x or y or both. Nov 30, 2023 at 15:45

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!)

• But I don't have a Numpy array with coordinates, but points from SHP ... I need to have a solution per point ... got lat and lng, not some Numpy array. I am really sorry, but I don't consider this as an answer to my my question. From lng and lat I need an x,y when giving transformation array from GDAL. Feb 20, 2022 at 15:23
• In my example, `x_g1` and `y_g1` would be your longitude and latitude. Feb 20, 2022 at 15:27