Получение высоты на широте / долготе от растра с использованием Python?

10

Мне было интересно, есть ли у кого-то опыт получения данных высот из растра без использования ArcGIS , но лучше получить информацию в виде питона listили dict?

Я получаю свои данные XY в виде списка кортежей:

xy34 =[perp_obj[j].CalcPnts(float(i.dist), orientation) for j in range (len(perp_obj))]

Я хотел бы пройтись по списку или передать его в функцию или метод класса, чтобы получить соответствующее повышение для xy-пар.

Я провел некоторое исследование по этой теме, и API-интерфейс gdal звучит многообещающе. Может ли кто-нибудь посоветовать мне, как идти о вещах, подводных камнях, пример кода?


GDAL не вариант, так как я не могу редактировать системную переменную пути на машине, на которой я работаю!

Кто-нибудь знает о другом подходе?

LarsVegas
источник
2
к сожалению, вам действительно нужно запустить GDAL в вашей системе, чтобы что-то делать с растром в Python. С помощью «невозможно изменить системную переменную пути на компьютере», вы ссылаетесь на эти инструкции ? Я нахожу этот метод установки очень плохим, и я не использую его и не рекомендую его. Если вы используете Windows, установите GDAL / Python простым способом .
Майк Т
Да, я был действительно. Я сейчас не на работе, но проверю ссылку, которую вы разместили. Выглядит многообещающе! Спасибо, что вернулись к моему вопросу!
LarsVegas
Я использовал установщик Кристофа Гольке (ссылка выше) на многих рабочих компьютерах, и это действительно просто. Вам нужно только убедиться, что вы соответствуете версии Python и 32- или 64-битной Windows. Пока вы в этом, вы также должны получить NumPy из того же места, так как это необходимо GDAL, как показано в ответах ниже.
Майк Т

Ответы:

15

Вот более программный способ использования GDAL, чем ответ @ Aragon. Я не проверял это, но в основном это код, который работал для меня в прошлом. Он опирается на привязки Numpy и GDAL, но это все.

import osgeo.gdal as gdal
import osgeo.osr as osr
import numpy as np
from numpy import ma

def maFromGDAL(filename):
    dataset = gdal.Open(filename, gdal.GA_ReadOnly)

    if dataset is None:
        raise Exception()

    # Get the georeferencing metadata.
    # We don't need to know the CRS unless we want to specify coordinates
    # in a different CRS.
    #projection = dataset.GetProjection()
    geotransform = dataset.GetGeoTransform()

    # We need to know the geographic bounds and resolution of our dataset.
    if geotransform is None:
        dataset = None
        raise Exception()

    # Get the first band.
    band = dataset.GetRasterBand(1)
    # We need to nodata value for our MaskedArray later.
    nodata = band.GetNoDataValue()
    # Load the entire dataset into one numpy array.
    image = band.ReadAsArray(0, 0, band.XSize, band.YSize)
    # Close the dataset.
    dataset = None

    # Create a numpy MaskedArray from our regular numpy array.
    # If we want to be really clever, we could subclass MaskedArray to hold
    # our georeference metadata as well.
    # see here: http://docs.scipy.org/doc/numpy/user/basics.subclassing.html
    # For details.
    masked_image = ma.masked_values(image, nodata, copy=False)
    masked_image.fill_value = nodata

    return masked_image, geotransform

def pixelToMap(gt, pos):
    return (gt[0] + pos[0] * gt[1] + pos[1] * gt[2],
            gt[3] + pos[0] * gt[4] + pos[1] * gt[5])

# Reverses the operation of pixelToMap(), according to:
# https://en.wikipedia.org/wiki/World_file because GDAL's Affine GeoTransform
# uses the same values in the same order as an ESRI world file.
# See: http://www.gdal.org/gdal_datamodel.html
def mapToPixel(gt, pos):
    s = gt[0] * gt[4] - gt[3] * gt[1]
    x = (gt[4] * pos[0] - gt[1] * pos[1] + gt[1] * gt[5] - gt[4] * gt[2]) / s
    y = (-gt[3] * pos[0] + gt[0] * pos[1] + gt[3] * gt[2] - gt[0] * gt[5]) / s
    return (x, y)

def valueAtMapPos(image, gt, pos):
    pp = mapToPixel(gt, pos)
    x = int(pp[0])
    y = int(pp[1])

    if x < 0 or y < 0 or x >= image.shape[1] or y >= image.shape[0]:
        raise Exception()

    # Note how we reference the y column first. This is the way numpy arrays
    # work by default. But GDAL assumes x first.
    return image[y, x]

