Wie erstelle ich Linien, um Unterschiede zwischen Polygon-Features in PostGIS zu visualisieren?

15

Ich habe eine PostGIS-Tabelle polygon_bmit einigen Polygon-Features. Es gibt auch eine Tabelle polygon_amit den gleichen Polygonen wie polygon_bmit geringfügigen Änderungen. Jetzt möchte ich Linien erstellen, um die Unterschiede zwischen den Polygon-Features zu visualisieren.

Bildbeschreibung hier eingeben Bildbeschreibung hier eingeben Bildbeschreibung hier eingeben

Das nehme ich an ST_ExteriorRingund ST_Differencewerde den Job machen, aber die WHERE-Klausel scheint ziemlich knifflig zu sein.

CREATE VIEW line_difference AS SELECT
row_number() over() AS gid,
g.geom::geometry(LineString, yourSRID) AS geom
FROM 
    (SELECT
    (ST_Dump(COALESCE(ST_Difference(ST_ExteriorRing(polygon_a.geom), ST_ExteriorRing(polygon_b.geom))))).geom AS geom
    FROM polygon_a, polygon_b
    WHERE 
    -- ?
    ) AS g;

Kann mir jemand helfen?

EDIT 1

Wie von 'tilt' gepostet, habe ich versucht, ST_Overlaps(polygon_a.geom, polygon_b.geom) AND NOT ST_Touches(polygon_a.geom, polygon_b.geom)aber das Ergebnis ist nicht wie erwartet.

CREATE VIEW line_difference AS SELECT
row_number() over() AS gid,
g.geom::geometry(LineString, your_SRID) AS geom
FROM 
    (SELECT
    (ST_Dump(COALESCE(ST_Difference(ST_ExteriorRing(polygon_a.geom), ST_ExteriorRing(polygon_b.geom))))).geom AS geom
    FROM polygon_a, polygon_b
    WHERE 
    ST_Overlaps(polygon_a.geom, polygon_b.geom) AND NOT ST_Touches(polygon_a.geom, polygon_b.geom))
     AS g;

Bildbeschreibung hier eingeben

BEARBEITEN 2

workupload.com/file/J0WBvRBb (Beispieldatensatz)


Ich habe versucht, die Polygone vor der Verwendung von ST_Difference in Mehrfachlinien umzuwandeln, aber die Ergebnisse sind immer noch seltsam.

CREATE VIEW multiline_a AS SELECT
row_number() over() as gid,
ST_Union(ST_ExteriorRIng(polygon_a.geom))::geometry(multilinestring, 4326) AS geom
FROM
polygon_a;

CREATE VIEW multiline_b AS SELECT
row_number() over() as gid,
ST_Union(ST_ExteriorRIng(polygon_b.geom))::geometry(multilinestring, 4326) AS geom
FROM
polygon_b;

CREATE VIEW line_difference AS SELECT
row_number() over() as gid,
g.geom
FROM
    (SELECT
    (ST_Dump(COALESCE(ST_Difference(multiline_a.geom, multiline_b.geom)))).geom::geometry(linestring, 4326) AS geom
    FROM
    multiline_a, multiline_b)
As g;

Bildbeschreibung hier eingeben

Mondmeer
quelle
Sieht eher nach einer Topologiefrage aus. Sie möchten Segmente identifizieren, die nicht von der anderen Ebene abgedeckt werden. Ich habe nicht viel mit der PostGIS-Topologie gearbeitet und kann Ihnen keine direkte Antwort geben, aber ich schlage vor, Sie gehen näher darauf ein.
Thomas
Interessant, haben Sie einen Beispieldatensatz zum Download?
Huckfinn

Antworten:

10

