Получить широту и долготу из файла GeoTIFF - PullRequest
Новая
       19

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

37 голосов
/ 27 мая 2010

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

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

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

Какая процедура требуется для этого?

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

Ответы [ 2 ]

76 голосов
/ 28 мая 2010

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

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, как в примере выше, вы имеете дело с спроецированной системой координат. Спроецированная система координат представляет собой представление сфероидальной земной поверхности, но сглажено и искажено на плоскости . Если вам нужны широта и долгота, вам нужно преобразовать координаты в нужную вам географическую систему координат.

К сожалению, не все пары широта / долгота созданы равными, основываясь на различных сфероидальных моделях Земли. В этом примере я конвертирую в WGS84 , географическую систему координат, предпочитаемую в GPS и используемую всеми популярными веб-сайтами. Система координат определяется четко определенной строкой. Их каталог доступен с пространственная ссылка , см., Например, 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)

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

14 голосов
/ 27 мая 2010

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

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

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

...