2

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?

2

2 Answers 2

3
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
5
  • 1
    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
    – nadya
    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
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!)

2
  • 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
  • 1
    In my example, x_g1 and y_g1 would be your longitude and latitude. Feb 20, 2022 at 15:27

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.