Obtendo elevação em lat / long de varredura usando python?

10

Eu queria saber se alguém tem alguma experiência em obter dados de elevação de um raster sem usar o ArcGIS , mas obter as informações como um python listou dict?

Eu recebo meus dados XY como uma lista de tuplas:

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

Gostaria de percorrer a lista ou passá-la para uma função ou método de classe para obter a elevação correspondente para os pares xy.

Eu fiz algumas pesquisas sobre o assunto e a API gdal parece promissora. Alguém pode me aconselhar sobre como fazer coisas, armadilhas, código de exemplo?


GDAL não é uma opção, pois não consigo editar a variável de caminho do sistema na máquina em que estou trabalhando!

Alguém sabe sobre uma abordagem diferente?

LarsVegas
fonte
2
infelizmente, você realmente precisa executar o GDAL em seu sistema para fazer qualquer coisa com uma varredura em Python. Com "não é possível editar a variável do caminho do sistema na máquina", você está se referindo a estas instruções ? Acho esse método de instalação muito ruim e não o uso nem o recomendo. Se você estiver usando o Windows, instale o GDAL / Python da maneira mais simples .
Mike T
Sim, eu estava mesmo. Não estou trabalhando agora, mas vou verificar o link que você postou. Parece promissor! Obrigado por voltar à minha pergunta!
LarsVegas
Eu usei o instalador de Christoph Gohlke (link acima) em muitos computadores de trabalho, e é realmente simples. Você só precisa garantir a correspondência da versão do Python e do Windows de 32 ou 64 bits. Enquanto você está nisso, você também deve obter o NumPy do mesmo local, pois isso é necessário para a GDAL, como mostrado nas respostas abaixo.
23712 Mike

Respostas:

15

Aqui está uma maneira mais programática de usar o GDAL do que a resposta da @ Aragon. Eu não testei, mas é principalmente o código da placa da caldeira que funcionou para mim no passado. Ele se baseia nas ligações Numpy e GDAL, mas é isso.

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
fonte
1
veja a edição da minha pergunta ... obrigado por postar de qualquer maneira! Eu votei nele.
LarsVegas
1
Ah droga! Bem, pelo menos, está aqui para a posteridade. TBH, a matemática em mapToPixel()e pixelToMap()são pouco importante, contanto que você pode criar uma matriz numpy (ou um Python regular de um, mas eles geralmente não são tão eficientes para este tipo de coisa), e obter caixa delimitadora geográfica da matriz.
MerseyViking
1
+1 no comentário (e código) sobre a reversão dos parâmetros para a matriz numpy. Eu estava procurando em todo lugar um bug no meu código, e essa troca o corrigiu!
aldo
1
Então eu sugiro que sua matriz ( gtno exemplo) esteja errada. Uma matriz afim conforme usada em CGAL (consulte: gdal.org/gdal_datamodel.html ) geralmente é invertível (caso contrário, você tem alguns valores de escala descolados). Então, onde temos g = p.A, também podemos fazer o p = g.A^-1Numpy.linalg é um pouco pesado para nossos propósitos - podemos reduzir tudo para duas equações simples.
MerseyViking
1
Reeditei o código para usar álgebra simples em vez de numpy linalg. Se a matemática estiver errada, corrija a página da Wikipedia.
MerseyViking
3

Confira minha resposta aqui ... e leia aqui para obter mais informações. As seguintes informações foram obtidas das geotips:

Com gdallocationinfo , podemos consultar a elevação em um ponto:

$ gdallocationinfo gmted/all075.vrt -geoloc 87360 19679

A saída do comando acima tem a forma:

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

Isso significa que o valor da elevação na geolocalização fornecida é 1418.

Aragão
fonte
Acabei de descobrir que não posso usar o GDAL, pois não consigo editar minha variável de sistema na máquina em que estou trabalhando. Obrigado pela contribuição embora.
LarsVegas
0

Veja, por exemplo, este código baseado em GDAL (e Python, não é necessário numpy): https://github.com/geometalab/retrieve-height-service

Stefan
fonte
É lamentável que o código não pareça ser licenciado em código aberto.
Ben Crowell
Agora tem :-).
Stefan
-1

O código python fornecido extrai os dados de valor de uma célula de varredura com base em determinadas cordas x, y. É uma versão ligeiramente alterada de um exemplo desta excelente fonte . Ele é baseado GDALe numpynão faz parte da distribuição python padrão. Agradecemos ao @ Mike Toews por apontar os binários não oficiais do Windows para pacotes de extensão Python para tornar a instalação e o uso rápidos e fáceis.

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
fonte
Parece que esta é apenas uma versão menos genérica e menos flexível do que o MerseyViking ofereceu acima?
WileyB 17/02