Obtenir l'élévation lat / long du raster en utilisant python?

10

Je me demandais si quelqu'un avait une certaine expérience dans l'obtention de données d'altitude à partir d'un raster sans utiliser ArcGIS , mais plutôt obtenir les informations sous forme de python listou dict?

J'obtiens mes données XY sous forme de liste de tuples:

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

Je voudrais parcourir la liste ou la passer à une fonction ou à une méthode de classe pour obtenir l'élévation correspondante pour les paires xy.

J'ai fait quelques recherches sur le sujet et l' API gdal semble prometteuse. Quelqu'un peut-il me conseiller sur la façon de procéder, les pièges, un exemple de code?


GDAL n'est pas une option car je ne peux pas modifier la variable de chemin d'accès système sur la machine sur laquelle je travaille!

Quelqu'un connaît-il une approche différente?

LarsVegas
la source
2
malheureusement, vous avez vraiment besoin de faire fonctionner GDAL sur votre système pour faire quoi que ce soit avec un raster en Python. Avec "impossible de modifier la variable de chemin d'accès système sur la machine", faites-vous référence à ces instructions ? Je trouve cette méthode d'installation très médiocre et je ne l'utilise pas et ne la recommande pas. Si vous utilisez Windows, installez GDAL / Python de manière simple .
Mike T
Oui, je l'étais en effet. Je ne suis pas au travail en ce moment, mais je vais consulter le lien que vous avez publié. Cela semble prometteur! Merci d'être revenu à ma question!
LarsVegas
J'ai utilisé l'installateur de Christoph Gohlke (lié ci-dessus) sur de nombreux ordinateurs de travail, et c'est vraiment simple. Vous devez uniquement vous assurer que vous correspondez à la version de Python et à Windows 32 ou 64 bits. Pendant que vous y êtes, vous devez également obtenir NumPy au même endroit, car cela est nécessaire à GDAL, comme indiqué dans les réponses ci-dessous.
Mike T

Réponses:

15

Voici une manière plus programmatique d'utiliser GDAL que la réponse de @ Aragon. Je ne l'ai pas testé, mais c'est surtout le code de plaque de chaudière qui a fonctionné pour moi dans le passé. Il s'appuie sur les liaisons Numpy et GDAL, mais c'est tout.

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
la source
1
voir la modification de ma question ... merci de poster quand même! Je l'ai voté.
LarsVegas
1
Ah putain! Eh bien au moins c'est ici pour la postérité. TBH, les calculs mapToPixel()et pixelToMap()sont le bit important, tant que vous pouvez créer un tableau numpy (ou un tableau Python normal, mais ils ne sont généralement pas aussi efficaces pour ce genre de chose), et obtenir le cadre de délimitation géographique du tableau.
MerseyViking
1
+1 pour le commentaire (et le code) sur l'inversion des paramètres dans le tableau numpy. Je cherchais partout un bogue dans mon code, et cet échange l'a corrigé!
aldo
1
Ensuite, je suggère que votre matrice ( gtdans l'exemple) est fausse. Une matrice affine telle qu'utilisée dans CGAL (voir: gdal.org/gdal_datamodel.html ) est généralement inversible (sinon vous avez des valeurs d'échelle funky en cours). Donc, là où nous avons, g = p.Anous pouvons également faire p = g.A^-1Numpy.linalg est un peu lourd pour nos besoins - nous pouvons tout réduire à deux équations simples.
MerseyViking
1
J'ai réédité le code pour utiliser l'algèbre simple plutôt que linalg numpy. Si les calculs sont erronés, corrigez la page Wikipedia.
MerseyViking
3

Consultez ma réponse ici ... et lisez ici pour quelques informations. Les informations suivantes proviennent de Geotips:

Avec gdallocationinfo , nous pouvons interroger l'élévation en un point:

$ gdallocationinfo gmted/all075.vrt -geoloc 87360 19679

La sortie de la commande ci-dessus a la forme:

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

Cela signifie que la valeur d'élévation à la géolocalisation fournie est 1418.

Aragon
la source
Je viens de découvrir que je ne peux pas utiliser GDAL car je ne suis pas en mesure de modifier ma variable système sur la machine sur laquelle je travaille. Merci pour votre contribution.
LarsVegas
0

Voir par exemple ce code basé sur GDAL (et Python, aucun numpy nécessaire): https://github.com/geometalab/retrieve-height-service

Stefan
la source
Il est regrettable que le code ne semble pas être sous licence open source.
Ben Crowell
Maintenant, il a :-).
Stefan
-1

Le code python fourni extrait les données de valeur d'une cellule raster sur la base de coordonnées x, y données. Il s'agit d'une version légèrement modifiée d'un exemple de cette excellente source . Il est basé sur GDALet numpyqui ne fait pas partie de la distribution standard de python. Merci à @Mike Toews d'avoir signalé les fichiers binaires non officiels de Windows pour les extensions Python pour rendre l'installation et l'utilisation rapides et faciles.

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
la source
Il semble que ce soit juste une version moins générique et moins flexible de ce que MerseyViking a offert ci-dessus?
WileyB