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

15

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.

entrez la description de l'image ici entrez la description de l'image ici entrez la description de l'image ici

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;

entrez la description de l'image ici

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;

entrez la description de l'image ici

Mer lunaire
la source
Ressemble plus à une question de topologie. Vous souhaitez identifier les segments qui ne sont pas couverts par l'autre couche. Je n'ai pas beaucoup travaillé avec la topologie PostGIS et je ne peux pas vous donner de réponse directe, mais je vous suggère de vous pencher davantage sur cela.
Thomas
Intéressant, avez-vous un exemple de jeu de données à télécharger?
huckfinn

Réponses:

10

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)

trois lignes


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: A_only 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:

Un cassé B cassé

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.

Mike T
la source
Bonne réponse. Je me demandais ce que vous utilisiez pour visualiser les résultats de vos requêtes intermédiaires. Cela ne ressemblait pas immédiatement à qgis, et peut-être que c'est quelque chose qui rend un peu plus rapide?
RoperMaps
1
J'utilise JTS Testbuilder pour visualiser et traiter les géométries. Il s'agit d'un moteur de géométrie associé à GEOS et Shapely, mais il possède une interface graphique basée sur Java.
Mike T
Existe-t-il un moyen d'ignorer / ignorer les problèmes d '«intersection sans signe de tête entre LINESTRING»? Tous les polygones en entrée semblent bien (vérifiés avec le vérificateur de géométrie QGIS).
eclipsed_by_the_moon
1
'ST_Boundary (ST_SnapToGrid (geom, 0.001))' au lieu de 'ST_Boundary (geom)' résout le problème.
eclipsed_by_the_moon
6

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

poly a

Noeuds du polygone de segments "différents" B

seg diff

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;

segs diff b

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;

segs diff a

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

version coupée

huckfinn
la source
Merci pour vos efforts. Eh bien, les résultats sont étranges. Je me demande simplement si ST_HausdorffDistance () peut aider à répondre à la question: gis.stackexchange.com/questions/180593/…
Lunar Sea
Hm, st_haudorffdistance vous donne une mesure de similitude et non les segments souhaités (flèches rouges pointant vers).
huckfinn
C'est juste une idée, ST_HausdorffDistance peut être utilisé pour comparer les géométries des deux tables. Si les polygones ne sont pas spatialement égaux, je dois créer des lignes. Je ne sais juste pas comment faire ça.
Lunar Sea
Cela semble être une question de précision et de topologie ( gis.stackexchange.com/a/182838/26213 et webhelp.esri.com/arcgisdesktop/9.2/… )
huckfinn
1

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.

inclinaison
la source