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

15

Tengo una tabla PostGIS polygon_b con algunas características de polígono. También hay una tabla polygon_a que contiene los mismos polígonos que polygon_b pero con cambios menores. Ahora quiero crear líneas para visualizar las diferencias entre las entidades poligonales.

SupongoqueST_ExteriorRingyST_DifferenceharáneltrabajoperolacláusulaWHEREpareceserbastantedifícil.

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;

¿Puede alguien ayudarme?

EDIT 1

Según lo publicado por 'tilt', probé 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;

EDIT2

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

He intentado convertir los polígonos en líneas múltiples antes de usar ST_Difference, pero los resultados siguen siendo 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;

    
pregunta Lunar Sea 14.02.2016 - 13:30

3 respuestas

9

Aquí hay algunos trucos nuevos, usando:

  • EXCEPT para eliminar las geometrías de cualquiera de las tablas que sean iguales, por lo que solo podemos enfocarnos en las geometrías que son únicas para cada tabla ( A_only y B_only ).
  • ST_Snap para obtener el nodo exacto para los operadores de superposición.
  • Use el operador de superposición ST_SymDifference para encontrar la diferencia simétrica entre los dos conjuntos de geometría para mostrar las diferencias. Actualizar : ST_Difference muestra el mismo resultado para este ejemplo. Puede probar cualquiera de las funciones para ver qué obtienen.

Esto debería obtener lo que esperas:

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

Paradesempacarunpocomásestarespuesta,elprimerpasoconST_Boundaryobtieneellímitedecadapolígono,enlugardesoloelexterior.Porejemplo,sihubieraagujeros,estosseríanrastreadosporellímite.

La cláusula EXCEPT se utiliza para eliminar geometrías de A que son parte de B, y las filas de B que forman parte de A. Esto reduce el número de filas que son solo parte de A y solo B. 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:

Acontinuación,seutilizaST_Union(DISTINCTA_only.geom)paracombinarlalíneaenunaúnicageometría,normalmenteunaMultiLineString.

ST_Snapseusaparaajustarnodosdeunageometríaaotra.Porejemplo,ST_Snap(A,B,tol)tomarálageometríaAyagregarámásnodosdelageometríaB,olosmoveráalageometríaB,siestándentrodeladistanciadetol.Probablementehayvariasformasdeusarestasfunciones,perolaideaesobtenercoordenadasdecadageometríaqueseanexactasentresí.Asíquelasdosgeometríasdespuésdeajustarsevenasí:

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

    
respondido por el Mike T 02.03.2016 - 21:45
6

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

Polígono de nodos A

Nodosdelos"diferentes" segmentos polígono B

Desafortunadamente,estosolomuestraladiferenciaenlaestructuradelsegmento,peroesperoqueseaunpuntodepartidayfuncioneasí:

Despuésdeunprocesodedescargaydescompresión,heimportadoelconjuntodedatosutilizandoPostgrSQL9.46,PostGIS2.1bajoDebianLinuxJessieconloscomandos.

$ 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 las cosas de SQL en varias vistas.

Correspondiente a esta publicación de GIS-SE , descompongo ambos polígonos en tablas de segmentos segments_a y segments_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;

El polígono A de la tabla de segmentos:

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;

El polígono B de la tabla de segmentos

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 crear una vista de tabla de diferencias llamada segments_diff_{a,b} . La diferencia viene dada por la no aparición de puntos de inicio o finalización 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;

Ylascosascomplementarias:

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;

Conclusión:paraobtenerunresultadoadecuadoparalospequeñossegmentospequeñosquemarcóconlaflecharoja,ambospolígonosdebentenerlamismaestructuradenodoyserequiereunpasodeintersecciónenunniveldenodo(insertarvérticesdelpolígonoAenB).Laintersecciónsepodríahacerpor:

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

    
respondido por el huckfinn 29.02.2016 - 03:16
1

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

  

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

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

    
respondido por el tilt 14.02.2016 - 13:57

Lea otras preguntas en las etiquetas