Hier sind ein paar neue Tricks:

  • EXCEPTUm Geometrien aus jeder Tabelle zu entfernen, die identisch sind, können wir uns nur auf Geometrien konzentrieren, die für jede Tabelle ( A_onlyund B_only) eindeutig sind .
  • ST_Snap um genaue Knoten für Overlay-Operatoren zu erhalten.
  • Verwenden Sie den ST_SymDifferenceOverlay-Operator, um den symmetrischen Unterschied zwischen den beiden Geometriesätzen zu ermitteln und die Unterschiede anzuzeigen. Update : ST_Differencezeigt das gleiche Ergebnis für dieses Beispiel. Sie können beide Funktionen ausprobieren, um zu sehen, was sie erhalten.

Dies sollte bekommen, was Sie erwarten:

-- CREATE OR REPLACE VIEW polygon_SymDifference AS
SELECT row_number() OVER () rn, *
FROM (
  SELECT (ST_Dump(ST_SymDifference(ST_Snap(A, B, tol), ST_Snap(B, A, tol)))).*
  FROM (
    SELECT ST_Union(DISTINCT A_only.geom) A, ST_Union(DISTINCT B_only.geom) B, 1e-5 tol
    FROM (
      SELECT ST_Boundary(geom) geom FROM polygon_a
      EXCEPT SELECT ST_Boundary(geom) geom FROM polygon_b
    ) A_only,
    (
      SELECT ST_Boundary(geom) geom FROM polygon_b
      EXCEPT SELECT ST_Boundary(geom) geom FROM polygon_a
    ) B_only
  ) s
) s;

 rn |                                        geom
----+-------------------------------------------------------------------------------------
  1 | LINESTRING(206.234028204842 -92.0360704110685,219.846021625456 -92.5340701703592)
  2 | LINESTRING(18.556700448873 -36.4496098325257,44.44438533894 -40.5104231486146)
  3 | LINESTRING(-131.974995802602 -38.6145334122719,-114.067738329597 -39.0215165366584)
(3 rows)

drei Zeilen


Um diese Antwort ein bisschen mehr auszupacken, ST_Boundaryerhält der erste Schritt mit die Grenze jedes Polygons und nicht nur das Äußere. Wenn es zum Beispiel Löcher gäbe, würden diese durch die Grenze verfolgt.

Mit dieser EXCEPTKlausel werden Geometrien aus A, die Teil von B sind, und Zeilen aus B, die Teil von A sind, entfernt. Dadurch wird die Anzahl der Zeilen verringert, die nur Teil von A und nur Teil von B sind. ZB um A_only zu bekommen:

SELECT ST_Boundary(geom) geom FROM polygon_a
EXCEPT SELECT ST_Boundary(geom) geom FROM polygon_b

Hier sind die 6 Zeilen von A_only und 3 Zeilen von B_only: A_only Nur B_only

Als nächstes ST_Union(DISTINCT A_only.geom)wird verwendet, um die Linien in einer einzigen Geometrie zu kombinieren, normalerweise einem MultiLineString.

Mit ST_Snap werden Knoten von einer Geometrie zu einer anderen gefangen. Zum Beispiel ST_Snap(A, B, tol)wird die A-Geometrie genommen und weitere Knoten aus der B-Geometrie hinzugefügt oder in die B-Geometrie verschoben, wenn sie sich innerhalb des tolAbstands befinden. Es gibt wahrscheinlich mehrere Möglichkeiten, diese Funktionen zu verwenden, aber die Idee ist, Koordinaten von jeder Geometrie zu erhalten, die genau zueinander sind. So sehen die beiden Geometrien nach dem Einrasten aus:

Ein schnappte B schnappte

Um Unterschiede anzuzeigen, können Sie entweder ST_SymDifferenceoder auswählen ST_Difference. Sie zeigen beide das gleiche Ergebnis für dieses Beispiel.

