Получение широты и долготы из файла GeoTIFF

Question

Использование GDAL в Python, как вы можете получить широту и долготу файла GeoTIFF?

GeoTIFF не хранит никакой информации о координатах. Вместо этого они сохраняют координаты XY Origin. Однако координаты XY не обеспечивают широту и долготу верхнего левого угла и нижнего левого угла.

Похоже, что для решения этой проблемы мне понадобится математика, но я не знаю, с чего начать.

Какая процедура должна быть выполнена?

Я знаю, что метод GetGeoTransform() важен для этого, однако я не знаю, что с ним делать.

Answer 1

Чтобы получить координаты углов вашего GeoTIFF сделайте следующее:

from osgeo import gdal 
ds = gdal.Open('path/to/file') 
width = ds.RasterXSize 
height = ds.RasterYSize 
gt = ds.GetGeoTransform() 
minx = gt[0] 
miny = gt[3] + width*gt[4] + height*gt[5] 
maxx = gt[0] + width*gt[1] + height*gt[2] 
maxy = gt[3] 

Однако, это не может быть в формате широта/долгота. Как отметил Джастин, ваш геотиф будет храниться с какой-то системой координат. Если вы не знаете, что система координат это, вы можете узнать, запустив gdalinfo:

gdalinfo ~/somedir/somefile.tif 

Какие выходы:

Driver: GTiff/GeoTIFF 
Size is 512, 512 
Coordinate System is: 
PROJCS["NAD27/UTM zone 11N", 
    GEOGCS["NAD27", 
     DATUM["North_American_Datum_1927", 
      SPHEROID["Clarke 1866",6378206.4,294.978698213901]], 
     PRIMEM["Greenwich",0], 
     UNIT["degree",0.0174532925199433]], 
    PROJECTION["Transverse_Mercator"], 
    PARAMETER["latitude_of_origin",0], 
    PARAMETER["central_meridian",-117], 
    PARAMETER["scale_factor",0.9996], 
    PARAMETER["false_easting",500000], 
    PARAMETER["false_northing",0], 
    UNIT["metre",1]] 
Origin = (440720.000000,3751320.000000) 
Pixel Size = (60.000000,-60.000000) 
Corner Coordinates: 
Upper Left ( 440720.000, 3751320.000) (117d38'28.21"W, 33d54'8.47"N) 
Lower Left ( 440720.000, 3720600.000) (117d38'20.79"W, 33d37'31.04"N) 
Upper Right ( 471440.000, 3751320.000) (117d18'32.07"W, 33d54'13.08"N) 
Lower Right ( 471440.000, 3720600.000) (117d18'28.50"W, 33d37'35.61"N) 
Center  ( 456080.000, 3735960.000) (117d28'27.39"W, 33d45'52.46"N) 
Band 1 Block=512x16 Type=Byte, ColorInterp=Gray 

Этот вывод может быть все, что вам нужно. Если вы хотите сделать это программно в python, так это то, как вы получаете ту же информацию.

Если система координат является PROJCS, как пример выше, вы имеете дело с проецируемой системой координат. Проецируемая система координат - это representation of the spheroidal earth's surface, but flattened and distorted onto a plane. Если вам нужна широта и долгота, вам необходимо преобразовать координаты в нужную географическую систему координат.

К сожалению, не все пары широты/долготы созданы равными, основываясь на разных сфероидальных моделях Земли. В этом примере я конвертирую в WGS84, географическую систему координат, предпочитаемую в GPS и используемую всеми популярными веб-сайтами картирования. Система координат определяется четко определенной строкой. Каталог из них можно получить от spatial ref, см., Например, WGS84.

from osgeo import osr, gdal 

# get the existing coordinate system 
ds = gdal.Open('path/to/file') 
old_cs= osr.SpatialReference() 
old_cs.ImportFromWkt(ds.GetProjectionRef()) 

# create the new coordinate system 
wgs84_wkt = """ 
GEOGCS["WGS 84", 
    DATUM["WGS_1984", 
     SPHEROID["WGS 84",6378137,298.257223563, 
      AUTHORITY["EPSG","7030"]], 
     AUTHORITY["EPSG","6326"]], 
    PRIMEM["Greenwich",0, 
     AUTHORITY["EPSG","8901"]], 
    UNIT["degree",0.01745329251994328, 
     AUTHORITY["EPSG","9122"]], 
    AUTHORITY["EPSG","4326"]]""" 
new_cs = osr.SpatialReference() 
new_cs .ImportFromWkt(wgs84_wkt) 

# create a transform object to convert between coordinate systems 
transform = osr.CoordinateTransformation(old_cs,new_cs) 

#get the point to transform, pixel (0,0) in this case 
width = ds.RasterXSize 
height = ds.RasterYSize 
gt = ds.GetGeoTransform() 
minx = gt[0] 
miny = gt[3] + width*gt[4] + height*gt[5] 

#get the coordinates in lat long 
latlong = transform.TransformPoint(x,y) 

Надеюсь, это сделает то, что вы хотите.

Answer 2

Я не знаю, если это полный ответ, но this site говорит:

X/Y размеры карты называются пищеблок и Northing. Для наборов данных в географической системе координат они будут содержать долготу и широту. Для проектируемых систем координат они обычно были бы восточным и северным в проектируемой системе координат. Для изображений без привязки восток и север должны быть просто смещением пикселов/линий каждого пикселя (что подразумевается единственной геотрансформацией).

поэтому они могут быть на самом деле долготой и широтой.