Trennen Sie Polygone anhand der Schnittmenge mit PostGIS

36

Ich habe eine PostGIS-Tabelle mit Polygonen, in denen sich einige überschneiden. Das versuche ich zu tun:

  • Geben Sie mir für ein bestimmtes, durch id ausgewähltes Polygon alle Polygone, die sich schneiden. Grundsätzlich gilt,select the_geom from the_table where ST_Intersects(the_geom, (select the_geom from the_table where source_id = '123'))
  • Aus diesen Polygonen muss ein neues Polygon erstellt werden, sodass die Schnittmenge zu einem neuen Polygon wird. Wenn sich also Polygon A mit Polygon B schneidet, erhalte ich drei neue Polygone: A minus AB, AB und B minus AB.

Irgendwelche Ideen?

atogle
quelle
1
Hmmm, ich denke, Sie sehen, worauf Sie hinaus wollen, aber eine einfache Grafik kann Wunder bewirken, dass ich (und andere) genau sehen, was Sie wollen.
Jason
Fügte einige in meiner Antwort hinzu.
Adam Matan

Antworten:

29

Da Sie angegeben haben, dass Sie für jedes Polygon, an dem Sie interessiert sind, eine Gruppe sich überschneidender Polygone erhalten, möchten Sie möglicherweise ein sogenanntes "Polygon-Overlay" erstellen.

Dies ist nicht genau das, was Adams Lösung tut. Schauen Sie sich dieses Bild einer ABC-Kreuzung an, um den Unterschied zu erkennen:

ABC-Kreuzung

Ich glaube, Adams Lösung wird ein "AB" -Polygon erzeugen, das sowohl den Bereich von "AB! C" und "ABC" abdeckt, als auch ein "AC" -Polygon, das "AC! B" und "ABC" abdeckt, und ein " BC "Polygon, das" BC! A "und" ABC "ist. Die Ausgabepolygone "AB", "AC" und "BC" überlappen also alle den Bereich "ABC".

Eine Polygon-Überlagerung erzeugt nicht überlappende Polygone, also wäre AB! C ein Polygon und ABC ein Polygon.

Das Erstellen eines Polygon-Overlays in PostGIS ist eigentlich ganz einfach.

Grundsätzlich gibt es drei Schritte.

Schritt 1 ist das Extrahieren der Linien [Beachten Sie, dass ich den äußeren Ring des Polygons verwende, es wird etwas komplizierter, wenn Sie Löcher richtig handhaben möchten]:

SELECT ST_ExteriorRing(polygon_col) AS the_geom FROM my_table) AS lines

Schritt 2 besteht darin, die Linien zu "knoten" (an jeder Kreuzung einen Knoten zu erzeugen). Einige Bibliotheken wie JTS haben "Noder" -Klassen, mit denen Sie dies tun können, aber in PostGIS erledigt die ST_Union- Funktion dies für Sie:

SELECT ST_Union(the_geom) AS the_geom FROM (...your lines...) AS noded_lines

Schritt 3 besteht darin, alle möglichen nicht überlappenden Polygone zu erstellen, die aus all diesen Zeilen stammen können. Dies geschieht mit der Funktion ST_Polygonize :

SELECT ST_Polygonize(the_geom) AS the_geom FROM (...your noded lines...)

Sie können die Ausgabe jedes dieser Schritte in einer temporären Tabelle speichern oder sie alle in einer einzigen Anweisung kombinieren:

CREATE TABLE my_poly_overlay AS
SELECT geom FROM ST_Dump((
    SELECT ST_Polygonize(the_geom) AS the_geom FROM (
        SELECT ST_Union(the_geom) AS the_geom FROM (
            SELECT ST_ExteriorRing(polygon_col) AS the_geom FROM my_table) AS lines
        ) AS noded_lines
    )
)

Ich verwende ST_Dump, weil die Ausgabe von ST_Polygonize eine Geometriesammlung ist und es (normalerweise) praktischer ist, eine Tabelle zu haben, in der jede Zeile eines der Polygone ist, aus denen die Polygonüberlagerung besteht.

Jeff
quelle
Beachten Sie, dass ST_ExteriorRingkeine Löcher fallen. ST_Boundarybewahrt die inneren Ringe, erzeugt aber auch ein Polygon in ihnen.
jpmc26
Die Vereinigung der Außenringe stürzt ab, wenn zu viele Polygone vorhanden sind. Leider funktioniert diese "unkomplizierte" Lösung nur bei kleinen Netzabdeckungen.
Pierre Racine
13

Wenn ich richtig verstehe, möchten Sie nehmen (A ist die linke Geometrie, B ist die rechte):

Bild von A∪B http://img838.imageshack.us/img838/3996/intersectab1.png

Und extrahieren:

A ∖ AB

Bild von A ∖ AB http://img830.imageshack.us/img830/273/intersectab2.png

AB

Bild von AB http://img828.imageshack.us/img828/7413/intersectab3.png

