¿Obteniendo elevación en lat / long desde la trama usando python?

10

Me preguntaba si alguien tiene alguna experiencia en obtener datos de elevación de un ráster sin usar ArcGIS , sino obtener la información como un pitón listo dict?

Recibo mis datos XY como una lista de tuplas:

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

Me gustaría recorrer la lista o pasarla a una función o método de clase para obtener la elevación correspondiente para los pares xy.

Investigué un poco sobre el tema y la API de gdal parece prometedora. ¿Alguien puede aconsejarme cómo hacer cosas, trampas, código de muestra?


¡GDAL no es una opción ya que no puedo editar la variable de ruta del sistema en la máquina en la que estoy trabajando!

¿Alguien sabe acerca de un enfoque diferente?

LarsVegas
fuente
2
desafortunadamente, realmente necesita que GDAL se ejecute en su sistema para hacer algo con un ráster en Python. Con "no se puede editar la variable de ruta del sistema en la máquina", ¿se refiere a estas instrucciones ? Este método de instalación me parece muy pobre, y no lo uso ni lo recomiendo. Si está utilizando Windows, instale GDAL / Python de la manera más sencilla .
Mike T
Sí, lo estaba de hecho. No estoy en el trabajo en este momento, pero revisaré el enlace que publicaste. ¡Parece prometedor! ¡Gracias por volver a mi pregunta!
LarsVegas
He utilizado el instalador de Christoph Gohlke (vinculado anteriormente) en muchas computadoras de trabajo, y es realmente simple. Solo necesita asegurarse de que coincida con la versión de Python y con Windows de 32 o 64 bits. Mientras lo hace, también debe obtener NumPy del mismo lugar, ya que GDAL lo necesita, como se muestra en las respuestas a continuación.
Mike T

Respuestas:

15

Aquí hay una forma más programática de usar GDAL que la respuesta de @ Aragon. No lo he probado, pero en su mayoría es un código de placa de caldera que me ha funcionado en el pasado. Se basa en enlaces Numpy y GDAL, pero eso es todo.

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
fuente
1
vea la edición de mi pregunta ... ¡gracias por publicar de todos modos! Yo lo voté.
LarsVegas
1
Ah maldición! Bueno, al menos está aquí para la posteridad. TBH, las matemáticas en mapToPixel()y pixelToMap()son la parte importante, siempre y cuando se puede crear una matriz numpy (o una pitón uno normal, pero por lo general no son tan eficientes para este tipo de cosas), y obtener cuadro de delimitación geográfica de la matriz.
MerseyViking
1
+1 para el comentario (y el código) sobre la inversión de los parámetros a la matriz numpy. Estaba buscando en todas partes un error en mi código, ¡y este intercambio lo solucionó!
aldo
1
Entonces sugiero que su matriz ( gten el ejemplo) esté equivocada. Una matriz afín como se usa en CGAL (ver: gdal.org/gdal_datamodel.html ) es generalmente invertible (de lo contrario, tiene algunos valores de escala funky en curso). Entonces, donde tenemos g = p.A, también podemos hacer p = g.A^-1Numpy.linalg es un poco pesado para nuestros propósitos: podemos reducir todo a dos ecuaciones simples.
MerseyViking
1
He reeditado el código para usar álgebra simple en lugar de linalg numpy. Si las matemáticas están mal, arregle la página de Wikipedia.
MerseyViking
3

Mira mi respuesta aquí ... y lee aquí para obtener más información. La siguiente información fue tomada de Geotips:

Con gdallocationinfo , podemos consultar la elevación en un punto:

$ gdallocationinfo gmted/all075.vrt -geoloc 87360 19679

La salida del comando anterior tiene la forma:

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

Esto significa que el valor de elevación en la geolocalización proporcionada es 1418.

Aragón
fuente
Me acabo de enterar que no puedo usar GDAL ya que no puedo editar mi variable de sistema en la máquina en la que estoy trabajando. Gracias por la entrada sin embargo.
LarsVegas
0

Consulte, por ejemplo, este código que se basa en GDAL (y Python, no se necesita numpy): https://github.com/geometalab/retrieve-height-service

Stefan
fuente
Es lamentable que el código no parezca ser de código abierto con licencia.
Ben Crowell
Ahora tiene :-).
Stefan
-1

El código python proporcionado extrae los datos de valor de una celda ráster en función de las coordenadas x, y. Es una versión ligeramente alterada de un ejemplo de esta excelente fuente . Se basa GDALy numpyno forma parte de la distribución estándar de Python. Gracias a @Mike Toews por señalar los binarios no oficiales de Windows para los paquetes de extensión de Python para hacer que la instalación y el uso sean rápidos y fáciles.

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
fuente
Parece que esta es solo una versión menos genérica y menos flexible de lo que MerseyViking ofreció anteriormente.
WileyB