Rasterwerte an Punkten mit Open Source GIS extrahieren?

21

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.

Vassilis
quelle
1
Sie haben also Punkte und müssen die Werte aus dem Raster unter diesen Punkten extrahieren oder die Rasterzellen in Punkte konvertieren. Ich überprüfe nur, bevor ich versuche, die Antwort zu finden.
Nathan W
Das erste, ich habe die Punkte und ich muss die Werte aus dem Raster unter diesen Punkten extrahieren. Danke !!
Vassilis

Antworten:

37

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/

Underdunkel
quelle
5
Die Leute finden den oben genannten Beitrag immer noch in diesem Q & A. Das Plugin ist jedoch nicht die einzige und nicht immer einfachste Lösung. Eine alternative Lösung ist die Saga-Funktion 'Rasterwerte zum Punkt hinzufügen' in der Verarbeitungs-Toolbox. Einzelheiten finden Sie in diesem Beitrag .
Ecodiv
Vorsicht. Ich habe gerade das Point Sampling Tool ausgeführt. 60.000 Punkte und 13 Raster. Die Ergebnisse haben meine 30 Stichproben für jedes Jahr nicht bestanden. Dieses Tool hat Probleme mit großen Datenmengen. Ich würde es nicht benutzen.
Wenn Sie nicht wissen, nur GIS
Trotz der erwähnten Probleme mit großen Datenmengen ist es sehr nützlich, um alle Multiband-Werte auf einmal zu extrahieren . Alle anderen QGIS-bezogenen Lösungen unterstützen nur die Extraktion eines Bands (wie GRASS r.what) oder verbieten die Verwendung von Multiband-Rastern (wie Saga - Rasterwerte zu Punkten).
EikeMike
7

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).

Pierre Racine
quelle
6

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.

Fernando MM
quelle
6

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))
  1. 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.

Ian
quelle
1
Dies ist die beste Antwort auf diese Frage
BC B.
4

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

jbaums
quelle
3

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

markusN
quelle
2
Für Grass fand ich Folgendes
Vassilis
Ein weiteres GRASS-Modul, mit dem ein Raster abgetastet werden kann, ist v.sample, das auch über die QGIS-Verarbeitungs-Toolbox verfügbar ist.
user55937
2

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')
spyrostheodoridis
quelle
1

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“ . )

Markieren Sie Irland
quelle
Nein, ich habe oder benutze FME nicht, ist es eine eigenständige Anwendung oder ein Plugin?
Vassilis
0

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

quelle
interessanter Ansatz, weil gestern anfangen zu lernen, wie man Postgis benutzt.
Vassilis
Danke, es ist ziemlich simpel, aber es funktioniert. Folgendes