Mendapatkan elevasi pada lat / long dari raster menggunakan python?

10

Saya bertanya-tanya apakah ada yang punya pengalaman dalam mendapatkan data elevasi dari raster tanpa menggunakan ArcGIS , tetapi lebih suka mendapatkan informasi sebagai python listatau dict?

Saya mendapatkan data XY saya sebagai daftar tupel:

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

Saya ingin mengulang daftar atau meneruskannya ke fungsi atau metode-kelas untuk mendapatkan elevasi yang sesuai untuk pasangan xy.

Saya melakukan riset pada topik dan API gdal terdengar menjanjikan. Adakah yang bisa menasihati saya cara mengerjakan sesuatu, perangkap, kode sampel?


GDAL bukan opsi karena saya tidak dapat mengedit variabel jalur sistem pada mesin yang sedang saya kerjakan!

Adakah yang tahu tentang pendekatan yang berbeda?

LarsVegas
sumber
2
Sayangnya, Anda benar-benar perlu menjalankan GDAL di sistem Anda untuk melakukan apa saja dengan raster di Python. Dengan "tidak dapat mengedit variabel jalur sistem pada mesin", apakah Anda merujuk pada instruksi ini ? Saya menemukan metode instalasi ini sangat buruk, dan saya tidak menggunakannya atau merekomendasikannya. Jika Anda menggunakan Windows, instal GDAL / Python dengan cara sederhana .
Mike T
Ya, memang. Saya tidak sedang bekerja sekarang tetapi saya akan memeriksa tautan yang Anda poskan. Terlihat menjanjikan! Terima kasih telah kembali ke pertanyaan saya!
LarsVegas
Saya telah menggunakan installer oleh Christoph Gohlke (ditautkan di atas) pada banyak komputer kerja, dan itu sangat sederhana. Anda hanya perlu memastikan bahwa Anda cocok dengan versi Python dan Windows 32- atau 64-bit. Ketika Anda berada di sana, Anda juga harus mendapatkan NumPy dari tempat yang sama, karena itu dibutuhkan oleh GDAL, seperti yang ditunjukkan dalam jawaban di bawah ini.
Mike T

Jawaban:

15

Inilah cara yang lebih terprogram untuk menggunakan GDAL daripada jawaban @ Aragon. Saya belum mengujinya, tetapi sebagian besar kode boiler-plate yang telah bekerja untuk saya di masa lalu. Itu bergantung pada ikatan Numpy dan GDAL, tapi itu saja.

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
sumber
1
lihat hasil edit pertanyaan saya ... terima kasih sudah memposting! Saya membatalkannya.
LarsVegas
1
Ah sial! Yah setidaknya itu di sini untuk anak cucu. TBH, matematika di dalam mapToPixel()dan pixelToMap()merupakan bit yang penting, selama Anda dapat membuat array numpy (atau yang Python biasa, tetapi mereka umumnya tidak seefisien untuk hal-hal semacam ini), dan dapatkan kotak batas geografis array.
MerseyViking
1
+1 untuk komentar (dan kode) tentang membalikkan parameter ke array numpy. Saya mencari bug di mana-mana di kode saya, dan swap ini memperbaikinya!
aldo
1
Maka saya sarankan matriks Anda ( gtdalam contoh) salah. Matriks affine seperti yang digunakan dalam CGAL (lihat: gdal.org/gdal_datamodel.html ) umumnya tidak dapat dibalik (jika tidak, Anda memiliki beberapa nilai penskalaan yang funky terjadi). Jadi di mana kita miliki, g = p.Akita juga dapat melakukan p = g.A^-1Numpy.linalg agak berat untuk tujuan kita - kita dapat mengurangi segalanya menjadi dua persamaan sederhana.
MerseyViking
1
Saya telah mengedit ulang kode untuk menggunakan aljabar sederhana daripada linalg numpy. Jika matematika salah, perbaiki halaman Wikipedia.
MerseyViking
3

Lihat jawaban saya di sini ... dan baca di sini untuk beberapa informasi. Info berikut diambil dari Geotips:

Dengan gdallocationinfo , kami dapat menanyakan ketinggian pada suatu titik:

$ gdallocationinfo gmted/all075.vrt -geoloc 87360 19679

Output dari perintah di atas memiliki bentuk:

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

Ini berarti, bahwa nilai ketinggian pada geolokasi yang disediakan adalah 1418.

Aragon
sumber
Baru diketahui saya tidak dapat menggunakan GDAL karena saya tidak dapat mengedit variabel sistem pada mesin yang sedang saya kerjakan. Terima kasih atas masukannya.
LarsVegas
0

Lihat misalnya kode ini yang didasarkan pada GDAL (dan Python, tidak perlu numpy): https://github.com/geometalab/retrieve-height-service

Stefan
sumber
Sangat disayangkan bahwa kode tersebut tampaknya tidak berlisensi sumber terbuka.
Ben Crowell
Sekarang memiliki :-).
Stefan
-1

Kode python yang disediakan mengekstrak data nilai sel raster berdasarkan pada x, y yang diberikan. Ini adalah versi contoh dari sumber yang sangat bagus ini . Ini didasarkan pada GDALdan numpyyang bukan bagian dari distribusi python standar. Terima kasih kepada @Mike Toews karena menunjukkan Binari Windows Tidak Resmi untuk Paket Ekstensi Python untuk membuat instalasi dan penggunaan yang cepat dan mudah.

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
sumber
Sepertinya ini hanya versi yang kurang umum, kurang fleksibel dari apa yang ditawarkan MerseyViking di atas?
WileyB