Rasterwerte an Punkten mit Open Source GIS extrahieren?

Wie kann ich Werte aus einem Raster nach Punkten extrahieren?

Ich bevorzuge nicht in Arcgis.

Ich bevorzuge in Qgis oder Mapwindow oder anderen Open Source Gis.

QGIS "Point Sampling Tool" sollte das gesuchte Plugin sein.

Hier finden Sie eine detaillierte Beschreibung der Verwendung: http://pvanb.wordpress.com/2010/02/15/sampling-raster-values-at-point-locations-in-qgis/

Update basierend auf Paolos Kommentar:

Das Plugin ist nicht die einzige und nicht immer die einfachste Lösung. Eine alternative Lösung ist die Saga-Funktion 'Rasterwerte zum Punkt hinzufügen' in der Verarbeitungs-Toolbox. Weitere Informationen finden Sie unter http://pvanb.wordpress.com/2014/07/01/sampling-raster-values-at-point-locations-in-qgis-an-update/

In PostGIS 2.0 können Sie:

SELECT ST_Value(rast, geom) val
FROM yourrastertabe, yourpointtable
WHERE ST_Intersects(rast, geom)

Stellen Sie sicher, dass Ihr Raster beim Laden sehr klein gekachelt ist (-t 10x10 mit dem Lader).

Versuchen Sie es mit QGIS 3.2.2 und SAGA (standardmäßig in QGIS installiert): Die Funktion "Rasterwerte in Punkte" erledigt alles für Sie: Sie nimmt eine Bilddatei und konvertiert sie in eine Punktvektorform, wobei die Informationen aus dem Rasterbild übernommen werden.

Ich hatte Probleme mit den in diesem Thread erwähnten QGIS- und SAGA-GUI-Tools ( Raster values to pointsaus irgendeinem Grund ist dies fehlgeschlagen und es wurden nicht hilfreiche Fehler ausgegeben, und GRASS v.samplehat eine völlig neue Ebene erstellt, die nicht hilfreich war). Nachdem ich einige Zeit mit den GUI-Tools gescheitert war, versuchte ich dies im Feldrechner. Es hat ganz gut funktioniert und ich konnte den Prozess ein wenig besser steuern, als es die GUIs zulassen, und einige andere Berechnungen auf dem Weg machen.

Angenommen, Sie haben einen Layer mit dem Namen ptsund einen anderen mit dem Namen rast, beide im selben Koordinatensystem. Sie möchten rastbei jedem in dargestellten X, Y-Paar eine Stichprobe machen pts.

Wenn Sie den Feldrechner noch nicht benutzt haben, ist es ziemlich einfach. Sie geben Ihre Berechnung in das Feld "Ausdruck" ein, und Q gibt Ihnen in der mittleren Spalte eine Reihe von Variablen und Operationen mit Hilfetext in der rechten Spalte an. Ich werde diesen Prozess in vier Schritte unterteilen:

1. Öffnen Sie die Attributtabelle des ptsLayers, mit dem Sie ein Sample erstellen möchten.

2. Wählen Sie im Dialogfeld "Feldberechnung" aus, ob Sie ein neues Feld erstellen oder ein vorhandenes Feld in Ihrem ptsLayer ändern möchten .

3. Erstellen Sie als Nächstes einen Ausdruck, um die neue oder vorhandene ptsAttributspalte zu füllen . Sie könnten damit beginnen, den für mich funktionierenden Ausdruckscode zu ändern:

raster_value('rast', 1, make_point($x, $y))

4. Sie müssen raster_value()einen Rasterebenennamen 'rast', eine Bandnummer 1und die Punktgeometrie unter angeben make_point(). $xund $ysind Geometrievariablen, die von der Position des Punkts in jeder Zeile der Attributtabelle abhängen.

Diese Methode ermöglicht auch arithmetische Operationen wie das Subtrahieren des Werts einer anderen aufgerufenen other_rastRasterebene rast, was mir gegenüber den GUI-Tools einige Zeit erspart hat. Beispiel unten:

raster_value('rast', 1, make_point($x, $y)) - raster_value('other_rast', 1, make_point($x, $y))

Man beachte erneut , dass die drei Lagen pts, rastund other_rastmüssen für diese Methode funktioniert in demselben Koordinatensystem sein.

