Как мне преобразовать аффинные координаты в широту / долготу?

Я чрезвычайно новичок в ГИС.

Я использую, gdalчтобы читать на карте землепользования / земельного покрова, и мне нужно выбрать широту / долготу определенных типов земного покрова, чтобы индексировать их в другой набор данных, который выражается только в широте / долготе. К сожалению, я не понимаю форму координат x и y, данных мне от геотрансформации, в частности, originXи originYниже:

geotransform = dataset.GetGeoTransform()
originX = geotransform[0]
originY = geotransform[3]

Печать этих значений дает мне координаты как (447466.693808, 4952570.40529). Как они связаны с первоначальной широтой и долготой?

Редактировать:

Вот простой пример Python, который дал мне то, что я искал:

srs = osr.SpatialReference()
srs.ImportFromWkt(dataset.GetProjection())

srsLatLong = srs.CloneGeogCS()
ct = osr.CoordinateTransformation(srs,srsLatLong)
print ct.TransformPoint(originX,originY)

Украдено из: tolatlong.py

gdal_translate перепроектирует ваши данные из любой проекции в любую другую (в этом случае вы хотите EPSG: 4326), используя:

gdal_translate -a_srs epsg:4326 srcfile new_file 

или вы можете использовать gdaltrasform для преобразования точек (и я уверен, что вы можете получить к нему доступ и из Python (?) тоже)

Геотрансформация документирована по адресу https://gdal.org/user/raster_data_model.html . Идея состоит в том, что вы берете (x, y) координаты прямо из набора данных, применяете линейное преобразование, чтобы получить (u, v) с

u = a*x + b*y
v = c*x + d*y

(вы можете принять это за определение линейного преобразования), затем сдвинуть источник, добавив геотрансформацию [0] к u и геотрансформирование [3] к v. Это дает «аффинное преобразование» (x, y). Он действительно предназначен для поворота, изменения масштаба, возможно, немного исправления для некоторых ошибок перекоса и изменения положения координат данных (x, y) в соответствии с известной системой координат. Предполагается, что в результате получаются спроецированные координаты. Это просто означает, что существует математическая процедура, принимающая (долгота, широта) и превращающая их в известные координаты: это называется «проекцией». «Непроектирование» делает обратное; поэтому, если вы знаете, какая проекция необходима, вы применяете ее к аффинно преобразованным (x, y) координатам получить широту и долготу.

Кстати, значения констант a, b, c, d задаются записями 1, 2, 4 и 5 в массиве геотрансформации.

Вы можете использовать следующее:

import gdal, numpy as n
def coord(file):
    padfTransform = file.GetGeoTransform()
    indices = n.indices(file.ReadAsArray().shape)
    xp = padfTransform[0] + indices[1]*padfTransform[1] + indices[1]*padfTransform[2]   
    yp = padfTransform[3] + indices[0]*padfTransform[4] + indices[0]*padfTransform[5]  
    return xp,yp
file = gdal.Open('Yourfile.tif') # A GeoTiff file
x,y = coord(file)

Координата вернет долготу (х) и широту (у) всех пикселей. Имейте в виду, что координаты левого угла пикселя