¿Cómo crear líneas para visualizar las diferencias entre las características de polígono en PostGIS?

15

Tengo una tabla PostGIS polygon_bcon algunas características de polígono. También hay una tabla polygon_aque contiene los mismos polígonos polygon_bpero con cambios menores. Ahora quiero crear líneas para visualizar las diferencias entre las características del polígono.

ingrese la descripción de la imagen aquí ingrese la descripción de la imagen aquí ingrese la descripción de la imagen aquí

Supongo que ST_ExteriorRingy ST_Differencea hacer el trabajo, pero la cláusula WHERE parece ser bastante complicado.

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;

¿Alguien puede ayudarme?

EDITAR 1

Según lo publicado por 'tilt', lo he intentado ST_Overlaps(polygon_a.geom, polygon_b.geom) AND NOT ST_Touches(polygon_a.geom, polygon_b.geom)pero el resultado no es el esperado.

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;

ingrese la descripción de la imagen aquí

EDITAR 2

workupload.com/file/J0WBvRBb (conjunto de datos de ejemplo)


Intenté convertir los polígonos en multilíneas antes de usar ST_Difference, pero los resultados aún son extraños.

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;

ingrese la descripción de la imagen aquí

Mar lunar
fuente
Parece más una pregunta de topología. Desea identificar segmentos que no están cubiertos por la otra capa. No he trabajado mucho con la topología de PostGIS y no puedo darle una respuesta directa, pero le sugiero que investigue más sobre esto.
Thomas
Interesante, ¿tiene un conjunto de datos de ejemplo para descargar?
huckfinn

Respuestas:

10

Aquí hay algunos trucos nuevos, usando:

  • EXCEPTpara eliminar geometrías de cualquier tabla que sean iguales, de modo que podamos centrarnos solo en geometrías que sean exclusivas de cada tabla ( A_onlyy B_only).
  • ST_Snap para obtener un nudo exacto para los operadores de superposición.
  • Use el ST_SymDifferenceoperador de superposición para encontrar la diferencia simétrica entre los dos conjuntos de geometría para mostrar las diferencias. Actualización : ST_Differencemuestra el mismo resultado para este ejemplo. Puede probar cualquiera de las funciones para ver qué obtienen.

Esto debería obtener lo que espera:

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

tres lineas


Para desempaquetar esta respuesta un poco más, el primer paso con ST_Boundaryobtiene el límite de cada polígono, en lugar de solo el exterior. Por ejemplo, si hubiera agujeros, estos serían trazados por el límite.

La EXCEPTcláusula se usa para eliminar las geometrías de A que son parte de B y las filas de B que son parte de A. Esto reduce el número de filas que son parte de A solamente y parte de B solamente. Por ejemplo, para obtener A_only:

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

Aquí están las 6 filas de A_only y 3 filas de B_only: A_solo B_solo

A continuación, ST_Union(DISTINCT A_only.geom)se usa para combinar el trabajo de línea en una sola geometría, generalmente una MultiLineString.

ST_Snap se usa para ajustar nodos de una geometría a otra. Por ejemplo ST_Snap(A, B, tol), tomará la geometría A y agregará más nodos de la geometría B, o los moverá a la geometría B, si están dentro de la toldistancia. Probablemente hay varias formas de usar estas funciones, pero la idea es obtener coordenadas de cada geometría que sean exactas entre sí. Entonces las dos geometrías después del ajuste se ven así:

Un roto B estalló

Y para mostrar las diferencias, puede elegir usar ST_SymDifferenceo ST_Difference. Ambos muestran el mismo resultado para este ejemplo.

Mike T
fuente
Buena respuesta. Me preguntaba qué usabas para visualizar los resultados de tus consultas intermedias. No parecía inmediatamente qgis, y ¿tal vez es algo que se procesa un poco más rápido?
RoperMaps
1
Yo uso JTS Testbuilder para ver y procesar geometrías. Es un motor de geometría relacionado con GEOS y Shapely, pero tiene una GUI basada en Java.
Mike T
¿Hay alguna manera de ignorar / omitir los problemas de 'intersección no noda entre LINESTRING'? Todos los polígonos de entrada parecen estar bien (marcados con el verificador de geometría QGIS).
eclipsed_by_the_moon
1
'ST_Boundary (ST_SnapToGrid (geom, 0.001))' en lugar de 'ST_Boundary (geom)' resuelve el problema.
eclipsado_por_la_luna
6

Creo que es un poco complicado, debido a los diferentes conjuntos de nodos de sus dos polígonos (polígono verde A, segmentos rojos diferentes del polón B). Al comparar los segmentos de ambos polígonos se obtiene una pista sobre qué segmentos del polígono B se modificarán.

Nodos polígono A

poli a

Nodos de los segmentos "diferentes" polígono B

seg diff

Desafortunadamente, esto muestra solo la diferencia en la estructura del segmento, pero espero que sea un punto de partida y funcione así:

Después de un proceso de descarga y descompresión, importé el conjunto de datos usando PostgrSQL 9.46, PostGIS 2.1 en Debian Linux Jessie con los comandos.

$ 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

Suponiendo que los segmentos del polígono A no están en B y viceversa, trato de construir la diferencia entre los segmentos de ambos conjuntos de polígonos, descuidando la pertenencia del segmento a los polígonos en cada grupo (A o B). Por razones didácticas, formulo el material SQL en varias vistas.

En correspondencia con esta publicación GIS-SE , descompongo los dos polígonos en tablas de segmentossegments_a ysegments_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;

La tabla de segmentos del polígono 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

El mismo procedimiento se aplicó al polígono 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;

La tabla de segmentos polígono 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)
...                        

Puedo construir una vista de tabla de diferencias llamada segments_diff_{a,b}. La diferencia está dada por la no ocurrencia de puntos de inicio o fin ordenados en el conjunto de segmentos A y 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

Y las cosas complementarias:

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

Conclusión: Para obtener un resultado adecuado para los pequeños segmentos pequeños que marcó con la flecha roja, ambos polígonos deben tener la misma estructura de nodo y un paso de intersección a nivel de nodo (se requiere insertar vértices del polígono A en B). La intersección podría hacerse por:

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;

Pero con resultados extraños ...

versión cortada

huckfinn
fuente
Gracias por tus esfuerzos. Bueno, los resultados son extraños. Me pregunto si ST_HausdorffDistance () puede ayudar a responder la pregunta: gis.stackexchange.com/questions/180593/…
Lunar Sea
Hm, st_haudorffdistance te da una medida de similitud, no los segmentos deseados (flechas rojas apuntando hacia).
huckfinn
Es solo una idea, ST_HausdorffDistance se puede usar para comparar las geometrías de ambas tablas. Si los polígonos no fueran espacialmente iguales, debo crear líneas. Simplemente no sé cómo hacer esto.
Mar Lunar
Parece ser una cuestión de precisión y topología ( gis.stackexchange.com/a/182838/26213 y webhelp.esri.com/arcgisdesktop/9.2/… )
huckfinn
1

Mirando el ejemplo, el cambio implica que las características de la nueva tabla que se han modificado siempre se superpondrán con las características de la tabla anterior. Por lo tanto, terminarías con

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

La negación de los toques se debe a que las características también se superponen si solo sus bordes comparten las mismas ubicaciones de vértices.

inclinación
fuente