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

Question

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

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

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

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

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

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

40

python math geolocation tiff gdal

задан Phanto 27 May 2010 в 15:36

2 ответа

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

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

, так что на самом деле они могут быть долготой и широтой.

16

ответ дан 27 November 2019 в 01:21

Другие вопросы по тегам:

python math geolocation tiff gdal

Похожие вопросы:

score 80 · Accepted Answer

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

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 и на всех популярных картографических сайтах. Система координат определяется четко определенной строкой. Их каталог можно найти на сайте 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(minx,miny)

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