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

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.

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;

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;

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)

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:

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:

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

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

Knoten des "unterschiedlichen" Segmentpolygons B

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;

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;

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

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.