und B ∖ AB

Bild von B ∖ AB http://img839.imageshack.us/img839/5458/intersectab4.png

Das heißt - drei verschiedene Geometrien für jedes sich überschneidende Paar.

Zunächst erstellen wir eine Ansicht aller sich überschneidenden Geometrien. Angenommen, Ihr Tabellenname lautet polygons_table:

CREATE OR REPLACE VIEW p_intersections AS    -- Create a view with the 
SELECT t1.the_geom as t1_geom,               -- intersecting geoms. Each pair
       t2.the_geom as t2_geom                -- appears once (t2.id<t2.id)
    FROM polygons_table t1, polygons_table t2  
         WHERE t1.id<t2.id AND t1.the_geom && t2.the_geom 
                           AND intersects t1.the_geom, t2.the_geom;

Jetzt haben wir eine Ansicht (praktisch eine Nur-Lese-Tabelle), in der Paare sich überschneidender Geoms gespeichert sind, wobei jedes Paar aufgrund der t1.id<t2.idBedingung nur einmal angezeigt wird .

Jetzt sammeln wir Ihre Kreuzungen - A∖AB, ABund B∖ABunter Verwendung von SQL ist UNIONauf allen drei Abfragen:

--AB:
SELECT ST_intersection(t1.the_geom, t2.the_geom) 
    AS geom 
    FROM p_intersections

UNION 

--AAB:
SELECT ST_Difference(t1.the_geom, t2.the_geom) 
    AS geom 
    FROM p_intersections

UNION

--BAB:
SELECT ST_Difference(t2.the_geom, t1.the_geom) 
    AS geom 
    FROM p_intersections;

Anmerkungen:

  1. Der &&Bediener wird als Filter vor dem intersectsBediener verwendet, um die Leistung zu verbessern.
  2. Ich habe mich dafür entschieden, eine VIEWstatt einer einzigen riesigen Abfrage zu erstellen . Dies dient nur der Vereinfachung.
  3. Wenn Sie gemeint haben, ABist dies die Vereinigung und nicht die Kreuzung von Aund B- Verwenden Sie ST_Union anstelle von st_intersection bei der UNIONAbfrage an den entsprechenden Stellen.
  4. Das Zeichen ist ein Unicode-Zeichen für Set difference. Entfernen Sie es aus dem Code, wenn es Ihre Datenbank verwirrt.
  5. Bilder mit freundlicher Genehmigung von Wikimedia Commons in der Kategorie nette Mengenlehre .
Adam Matan
quelle
Mein Bug Ticket auf meta: meta.gis.stackexchange.com/questions/79/…
Adam Matan
Schöne Erklärung! Die Ergebnisse sind die gleichen wie in der SCW-Lösung, aber sein Weg sollte schneller sein (berechnet / oder speichert / speichert keine zusätzlichen Schnittpunkte von A und B)
Stachu
Vielen Dank! Ich glaube, ich speichere keine zusätzlichen Informationen, da ich nur SQL VIEWs und keine Tabellen erstelle.
Adam Matan
Ja, das ist wahr, aber Sie berechnen zusätzliche Schnittmenge von A und B, die nicht erforderlich ist
Stachu
5
Die Bilder in dieser Antwort funktionieren nicht mehr.
Fezter
8

Was Sie beschreiben, ist die Art und Weise, wie der Union-Operator in ArcGIS arbeitet, jedoch ein wenig anders als Union oder Intersection in der GEOS-Welt. Das Shapely-Handbuch enthält Beispiele für die Funktionsweise von Sets in GEOS . Das PostGIS-Wiki bietet jedoch ein gutes Beispiel für die Verwendung von Linien, die den Trick für Sie tun sollten.

Alternativ können Sie drei Dinge berechnen:

  1. ST_Interection (geom_a, geom_b)
  2. ST_Differenz (geom_a, geom_b)
  3. ST_Differenz (geom_b, geom_a)

Dies sollten die drei Polygone sein, die Sie in Ihrem zweiten Aufzählungspunkt erwähnt haben.

scw
quelle
2
Das PostGIS Wiki Beispiel ist gut
fmark 23.07.10
Würde ST_Intersects nicht einen Booleschen Wert zurückgeben, wenn sie sich schneiden oder nicht? Ich denke, ST_Intersection würde funktionieren.
Jason
Ja, Tippfehler meinerseits - jetzt im Original behoben, danke Jason!
scw
-2

So etwas wie:

INSERT INTO new_table VALUES ((wähle id, the_geom aus old_table wo st_intersects (the_geom, (wähle the_geom aus old_table wo id = '123')) = true

BEARBEITEN: Sie benötigen die tatsächliche Schnittmenge des Polygons.

INSERT INTO new_table-Werte ((select id, ST_Intersection (the_geom, (select the_geom from old where id = 123))

Mal sehen, ob das klappt.

George Silva
quelle