Mike T
quelle
Gute Antwort. Ich habe mich gefragt, was Sie verwendet haben, um die Ergebnisse Ihrer Zwischenabfragen zu visualisieren. Es sah nicht sofort aus wie QGIS, und ich denke, es ist etwas, das etwas schneller macht.
RoperMaps
1
Ich benutze JTS Testbuilder zum Anzeigen und Verarbeiten von Geometrien. Es ist eine verwandte Geometrie-Engine zu GEOS und Shapely, hat aber eine Java-basierte GUI.
Mike T
Gibt es eine Möglichkeit zum Ignorieren / Überspringen von Problemen mit nicht geknoteten Schnittpunkten zwischen LINESTRING? Alle Eingabepolygone scheinen in Ordnung zu sein (mit dem QGIS-Geometrieprüfer überprüft).
Eclipsed_by_the_moon
1
'ST_Boundary (ST_SnapToGrid (geom, 0.001)' anstelle von 'ST_Boundary (geom)' löst das Problem.
Eclipsed_by_the_moon
6

Ich denke, es ist ein bisschen knifflig, weil die beiden Polygone unterschiedliche Knotensätze haben (grünes Polygon A, rote unterschiedliche Segmente von Polyon B). Der Vergleich der Segmente beider Polygone gibt einen Hinweis darauf, welche Segmente von Polygon B modifiziert werden.

Knotenpolygon A

poly a

Knoten des "unterschiedlichen" Segmentpolygons B

seg diff

Leider zeigt dies nur den Unterschied in der Segmentstruktur, aber ich hoffe, es ist ein Ausgangspunkt und es funktioniert so:

Nach einem Download- und Entpackvorgang habe ich den Datensatz mit PostgrSQL 9.46, PostGIS 2.1 unter Debian Linux Jessie mit den Befehlen importiert.

$ createdb gis-se
$ psql gis-se < /usr/share/postgis-2.1/postgis.sql
$ psql gis-se < /usr/share/postgis-2.1/spatial_ref_sys.sql
$ shp2pgsql -S polygon_a | psql gis-se
$ shp2pgsql -S polygon_b | psql gis-se

Unter der Annahme, dass die Segmente von Polygon A nicht in B und umgekehrt sind, versuche ich, den Unterschied zwischen den Segmenten beider Polygonsätze zu bilden, wobei ich die Segmentzugehörigkeit zu den Polygonen in jeder Gruppe (A oder B) vernachlässige. Aus didaktischen Gründen formuliere ich das SQL-Zeug in mehreren Ansichten.

Entsprechend diesem GIS-SE-Beitrag zerlege ich die beiden Polygone in Segmenttabellen segments_aundsegments_b

-- Segments of the polygon A
CREATE VIEW segments_a AS SELECT sp, ep
FROM
   -- extract the endpoints for every 2-point line segment for each linestring
   (SELECT
      ST_PointN(geom, generate_series(1, ST_NPoints(geom)-1)) as sp,
      ST_PointN(geom, generate_series(2, ST_NPoints(geom)  )) as ep
    FROM
    -- extract the individual linestrings
     (SELECT (ST_Dump(ST_Boundary(geom))).geom
      FROM polygon_a
     ) AS linestrings
    -- be sure that nothing is scrambled
    ORDER BY sp, ep
) AS segments;

Das Segmenttabellenpolygon A:

SELECT 
  st_astext(sp) AS sp, 
  st_astext(ep) AS ep 
FROM segments_a 
LIMIT 3;
                    sp                     |                 ep
