Verwenden Sie ST_Difference, um überlappende Features zu entfernen?

11

Ich versuche, ST_Difference zu verwenden, um mit PostGis 2.1 (und Postgres SQL 9.3) einen Satz von Polygonen (process.trimmedparcelsnew) zu erstellen, die keinen der Bereiche enthalten, die von einem anderen Satz von Polygonen (test.single_geometry_1) abgedeckt werden. Hier ist meine Frage:

CREATE TABLE processing.trimmedparcelsnew AS
SELECT
    orig.id, ST_Difference(orig.geom, cont.geom) AS difference
FROM 
    test.single_geometry_1 cont,
    test.multi_geometry_1 orig;

Die resultierenden Polygone wurden jedoch nicht zugeschnitten, sondern scheinen dort geteilt worden zu sein, wo sie sich mit der anderen Schicht schneiden. Ich habe versucht, nur die Auswahl auszuführen, ohne das Ergebnis in eine Tabelle und alles andere, was mir einfällt, einzutragen, aber ich kann diese Funktion anscheinend nicht zum Laufen bringen.

Ich habe ein Bild des Ergebnisses angehängt

Geben Sie hier die Bildbeschreibung ein


Nach Kommentaren habe ich versucht, eine WHERE-Klausel hinzuzufügen. Ich möchte, dass die Pakete, die keine Schnittpunkte haben, und die Schnittbereiche der anderen Pakete entfernt werden (die Ebene test.single_geometry repräsentiert die Kontamination, die ich von meinen Paketen entfernen möchte). Ich habe eine Kreuzung versucht, aber natürlich möchte ich eigentlich die Nicht-Kreuzungen, also versuche ich jetzt eine Disjunktion. Ich habe auch versucht, das orig zu meiner Tabelle hinzuzufügen, aber die Dokumentation für ST_Difference ( http://postgis.net/docs/ST_Difference.html ) besagt, dass es genau die Geometrie zurückgibt, die ich brauche (eine Geometrie, die den Teil der Geometrie A darstellt, der schneidet sich nicht mit der Geometrie B), daher bin ich verwirrt, warum ich stattdessen das ursprüngliche Polygon in meiner Tabelle haben möchte. Wie auch immer, hier ist mein geänderter Code:

CREATE TABLE processing.trimmedparcelsnew AS
SELECT
    orig.id, ST_Difference(orig.geom, cont.geom) AS difference, orig.geom AS geom
FROM 
    test.single_geometry_1 cont,
    test.multi_geometry_1 orig
WHERE ST_Disjoint(orig.geom, cont.geom);

Nach der Antwort von dbaston habe ich jetzt versucht:

CREATE TABLE processing.parcels_trimmed AS
SELECT id, COALESCE(ST_Difference(geom, (SELECT ST_Union(b.geom) 
                                         FROM test.single_geometry_1 b
                                         WHERE ST_Intersects(a.geom, b.geom)
                                         AND a.id != b.id)), a.geom)
FROM test.multi_geometry_1 a;

Das Ergebnis ist nur eine Kopie von test.multi_geometry_1. Obwohl jetzt die Aufteilung nicht mehr stattfindet.

Ich habe die frühere Version ausprobiert, erhalte aber erneut eine Kopie von test.multi_geometry_1:

CREATE TABLE processing.parcels_trimmed_no_coalesce AS
SELECT id, COALESCE(ST_Difference(geom, (SELECT ST_Union(b.geom) 
                                         FROM test.single_geometry_1 b
                                         WHERE ST_Intersects(a.geom, b.geom)
                                         AND a.id != b.id)), a.geom)
FROM test.multi_geometry_1 a;

Ich frage mich, ob ich noch etwas falsch mache. Die folgende Erklärung lautet:

DROP TABLE IF EXISTS processing.parcels_trimmed_no_coalesce;

Und ich führe die Abfragen aus dem PostgreSQL SQL-Abfragefenster und Openjump aus.

Die Aussage, mit der ich die Tabelle sehe, lautet:

SELECT * FROM processing.parcels_trimmed_no_coalesce;

Im Interesse der Vereinfachung habe ich diese Abfrage jetzt auf nur Folgendes reduziert:

SELECT id, COALESCE(ST_Difference(geom, (SELECT ST_Union(b.geom) 
                                         FROM test.geometriestocutagainst b
                                         WHERE ST_Intersects(a.geom, b.geom)
                                         AND a.id != b.id)), a.geom)
FROM test.geometriestocut a;

Dies führt immer noch nur zu den ursprünglichen Polygonen (test.geometriestocut), wenn das gewünschte Ergebnis das Original ist, das gegen test.geometriestocut gegen getrimmt wurde.

Mart
quelle
Sie haben keine WHEREKlausel angegeben, sodass die resultierende Tabelle möglicherweise eine Polynomerweiterung enthält. Wie viele Zeilen sind in trimmedparcelsnew?
Vince
Wenn Sie nur den Unterschied dort möchten, wo sie sich schneiden, können Sie versuchen, WHERE ST_Intersects (orig.geom, cont.geom) hinzuzufügen. Andernfalls ist der Unterschied zwischen zwei Polygonen, die sich nicht schneiden, das ursprüngliche Polygon.
John Powell
Es gibt 24 Zeilen in beschnittenem Paket neu. Ich möchte den Unterschied auch dann, wenn sie sich nicht überschneiden. Würde ich also korrigieren, dass ich orig.geom in der Tabelle anstelle von Unterschied verwenden muss?
Mart
Ein disjunkter Test sollte eine Polynomerweiterung erzeugen - jedes Merkmal in cont wird einmal für jedes Merkmal in orig angezeigt , das sich nicht überlappt, und der Unterschied wird niemals eine Eingabegeometrie ändern
Vince
Ok, danke für die Klarstellung, aber ich bin mir immer noch nicht sicher, warum der ursprüngliche Code nicht funktioniert. Wenn ST_Difference (orig.geom, cont.geom) die Geometrien in a zurückgibt, die sich nicht mit b schneiden, warum enthält die Tabelle dann die geteilten Geometrien und nicht die Geometrien in a, die b nicht schneiden.
Mart

Antworten:

14

Mit einem Self-Join können Sie die Beziehung zwischen Paaren zweier Features bearbeiten. Ich glaube jedoch nicht, dass Sie an Paaren interessiert sind: Für jedes Feature möchten Sie die Beziehung zwischen diesem Feature und allen anderen Features in Ihrem Datensatz bearbeiten. Sie können dies mit einem Unterabfrageausdruck erreichen:

CREATE TABLE parcels_trimmed AS
SELECT id, ST_Difference(geom, (SELECT ST_Union(b.geom) 
                                FROM parcels b
                                WHERE ST_Intersects(a.geom, b.geom)
                                  AND a.id != b.id))
FROM parcels a;

Möglicherweise sehen Sie jedoch etwas Seltsames in den Ergebnissen. Pakete ohne Überlappungen werden vollständig verworfen! Das liegt daran, dass das ST_UnionAggregat in einem leeren Recordset sein NULLwird und ST_Difference(geom, NULL)ist NULL. Um dies zu glätten, müssen Sie Ihren ST_DifferenceAnruf in Folgendes einschließen COALESCE:

CREATE TABLE parcels_trimmed AS
SELECT id, COALESCE(ST_Difference(geom, (SELECT ST_Union(b.geom) 
                                         FROM parcels b
                                         WHERE ST_Intersects(a.geom, b.geom)
                                         AND a.id != b.id)), a.geom)
FROM parcels a;

Dies bedeutet , dass , wenn das Ergebnis ST_Differenceist NULL, der zusammengefügte Ausdruck auf die ursprüngliche Geometrie bewerten wird.

Mit der obigen Abfrage werden überlappende Bereiche vollständig aus Ihrer Domain entfernt. Wenn Sie stattdessen einen Gewinner auswählen möchten, können Sie dies a.id < b.idanstelle eines anderen Kriteriums tun a.id != b.id.

dbaston
quelle
Vielen Dank für die Antwort. Leider habe ich Probleme, dies für mich zum Laufen zu bringen. Stattdessen habe ich nur das ursprüngliche Polygon (a). Ich werde meine Frage mit weiteren Informationen bearbeiten. Danke noch einmal.
Mart
2

Ich hatte das gleiche Problem wie du. Ich weiß nicht, ob Sie bereits eine Lösung für Ihr Problem gefunden haben, aber ich habe die oben akzeptierte Antwort geändert und das bekommen, was ich wollte.

CREATE TABLE parcels_trimmed AS
SELECT id, COALESCE(ST_Difference(geom, (SELECT ST_Collect(b.geom) 
                                         FROM parcels b
                                         WHERE ST_Intersects(a.geom, b.geom)
                                         )), a.geom)
FROM parcels a;
Dag
quelle
1

Ich verwende ST_DifferenceAgg () aus den PostGIS-Addons . Sie müssen die beiden Tabellen zusammenführen, einen eindeutigen Bezeichner und einen Index für die Geometriespalte haben. Hier ist ein kurzes Beispiel:

WITH overlappingtable AS (
  SELECT 1 id, ST_GeomFromText('POLYGON((0 1, 3 2, 3 0, 0 1), (1.5 1.333, 2 1.333, 2 0.666, 1.5 0.666, 1.5 1.333))') geom
  UNION ALL
  SELECT 2 id, ST_GeomFromText('POLYGON((1 1, 3.8 2, 4 0, 1 1))')
  UNION ALL
  SELECT 3 id, ST_GeomFromText('POLYGON((2 1, 4.6 2, 5 0, 2 1))')
  UNION ALL
  SELECT 4 id, ST_GeomFromText('POLYGON((3 1, 5.4 2, 6 0, 3 1))')
  UNION ALL
  SELECT 5 id, ST_GeomFromText('POLYGON((3 1, 5.4 2, 6 0, 3 1))')
)
SELECT a.id, ST_DifferenceAgg(a.geom, b.geom) geom
FROM overlappingtable a,
     overlappingtable b
WHERE a.id = b.id OR -- Make sure to pass at least once the polygon with itself
      ((ST_Contains(a.geom, b.geom) OR -- Select all the containing, contained and overlapping polygons
        ST_Contains(b.geom, a.geom) OR
        ST_Overlaps(a.geom, b.geom)) AND
       (ST_Area(a.geom) < ST_Area(b.geom) OR -- Make sure bigger polygons are removed from smaller ones
        (ST_Area(a.geom) = ST_Area(b.geom) AND -- If areas are equal, arbitrarily remove one from the other but in a determined order so it's not done twice.
         a.id < b.id)))
GROUP BY a.id
HAVING ST_Area(ST_DifferenceAgg(a.geom, b.geom)) > 0 AND NOT ST_IsEmpty(ST_DifferenceAgg(a.geom, b.geom));

Dadurch werden die überlappenden Teile mit dem größten überlappenden Polygon zusammengeführt. Wenn Sie den überlappenden Teil getrennt halten möchten, sehen Sie sich das Beispiel ST_splitAgg () an.

Pierre Racine
quelle