使用python和gd从.geotiff中的lat/long点查找像素坐标

Question

我有（lat，long）坐标来描述一个点在.geotiff图像中的位置。

我希望找到图像中长的lat的等效像素坐标。

我成功地从命令行使用gdaltransform并执行以下指令：

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

但是我想从python代码中检索这种类型的等价性。我尝试了以下方法：

^{pr2}$

但它返回以下错误：

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)

我做错什么了？我试图克服这个错误，但没有成功。有别的办法吗？

编辑1:

我通过终端中的gdaltransform命令实现了一个单一的转换：

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

由于我需要以pythonic的方式检索像素，所以我尝试使用如下子进程调用命令：

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

但它不起作用。可能是因为命令本身的行为方式，比如不返回结果并结束自己，而是显示结果并保持忙碌。

Answer 1

根据the cookbook你正在翻转变换和点的用法。您应该在给定变换的点上调用transform，而不是相反。看起来你好像在翻转源和目标，但是你做了两次，所以它会起作用。在

但是我相信target.ImportFromUrl(path + TIFFFilename)不会起作用。相反，您可以使用gdal从geotiff中提取空间引用。在

下面这样的方法应该行得通

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引用中的坐标，但这不是像素坐标。在

要将点转换为像素，可以使用如下方法（线条的减号可能需要根据您所在的位置进行翻转）

^{pr2}$

所以你的最终代码看起来像

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)