Find pixel coordinates from lat/long point in .geotiff using python and gdal

Question

I have (lat,long) coordinate describing the position of a point in a .geotiff image.

I wish to find the equivalent pixel coordinates of the lat,long ones inside the image.

I succeded using gdaltransform from the command line with the following instruction :

gdaltransform -i -t_srs epsg:4326 /path/imagename.tiff
-17.4380493164062 14.6951949085676

But i would like to retrieve such type of equivalence from python code. I tried the following :

from osgeo import osr

source = osr.SpatialReference()
source.ImportFromUrl(path + TIFFFilename)

target = osr.SpatialReference()
target.ImportFromEPSG(4326)

transform = osr.CoordinateTransformation(target,source )

point_xy = np.array(transform.TransformPoint(-17.4380493164062,14.6951949085676))

But it returns this error :

NotImplementedError: Wrong number or type of arguments for overloaded function 'CoordinateTransformation_TransformPoint'.
  Possible C/C++ prototypes are:
    OSRCoordinateTransformationShadow::TransformPoint(double [3])
    OSRCoordinateTransformationShadow::TransformPoint(double [3],double,double,double)

What am i doing wrong ? I tried to work around this error but without any success. Is there an other way to do it ?

EDIT 1 :

I achieved a single transformation via gdaltransform commands in terminal :

gdaltransform -i -t_srs epsg:4326 /path/image.tiff
-17.4380493164062 14.6951949085676

As i need to retrieve the pixel in a pythonic way, i tried calling the command using subprocess like :

# TRY 1:
subprocess.run(['gdaltransform','-i',' -t_srs','epsg:4326','/pat/img.tiff\n'], stdout=subprocess.PIPE)
# TRY 2 :
cmd = '''gdaltransform -i -t_srs epsg:4326 /home/henri/Work/imdex_visio/AllInt/Dakar_X118374-118393_Y120252-120271_PHR1A_2016-03-10T11_45_39.781Z_Z18_3857.tiff
-17.4380493164062 14.6951949085676'''
subprocess.Popen(cmd,stdout=subprocess.PIPE, shell=True)

But it does not work. Maybe because of the way the command itself behaves, like not actually returning a result and ending itself, but displaying the result and staying busy.

Have your tried just parsing a third argument to TransformPoint? Like `transform.TransformPoint(-17.4380493164062, 14.6951949085676, 0)` — the_cheff, Oct 31 '19 at 06:39
Can you check GDAL version? It might be that in an older implementation of python bindings, the method has a different signature. — cosmarc, Nov 07 '19 at 08:39
@cosmarc Here is my version : ```$ ogrinfo --version GDAL 2.4.2, released 2019/06/28 ``` — drheinrich940, Nov 12 '19 at 07:40
@the_cheff I had already tried this, and just did it again to make sure, but it returns me the same error : ```--------------------------------------------------------------------------- NotImplementedError Traceback (most recent call last) in ----> 1 transform.TransformPoint(-17.4380493164062, 14.6951949085676, 0) Wrong number or type of arguments for overloaded function 'CoordinateTransformation_TransformPoint ``` — drheinrich940, Nov 12 '19 at 07:42

the_cheff · Accepted Answer · 2019-11-12T11:14:55.673

6

According to the cookbook you are flipping the use of transform and point. You should call transform on the point given the transform, not the other way around. It also seems like you are flipping source and target, but you do that two times, so it will work.

However I believe that target.ImportFromUrl(path + TIFFFilename) will not work. Instead you can extract the spatial reference from the geotiff using gdal.

Something like the following should work

from osgeo import osr, ogr, gdal

# Extract target reference from the tiff file
ds = gdal.Open(path + TIFFFilename)
target = osr.SpatialReference(wkt=ds.GetProjection())

source = osr.SpatialReference()
source.ImportFromEPSG(4326)

transform = osr.CoordinateTransformation(source, target)

point = ogr.Geometry(ogr.wkbPoint)
point.AddPoint(-17.4380493164062, 14.6951949085676)
point.Transform(transform)
print(point.GetX(), point.GetY())

This provides you with the coordinates in your geotiffs reference, however this is not pixel coordinates.

To convert the point to pixels you could use something like the following (the minus for the line might have to be flipped, based on where in the world you are)

def world_to_pixel(geo_matrix, x, y):
    """
    Uses a gdal geomatrix (gdal.GetGeoTransform()) to calculate
    the pixel location of a geospatial coordinate
    """
    ul_x= geo_matrix[0]
    ul_y = geo_matrix[3]
    x_dist = geo_matrix[1]
    y_dist = geo_matrix[5]
    pixel = int((x - ul_x) / x_dist)
    line = -int((ul_y - y) / y_dist)
    return pixel, line

So your final code would look something like

from osgeo import osr, ogr, gdal


