Ich verwende PostGIS2.0 , um einige Raster- / Polygonschnitte zu erstellen . Ich habe Schwierigkeiten zu verstehen, welche Operation ich verwenden soll und wie dies am schnellsten durchgeführt werden kann. Mein Problem ist wie folgt:
- Ich habe ein Polygon und ein Raster
- Ich möchte alle Pixel finden, die in das Polygon fallen, und die Summe der Pixelwerte erhalten
- Und (aktualisiertes Problem): Ich erhalte massive Werte für einige Pixel, die im ursprünglichen Raster nicht vorhanden sind, wenn ich die Abfrage durchführe
Ich habe Schwierigkeiten zu verstehen, ob ich ST_Intersects()
oder verwenden soll ST_Intersection()
. Ich weiß auch nicht, wie ich meine Pixel am besten summieren kann. Hier ist der erste Ansatz, den ich versucht habe (# 1):
SELECT
r.rast
FROM
raster as r,
polygon as p
WHERE
ST_Intersects(r.rast, p.geom)
Dies gibt eine Liste von rast
Werten zurück, mit denen ich nicht sicher bin, was ich tun soll. Ich habe versucht, die Auswertungsstatistik mit zu berechnen ST_SummaryStats()
, bin mir aber nicht sicher, ob dies die gewichtete Summe aller Pixel ist, die innerhalb des Polygons liegen.
SELECT
(result).count,
(result).sum
FROM (
SELECT
ST_SummaryStats(r.rast) As result
FROM
raster As r,
polygon As p
WHERE
ST_Intersects(r.rast, p.geom)
) As tmp
Der andere Ansatz, den ich ausprobiert habe (# 2), verwendet ST_Intersection()
:
SELECT
(gv).geom,
(gv).val
FROM
(
SELECT
ST_Intersection(r.rast, p.geom) AS gv
FROM
raster as r,
polygon as p
WHERE
ST_Intersects(r.rast, p.geom)
) as foo;
Dies gibt eine Liste von Geometrien zurück, die ich weiter analysiere, aber ich gehe davon aus, dass dies weniger effizient ist.
Ich bin mir nicht sicher, in welcher Reihenfolge das auch am schnellsten geht. Soll ich das Polygon immer auswählen raster, polygon
oder polygon, raster
in ein Raster umwandeln, damit es ist raster, raster
?
BEARBEITEN: Ich habe Ansatz 2 mit einigen Details aus R.K.
dem Link aktualisiert .
Unter Verwendung von Ansatz 2 habe ich den folgenden Fehler in den Ergebnissen festgestellt, der ein Teil des Grundes ist, warum ich die Ausgabe nicht verstanden habe. Hier ist das Bild meines ursprünglichen Rasters und ein Umriss des Polygons, das zum Schneiden verwendet wird, überlagert:
Und hier ist das Ergebnis der Schnittmenge mit PostGIS:
Das Problem mit dem Ergebnis ist, dass Werte von 21474836 zurückgegeben werden, die nicht im ursprünglichen Raster enthalten sind. Ich habe keine Ahnung, warum dies geschieht. Ich vermute, dass es irgendwo mit kleinen Zahlen zusammenhängt (dividiert durch fast 0), aber es gibt das falsche Ergebnis zurück.
ST_SummaryStats()
für # 1 verwendet, bin mir aber nicht sicher, wie ich es für # 2 machen soll.Antworten:
Ich habe ein Tutorial zum Schneiden von Vektorpolygonen mit einer großen Rasterabdeckung mit PostGIS WKT Raster gefunden . Könnte die Antworten haben, die Sie suchen.
Das Lernprogramm verwendete zwei Datensätze, eine Punktformdatei, die zur Erzeugung von Polygonen gepuffert wurde, und eine Reihe von 13 SRTM-Höhenrastern. Zwischendurch gab es viele Schritte, aber die Abfrage, mit der Raster und Vektor geschnitten wurden, sah folgendermaßen aus:
Die Werte wurden dann folgendermaßen zusammengefasst:
Ich kenne nicht wirklich genug PostGIS, um dies zu erklären, aber es klingt nach dem, was Sie erreichen wollten. Das Tutorial soll Aufschluss über die Zwischenschritte geben. Viel Glück :)
quelle
In Bezug auf Punkt 2 in der ursprünglichen Frage verwendeten einige der Postgis 2.0-Entwicklungsversionen eine Version der GDAL-Bibliothek, die Floats in Ints umwandelt. Wenn Ihr Raster Float-Werte enthält und Sie eine Version von GDAL niedriger als 1.9.0 oder eine PostGIS 2.0-Vorabversion verwendet haben, die GDALFPolygonize () nicht ordnungsgemäß aufgerufen hat, ist möglicherweise dieser Fehler aufgetreten. Tickets in den PostGIS und GDAL Bug Trackern wurden abgelegt und geschlossen. Dieser Fehler war zum Zeitpunkt der ursprünglichen Frage aktiv.
In Bezug auf die Leistung werden Sie feststellen, dass die Verwendung
ST_Intersects(raster, geom)
viel schneller ist als die VerwendungST_Intersects(geom, raster)
. In der ersten Version wird die Geometrie gerastert und eine Schnittmenge zwischen Raster und Raum erstellt. Die zweite Version vektorisiert die Geometrie und führt eine Vektorraum-Überschneidung durch, was ein weitaus teurerer Prozess sein kann.quelle
Ich hatte auch seltsame Probleme
ST_SummaryStats
mitST_Clip
. Eine andere Abfrage der DatenST_SummaryStats
ergab, dass der minimale Wert meines Rasters 32 und der maximale Wert 300 betrug. Für die Pixelwerte in meinem Zielpolygon wurde jedoch -32700 zurückgegeben.Am Ende habe ich das Problem so gehackt:
quelle