LINUX.ORG.RU
Форум — Development  

GDAL: как из содержимого Geotiffки получить XYZ-координаты?

 ,


0

1

I want to get an XYZ-like surface from a Geotiff with GDAL in Python. Собственно, делаю я вот что:

1) Качаю Geotiff отсюда: http://srtm.csi.cgiar.org/SELECTION/inputCoord.asp

2) Пишу вот такой код:

from osgeo import ogr, osr, gdal

file_path = os.path.join(os.path.dirname("/home/yak/dump/srtm_62_05.tif")
gdal_data = gdal.Open(file_path, gdal.GA_ReadOnly)
raster_band = gdal_data.GetRasterBand(1)
raster = raster_band.ReadAsArray()

Всё круто, только, во-первых, ни разу не понятно, в каких попугаях измеряется возвышение (написано, что в метрах, но относительно каких единиц? Да и в тех местах, где море, стоят какие-то страшные -32000 с гаком) — а во-вторых, мне бы это сконвертировать в обычные декартовы XYZ-координаты.

Собственно, имею спросить: а как бы это проделать-то? Явно же есть какие-то стандартные средства в самом GDAL?

★

Последнее исправление: Yak (всего исправлений: 2)

Давай я погуглю за тебя: http://google.com/search?q=qet+xyz+from+geotiff

https://gis.stackexchange.com/questions/220561/how-to-get-long-lat-and-altitu...

https://stackoverflow.com/questions/2922532/obtain-latitude-and-longitude-fro...

https://gis.stackexchange.com/questions/255596/gdal-for-python-how-can-i-get-...

https://gis.stackexchange.com/questions/217890/converting-tif-to-xyz-in-qgis

https://gis.stackexchange.com/questions/59323/convert-geotiff-to-simple-xyz-e...

atsym ★★★★★
()
Ответ на: Давай я погуглю за тебя: http://google.com/search?q=qet+xyz+from+geotiff от atsym

Да я гуглил уже, но там в половине случаев получают [широту, долготу, возвышение], причём возвышение от какой-то странной отметки, а не уровнем-моря. А в другой половине случаев конвертят не через библиотеку GDAL, а с помощью утилит каких-то.

Yak ★
() автор топика

стоят какие-то страшные -32000 с гаком

это nodata-значения.

Берешь программку gdalinfo, натравливаешь на свой файл, там выясняешь что srtm dem создан в системе координат UTM (нагуглишь сам дальше что такое UTM).

У твоего файлика, кроме датума есть еще и transformation (6 double коэфициентов афинного преобразования), который как раз и указывает как положение пикселя в растре перевести в значения в UTM, для этого смотрим в доку по GDAL и находим (http://www.gdal.org/classGDALDataset.html#a5101119705f5fa2bc1344ab26f66fd1d): 

Fetches the coefficients for transforming between pixel/line (P,L) raster space, and projection coordinates (Xp,Yp) space.

Xp = padfTransform[0] + P*padfTransform[1] + L*padfTransform[2];
Yp = padfTransform[3] + P*padfTransform[4] + L*padfTransform[5];

Вообщем положение пикселя даст тебе x, y, а значение интенсивности - z.

Если значение Z совпадает с значением nodata - значит пропускаешь этот пиксель.

Загугли gdal api tutorial, там все есть.

nikitos ★★★
()

А ещё можно взять gdal2xyz и не изобретать велосипед. На входе растр, на выходе CSV c искомыми значениями.

Единицы измерения высоты зависят от системы координат исходного растра. Если растр в географических координатах (широта-долгота), то предварительно его надо трансформировать либо в псевдо-меркатор, либо в местную СК. Трасформацию можно провести либо через gdalwarp, либо через API.

P.S.: лучше не использовать ReadAsArray(). Конечно, если у вас «игрушечные» растры малого размера, то ничего страшного. А при работе с настоящими данными легко получите MemoryError, более-менее приличный растр целиком в память, скорее всего, не поместится. Читайте блоками, оптимальный размер блока есть в метаданных растра.

voltron ★
()
Ответ на: комментарий от voltron

Ну, у меня растры 6001x6001, там не всё так страшно.

Yak ★
() автор топика
Вы не можете добавлять комментарии в эту тему. Тема перемещена в архив.
Development