def world_to_pixel(geo_matrix, x, y):
    """
    Uses a gdal geomatrix (gdal.GetGeoTransform()) to calculate
    the pixel location of a geospatial coordinate
    """
    ul_x= geo_matrix[0]
    ul_y = geo_matrix[3]
    x_dist = geo_matrix[1]
    y_dist = geo_matrix[5]
    pixel = int((x - ul_x) / x_dist)
    line = -int((ul_y - y) / y_dist)
    return pixel, line

# Extract target reference from the tiff file
ds = gdal.Open(path + TIFFFilename)
target = osr.SpatialReference(wkt=ds.GetProjection())

source = osr.SpatialReference()
source.ImportFromEPSG(4326)

transform = osr.CoordinateTransformation(source, target)

point = ogr.Geometry(ogr.wkbPoint)
point.AddPoint(-17.4380493164062, 14.6951949085676)
point.Transform(transform)

x, y = world_to_pixel(ds.GetGeoTransform(), point.GetX(), point.GetY())
print(x, y)

edited Nov 12 '19 at 11:14

answered Nov 12 '19 at 08:17

the_cheff

4,690
3
15
23

Thanks a lot for your reply and explanation, i tried to run your code and everything was smooth until the ```point.Transform(transform)``` call. The error is the following : ```~/anaconda3/lib/python3.7/site-packages/osgeo/ogr.py in Transform(self, *args) 7369 OGRERR_NONE on success or an error code. 7370 """ -> 7371 return _ogr.Geometry_Transform(self, *args) 7372 7373 TypeError: in method 'Geometry_Transform', argument 2 of type 'OSRCoordinateTransformationShadow *' ``` ....I'm currently trying an alternative with the gdal.Transformer – drheinrich940 Nov 12 '19 at 10:04
1

It is probably due to the fact that your `target` is not correctly imported. For instance try to replace `target.ImportFromUrl(path + TIFFFilename)` with `target.ImportFromEPSG(4088)` and verify that it works (with wrong transformation though). If you know the target EPSG you can simply insert it. I do not think that importfromurl is doing what you think it does ;-) – the_cheff Nov 12 '19 at 10:22
I have updated the code in my answer to extract the Spatial Reference from your geotiff, which is what I think you are trying to do – the_cheff Nov 12 '19 at 10:31
I just achieved the same thing after doing some researches ! Thanks a lot for your efforts in helping me. Now i am quite confused because what i am trying to do is find the pixel coordinates of a given position. My method and your do return the same thing but it is definitely not the pixel coordinates of the point. – drheinrich940 Nov 12 '19 at 10:40
1

That is correct you have just got coordinates in the rasters coordinate system. I have updated my answer, to explain how you can convert that into pixels – the_cheff Nov 12 '19 at 11:15
you sir managed to eventually make me understand all this, thanks a lot this is perfect. I was missing some bricks in my understanding. – drheinrich940 Nov 13 '19 at 07:03
Hello, i am currently trying to perform the inverse operation based on the method you described, but the final coordinated are wrong. Could you take a look ? I could edit my initial post with a new section dedicated to pixel to long lat conversion if it is relevant. – drheinrich940 Nov 26 '19 at 10:30
In order to help other readers I think you should ask that as a new question, and post a link to it as a comment here. Then i will be happy to answer it ;-) – the_cheff Nov 26 '19 at 11:29
of course ! Here is the question https://stackoverflow.com/questions/59052516/find-lat-long-coordinates-from-pixel-point-in-geotiff-using-python-and-gdal – drheinrich940 Nov 26 '19 at 13:56
This code worked as is when using a dem in EPSG:4326, however to get it to work with a dem in EPSG:21781+5728, I had to flip `point.AddPoint(long, lat)` to `point.AddPoint(lat, long)`. No idea why... – fariadantes Apr 22 '23 at 18:09
1

@fariadantes In GDAL 3.0 special ordering was introduced for 4326. this issue is for instance mentioned here https://github.com/OSGeo/gdal/issues/1546 I do find it to be more confusing than helpful though. – the_cheff May 01 '23 at 12:27

score 2 · Answer 2 · edited Jul 04 '22 at 05:50

The proposed solution might work in most of the cases as row/column rotation is typically zero, but it should be at least checked or better included:

import numpy as np
from osgeo import gdal

def world_to_pxl(gt, x, y):
    # 'Affine transformation': W = A * pxl + ul
    # world[2, 1] = a[2, 2] * pxl[2, 1] + upper_left[2,1]
    world = np.array([[x], [y]])  # world coordinates
    upper_left = np.array(
        [
            [gt[0]], [gt[3]]    # upper left corner of image
        ]
    )
    a = np.array([[gt[1], gt[2]],
                  [gt[4], gt[5]]])
    
    # Reformulate: A^-1 * (W - ul) = pxl
    pxl = np.matmul(np.linalg.inv(a), (world - upper_left))
    row = pxl[0]  # x_pixel
    col = pxl[1]  # y_line
    return row, col

Find pixel coordinates from lat/long point in .geotiff using python and gdal

2 Answers2

Linked