-------------------------------------------+--------------------------------------------
POINT(-292.268907321861 95.0342877387557)  | POINT(-287.118411917425 99.4165242769195)
POINT(-287.118411917425 99.4165242769195)  | POINT(-264.62129248575 93.2470010145007)
POINT(-277.459563916327 -44.5629543976138) | POINT(-292.268907321861 95.03428773875

Das gleiche Verfahren wurde auf Polygon B angewendet.

-- Segments of the polygon B
CREATE VIEW segments_b AS SELECT sp, ep
FROM
   -- extract the endpoints for every 2-point line segment for each linestring
   (SELECT
      ST_PointN(geom, generate_series(1, ST_NPoints(geom)-1)) as sp,
      ST_PointN(geom, generate_series(2, ST_NPoints(geom)  )) as ep
    FROM
    -- extract the individual linestrings
     (SELECT (ST_Dump(ST_Boundary(geom))).geom
      FROM polygon_b
     ) AS linestrings
    -- be sure that nothing is scrambled
    ORDER BY sp, ep
) AS segments;

Das Segmenttabellenpolygon B

SELECT
  st_astext(sp) AS sp, 
  st_astext(ep) AS ep 
FROM segments_b 
LIMIT 3;
                    sp                     |                    ep
-------------------------------------------+-------------------------------------------
POINT(-292.268907321861 95.0342877387557)  | POINT(-287.118411917425 99.4165242769195)
POINT(-287.118411917425 99.4165242769195)  | POINT(-264.62129248575 93.2470010145007)
POINT(-277.459563916327 -44.5629543976138) | POINT(-292.268907321861 95.0342877387557)
...                        

Ich kann eine Differenztabellenansicht namens erstellen segments_diff_{a,b}. Die Differenz ergibt sich aus dem Nichtauftreten sortierter Start- oder Endpunkte in der Segmentmenge A und B.

CREATE VIEW segments_diff_a AS
SELECT st_makeline(b.sp, b.ep) as geom
FROM segments_b as b
LEFT JOIN segments_a as a ON (a.sp=b.sp and a.ep = b.ep)
-- filter segments without corresponding stuff in polygon A
WHERE a.sp IS NULL;

segs diff b

Und das Komplementäre:

CREATE VIEW segments_diff_b AS
SELECT st_makeline(a.sp, a.ep) as geom
FROM segments_a as a
LEFT JOIN segments_b as b ON (a.sp=b.sp and a.ep = b.ep)
-- filter segments without corresponding stuff in polygon B
WHERE b.sp IS NULL;

segs diff a

Fazit: Um ein korrektes Ergebnis für die kleinen Segmente zu erhalten, die Sie mit dem roten Pfeil markiert haben, müssen beide Polygone die gleiche Knotenstruktur aufweisen und ein Schnittschritt auf Knotenebene (Einfügen von Eckpunkten von Polygon A in B) ist erforderlich. Die Überschneidung könnte erfolgen durch:

CREATE VIEW segments_bi AS 
SELECT distinct sp, ep
FROM (
 SELECT
      ST_PointN(geom, generate_series(1, ST_NPoints(geom)-1)) as sp,
      ST_PointN(geom, generate_series(2, ST_NPoints(geom)  )) as ep
 FROM (
   SELECT st_difference(b.seg, a.seg) as geom FROM 
      segments_diff_a as a, segments_diff_b as b 
      WHERE st_intersects(a.seg, b.seg)
    ) as cut
  ) as segments
  WHERE sp IS NOT NULL AND ep IS NOT NULL
ORDER BY sp, ep;

Aber mit seltsamen Ergebnissen ...

Schnittversion

huckfinn
quelle
Danke für Ihr Bemühen. Nun, die Ergebnisse sind seltsam. Ich frage mich nur, ob ST_HausdorffDistance () helfen kann, die Frage zu beantworten: gis.stackexchange.com/questions/180593/…
Lunar Sea
Hm, st_haudorffdistance gibt Ihnen ein Ähnlichkeitsmaß, nicht die gewünschten Segmente (rote Pfeile zeigen auf).
Huckfinn
Es ist nur eine Idee, mit ST_HausdorffDistance können die Geometrien beider Tabellen verglichen werden. Waren die Polygone räumlich nicht gleich, muss ich Linien erzeugen. Ich weiß nur nicht, wie ich das machen soll.
Lunar Sea
Es scheint eine Frage der Präzision und Topologie zu sein ( gis.stackexchange.com/a/182838/26213 und webhelp.esri.com/arcgisdesktop/9.2/… )
huckfinn
1

Im Beispiel bedeutet die Änderung, dass sich geänderte Features aus der neuen Tabelle immer mit Features aus der alten Tabelle überschneiden. Damit wären Sie erledigt

ST_Overlaps (Geoma, Geomb) AND! ST_Touches (Geoma, Geomb)

Die Negation bei Berührungen besteht darin, dass sich Features auch überlappen, wenn nur ihre Ränder die gleichen Eckpunkte aufweisen.

Neigung
quelle