Comment créer des lignes pour visualiser les différences entre les entités surfaciques dans PostGIS?

J'ai une table PostGIS polygon_bavec quelques entités polygonales. Il y a aussi une table polygon_aqui contient les mêmes polygones polygon_bmais avec des changements mineurs. Maintenant, je veux créer des lignes pour visualiser les différences entre les entités surfaciques.

Je suppose que ST_ExteriorRinget ST_Differencefera le travail, mais la clause WHERE semble être assez délicate.

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;

Quelqu'un peut-il m'aider?

EDIT 1

Comme indiqué par «tilt», j'ai essayé ST_Overlaps(polygon_a.geom, polygon_b.geom) AND NOT ST_Touches(polygon_a.geom, polygon_b.geom)mais le résultat n'est pas comme prévu.

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;

EDIT 2

workupload.com/file/J0WBvRBb (exemple d'ensemble de données)

J'ai essayé de transformer les polygones en multilignes avant d'utiliser ST_Difference, mais les résultats sont toujours étranges.

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;

Voici quelques nouvelles astuces, en utilisant:

• EXCEPTpour supprimer les géométries des deux tables qui sont identiques, nous pouvons donc nous concentrer uniquement sur les géométries uniques à chaque table ( A_onlyet B_only).
• ST_Snap pour obtenir un signe de tête exact pour les opérateurs de superposition.
• Utilisez l' ST_SymDifferenceopérateur de superposition pour trouver la différence symétrique entre les deux jeux de géométrie pour afficher les différences. Mise à jour : ST_Differenceaffiche le même résultat pour cet exemple. Vous pouvez essayer l'une ou l'autre fonction pour voir ce qu'ils obtiennent.

Cela devrait obtenir ce que vous attendez:

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

Pour décompresser un peu plus cette réponse, la première étape avec ST_Boundaryobtient la limite de chaque polygone, plutôt que juste l'extérieur. Par exemple, s'il y avait des trous, ceux-ci seraient tracés par la frontière.

La EXCEPTclause est utilisée pour supprimer des géométries de A qui font partie de B et des lignes de B qui font partie de A. Cela réduit le nombre de lignes qui font partie de A uniquement et une partie de B uniquement. Par exemple, pour obtenir A_only:

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

Voici les 6 lignes de A_only et 3 lignes de B_only:

Ensuite, ST_Union(DISTINCT A_only.geom)est utilisé pour combiner le dessin au trait en une seule géométrie, généralement une chaîne MultiLineString.

ST_Snap est utilisé pour accrocher les nœuds d'une géométrie à une autre. Par exemple ST_Snap(A, B, tol), prendra la géométrie A et ajoutera plus de nœuds de la géométrie B, ou les déplacera vers la géométrie B, s'ils sont à toldistance. Il existe probablement plusieurs façons d'utiliser ces fonctions, mais l'idée est d'obtenir des coordonnées exactes les unes des autres à partir de chaque géométrie. Ainsi, les deux géométries après la capture ressemblent à ceci:

Et pour montrer les différences, vous pouvez choisir d'utiliser soit ST_SymDifferenceou ST_Difference. Ils montrent tous les deux le même résultat pour cet exemple.

Je pense que c'est un peu délicat, en raison des différents ensembles de nœuds de vos deux polygones (polygone vert A, segments rouges différents du polyon B). La comparaison des segments des deux polygones donne un indice sur les segments du polygone B qui seront modifiés.

Nœuds polygone A

Noeuds du polygone de segments "différents" B

Malheureusement, cela ne montre que la différence dans la structure des segments, mais j'espère que c'est un point de départ et que cela fonctionne comme ceci:

Après un processus de téléchargement et de décompression, j'ai importé l'ensemble de données à l'aide de PostgrSQL 9.46, PostGIS 2.1 sous Debian Linux Jessie avec les commandes.

$ 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

En supposant que les segments du polygone A ne sont pas en B et vice-versa, j'essaie de créer la différence entre les segments des deux ensembles de polygones, en négligeant l'appartenance des segments aux polygones de chaque groupe (A ou B). Pour des raisons didactiques, je formule le truc SQL dans plusieurs vues.

Correspondant à ce post GIS-SE , je décompose les deux polygones en tables de segments segments_aetsegments_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;

Le polygone de table de segments 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

La même procédure a été appliquée au polygone B.

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

Le polygone de table de segments 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)
...                        

Je peux créer une vue de table des différences nommée segments_diff_{a,b}. La différence est donnée par la non-occurrence de points de début ou de fin triés dans l'ensemble de segments A et 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;

Et les trucs complémentaires:

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;

Conclusion: Pour obtenir un résultat correct pour les petits petits segments que vous avez marqués avec la flèche rouge, les deux polygones doivent avoir la même structure de nœud et une étape d'intersection au niveau du nœud (insertion de sommets du polygone A dans B) est requise. L'intersection pourrait se faire par:

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;

Mais avec des résultats étranges ...

En regardant l'exemple, la modification implique que les entités de la nouvelle table qui ont été modifiées seront toujours des entités se chevauchant de l'ancienne table. Par conséquent, vous en auriez fini avec

ST_Overlaps (geoma, geomb) ET! ST_Touches (geoma, geomb)

La négation au toucher est due au fait que les entités se chevauchent également si seules leurs bordures partagent les mêmes emplacements de sommets.