try:
    image, geotransform = maFromGDAL('myimage.tif')
    val = valueAtMapPos(image, geotransform, (434323.0, 2984745.0))
    print val
except:
    print('Something went wrong.')
MerseyViking
источник
1
см. редактирование на мой вопрос ... спасибо за публикацию в любом случае! Я проголосовал за это.
LarsVegas
1
Ах, черт! Ну, по крайней мере, это здесь для потомков. TBH, математика mapToPixel()и pixelToMap()важный бит, если вы можете создать простой массив (или обычный Python, но они обычно не так эффективны для такого рода вещей) и получить географическую ограничивающую рамку массива.
MerseyViking
1
+1 за комментарий (и код) об обращении параметров к массиву numpy. Я искал повсюду ошибку в моем коде, и этот обмен исправил ее!
Альдо
1
Тогда я предполагаю, что ваша матрица ( gtв примере) неверна. Аффинная матрица, используемая в CGAL (см. Gdal.org/gdal_datamodel.html ), обычно является обратимой (в противном случае у вас есть несколько весёлых значений масштабирования). Поэтому, где у нас есть, g = p.Aмы также можем сделать p = g.A^-1Numpy.linalg немного тяжелым для наших целей - мы можем свести все к двум простым уравнениям.
MerseyViking
1
Я заново отредактировал код, чтобы использовать простую алгебру, а не пустую строчку. Если математика неверна, исправьте страницу Википедии.
MerseyViking
3

Проверьте мой ответ здесь ... и прочитайте здесь для получения дополнительной информации. Следующая информация была взята из Geotips:

С gdallocationinfo мы можем запросить высоту в точке:

$ gdallocationinfo gmted/all075.vrt -geoloc 87360 19679

Вывод вышеуказанной команды имеет вид:

Report:
   Location: (87360P,19679L)
Band 1:
   Value: 1418

Это означает, что значение высоты при заданном геолокации составляет 1418.

Арагон
источник
Только что обнаружил, что не могу использовать GDAL, так как не могу редактировать системную переменную на машине, на которой я работаю. Спасибо за вклад, хотя.
LarsVegas
0

Посмотрите, например, этот код, который основан на GDAL (и Python, без нужды): https://github.com/geometalab/retrieve-height-service

Стефан
источник
К сожалению, код, похоже, не имеет лицензии на открытый код.
Бен Крауэлл
Теперь это имеет :-).
Стефан
-1

Предоставленный код Python извлекает данные о значении растровой ячейки на основе заданных координат x, y. Это слегка измененная версия примера из этого превосходного источника . Он основан GDALи numpyне является частью стандартного дистрибутива Python. Спасибо @Mike Toews за то, что они указали на неофициальные бинарные файлы Windows для пакетов расширения Python, чтобы сделать установку и использование быстрой и простой.

import os, sys, time, gdal
from gdalconst import *


# coordinates to get pixel values for
xValues = [122588.008]
yValues = [484475.146]

# set directory
os.chdir(r'D:\\temp\\AHN2_060')

# register all of the drivers
gdal.AllRegister()
# open the image
ds = gdal.Open('i25gn1_131.img', GA_ReadOnly)

if ds is None:
    print 'Could not open image'
    sys.exit(1)

# get image size
rows = ds.RasterYSize
cols = ds.RasterXSize
bands = ds.RasterCount

# get georeference info
transform = ds.GetGeoTransform()
xOrigin = transform[0]
yOrigin = transform[3]
pixelWidth = transform[1]
pixelHeight = transform[5]

# loop through the coordinates
for xValue,yValue in zip(xValues,yValues):
    # get x,y
    x = xValue
    y = yValue

    # compute pixel offset
    xOffset = int((x - xOrigin) / pixelWidth)
    yOffset = int((y - yOrigin) / pixelHeight)
    # create a string to print out
    s = "%s %s %s %s " % (x, y, xOffset, yOffset)

    # loop through the bands
    for i in xrange(1,bands):
        band = ds.GetRasterBand(i) # 1-based index
        # read data and add the value to the string
        data = band.ReadAsArray(xOffset, yOffset, 1, 1)
        value = data[0,0]
        s = "%s%s " % (s, value) 
    # print out the data string
    print s
    # figure out how long the script took to run
LarsVegas
источник
Похоже, что это просто менее универсальная, менее гибкая версия того, что MerseyViking предлагал выше?
WileyB