使用python和gdal从.geotiff中的纬度/经度点查找像素坐标 [英] Find pixel coordinates from lat/long point in .geotiff using python and gdal

问题描述

我有(纬长)坐标来描述.geotiff图像中点的位置.

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.

我在命令行中使用gdaltransform成功执行了以下指令:

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

但是我想从python代码中检索这种类型的等效项.我尝试了以下方法:

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

但是它返回此错误:

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 ?

我通过终端中的gdaltransform命令实现了一次转换:

I achieved a single transformation via gdaltransform commands in terminal :

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

由于我需要以pythonic方式检索像素，因此我尝试使用像这样的子进程来调用命令:

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.

推荐答案

根据食谱，您正在翻转使用transform和point.您应该在给定转换的点上调用转换，而不是相反.似乎您也在翻转源和目标，但是您做了两次，所以它可以工作.

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.

但是，我相信target.ImportFromUrl(path + TIFFFilename)无法正常工作.相反，您可以使用gdal从Geotiff中提取空间参考.

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

这将为您提供Geotiffs参考中的坐标，但这不是像素坐标.

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)

这篇关于使用python和gdal从.geotiff中的纬度/经度点查找像素坐标的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！