Найти координаты пикселей от широты / долготы в .geotiff, используя python и gdal - PullRequest
Новая
       145

Найти координаты пикселей от широты / долготы в .geotiff, используя python и gdal

0 голосов
/ 30 октября 2019

У меня есть координаты (широта, долгота), описывающие положение точки на изображении .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

Поскольку мне нужно получить пиксель вПитонически, я пытался вызвать команду с помощью подпроцесса, как: 

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

Но это не работает. Возможно, из-за того, как ведет себя сама команда, например, не возвращает результат и не завершает себя, а отображает результат и остается занятым.

1 Ответ

1 голос
/ 12 ноября 2019

Согласно поваренной книге вы переворачиваете использование преобразования и точки. Вы должны вызывать transform в точке, заданной для transform, а не наоборот. Также кажется, что вы переключаете источник и цель, но вы делаете это два раза, так что это будет работать.

Однако я считаю, что target.ImportFromUrl(path + TIFFFilename) не будет работать. Вместо этого вы можете извлечь пространственную привязку из геотиффа, используя gdal.

Должно работать что-то вроде следующего 

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

Это дает вам координаты в вашей привязке геотипов, однако это не пикселькоординаты.

Чтобы преобразовать точку в пиксели, вы можете использовать что-то вроде следующего (минус для линии, возможно, придется перевернуть, основываясь на том, где вы находитесь в мире) 

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