0

我有（纬度，经度）坐标来描述 .geotiff 图像中一个点的位置。

我希望找到图像内经纬度的等效像素坐标。

我通过以下指令从命令行使用 gdaltransform 成功：

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

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

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)

我究竟做错了什么 ？我试图解决这个错误，但没有任何成功。还有其他方法吗？

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

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

4

2 回答 2

5

根据食谱，您正在翻转使用变换和点。您应该在给定变换的点上调用变换，而不是相反。似乎您正在翻转源和目标，但是您这样做了两次，所以它会起作用。

但是我相信这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 参考中的坐标，但这不是像素坐标。

要将点转换为像素，您可以使用以下内容（根据您在世界上的位置，可能必须翻转线的减号）

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

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

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)
于 2019-11-12T08:17:33.113 回答
0

建议的解决方案可能在大多数情况下都有效，因为行/列旋转通常为零，但至少应该检查或更好地包括：

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
于 2021-03-05T11:41:19.650 回答