Die GME-Tools von Hawthorne Beyer erledigen dies problemlos über die Befehlszeile und ermöglichen ein einfaches Batching mit For-Schleifen.

isectpntrst(in="path/to/shapefile", raster="path/to/raster", field="fieldname")

GME isectpntrst Befehlsreferenz

In GRASS GIS können Sie entweder die Karte in der GUI abfragen oder http://grass.osgeo.org/gdp/html_grass64/r.what.html verwenden

http://gis-techniques.blogspot.com/2012/10/extract-raster-values-from-points.html enthält eine schrittweise Anleitung zur Verwendung des R-Raster-Pakets zum Extrahieren von Rasterwerten aus Punkten.

Sie können dies verwenden: http://www.saga-gis.org/saga_module_doc/2.1.3/shapes_grid_3.html

Es ist in der SAGA Toolbox von Qgis! Es macht alles in einem Schritt :)

Hier ist eine Funktion, die ich mit Python und GDAL geschrieben habe. Die Funktion erstellt eine Liste von Rastern und einen Pandas-Datenrahmen mit den Punktkoordinaten und gibt einen Pandas-Datenrahmen mit den Punktkoordinaten, den Schwerpunkten für die jeweiligen Rasterzellen und den jeweiligen Zellenwerten zurück. Die Funktion ist Teil des in Entwicklung befindlichen Pakets Chorospy ( hier zu finden ).

import pandas
import numpy
from osgeo import gdal

def getValuesAtPoint(indir, rasterfileList, pos, lon, lat):
    #gt(2) and gt(4) coefficients are zero, and the gt(1) is pixel width, and gt(5) is pixel height.
    #The (gt(0),gt(3)) position is the top left corner of the top left pixel of the raster.
    for i, rs in enumerate(rasterfileList):

        presValues = []
        gdata = gdal.Open('{}/{}.tif'.format(indir,rs))
        gt = gdata.GetGeoTransform()
        band = gdata.GetRasterBand(1)
        nodata = band.GetNoDataValue()

        x0, y0 , w , h = gt[0], gt[3], gt[1], gt[5]

        data = band.ReadAsArray().astype(numpy.float)
        #free memory
        del gdata

        if i == 0:
            #iterate through the points
            for p in pos.iterrows():
                x = int((p[1][lon] - x0)/w)
                Xc = x0 + x*w + w/2 #the cell center x
                y = int((p[1][lat] - y0)/h)
                Yc = y0 + y*h + h/2 #the cell center y
                try:
                    if data[y,x] != nodata:
                        presVAL = [p[1][lon],p[1][lat], '{:.6f}'.format(Xc), '{:.6f}'.format(Yc), data[y,x]]
                        presValues.append(presVAL)
                except:
                    pass
            df = pandas.DataFrame(presValues, columns=['x', 'y', 'Xc', 'Yc', rs])
        else:
            #iterate through the points
            for p in pos.iterrows():
                x = int((p[1][lon] - x0)/w)
                y = int((p[1][lat] - y0)/h)
                try:
                    if data[y,x] != nodata:
                        presValues.append(data[y,x])
                except:
                    pass
            df[rs] = pandas.Series(presValues)
    del data, band
    return df

Beispiel für die Ausführung unter der Voraussetzung, dass sich die Raster in Ihrem aktuellen Arbeitsverzeichnis befinden:

rasDf = getValuesAtPoint('.', ['raster1', 'raster2'], inPoints, 'x', 'y')

Wenn Sie Zugriff auf FME haben , können Sie einen von zwei Transformatoren in FME Workbench verwenden.

Der RasterCellCoercer (" Zerlegt alle numerischen Eingabe-Raster-Features in einzelne Punkte oder Polygone. Für jede Zelle im Raster wird ein Vektor-Feature ausgegeben.")

Die PointOnRasterValueExtractor ( „führt in Punktmerkmalen und ein einzigen Referenzraster. Der Ausgang besteht aus dem Band und Palette Wert (e) an der Stelle jeden Punktes“ . )

Schneller Gedanke:

1. gdal_polygonize.py - polygonisiert Ihr Raster-Feature
2. Fügen Sie Ihre Punktobjekte und Polygone in die PostGIS-Datenbank ein
3. Verwenden Sie die Funktion st_intersects, um alle Höhenwerte dort abzurufen, wo sich Features schneiden