Como criar linhas para visualizar as diferenças entre os recursos de polígono no PostGIS?

15

Eu tenho uma tabela PostGIS polygon_bcom alguns recursos de polígono. Há também uma tabela polygon_aque contém os mesmos polígonos, polygon_bmas com pequenas alterações. Agora, quero criar linhas para visualizar as diferenças entre os recursos de polígono.

insira a descrição da imagem aqui insira a descrição da imagem aqui insira a descrição da imagem aqui

Suponho que ST_ExteriorRinge ST_Differencefaça o trabalho, mas a cláusula WHERE parece ser bastante complicada.

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;

Alguém pode me ajudar?

EDIT 1

Conforme publicado por 'tilt', tentei, ST_Overlaps(polygon_a.geom, polygon_b.geom) AND NOT ST_Touches(polygon_a.geom, polygon_b.geom)mas o resultado não é o 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;

insira a descrição da imagem aqui

EDIT 2

workupload.com/file/J0WBvRBb (exemplo de conjunto de dados)


Tentei transformar os polígonos em multilinhas antes de usar ST_Difference, mas os resultados ainda são estranhos.

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;

insira a descrição da imagem aqui

Mar Lunar
fonte
Parece mais uma pergunta de topologia. Você deseja identificar segmentos que não são cobertos pela outra camada. Eu não trabalhei muito com a topologia do PostGIS e não posso dar uma resposta direta, mas sugiro que você analise mais isso.
Thomas
Interessante, você tem um exemplo de conjunto de dados para download?
huckfinn

Respostas:

10

Aqui estão alguns novos truques, usando:

  • EXCEPTpara remover geometrias das tabelas iguais, para que possamos focar apenas nas geometrias exclusivas de cada tabela ( A_onlye B_only).
  • ST_Snap para obter acenos exatos para operadores de sobreposição.
  • Use o ST_SymDifferenceoperador de sobreposição para encontrar a diferença simétrica entre os dois conjuntos de geometria para mostrar as diferenças. Atualização :ST_Difference mostra o mesmo resultado para este exemplo. Você pode tentar qualquer uma das funções para ver o que elas recebem.

Isso deve obter o que você 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)

três linhas


Para descompactar um pouco mais esta resposta, o primeiro passo com ST_Boundary obtém o limite de cada polígono, e não apenas o exterior. Por exemplo, se houvesse buracos, eles seriam rastreados pelo limite.

A EXCEPTcláusula é usada para remover geometrias de A que fazem parte de B e linhas de B que fazem parte de A. Isso reduz o número de linhas que fazem parte apenas de A e apenas parte de B. Por exemplo, para obter A_only:

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

Aqui estão as 6 linhas de A_only e 3 linhas de B_only: A_only B_only

Próximo, ST_Union(DISTINCT A_only.geom) é usado para combinar o trabalho de linha em uma única geometria, geralmente uma MultiLineString.

ST_Snap é usado para ajustar nós de uma geometria para outra. Por exemplo ST_Snap(A, B, tol), pegue a geometria A e adicione mais nós da geometria B ou mova-os para a geometria B, se estiverem dentro da toldistância. Provavelmente, existem várias maneiras de usar essas funções, mas a idéia é obter coordenadas de cada geometria que sejam exatas entre si. Portanto, as duas geometrias após o snap são assim:

Um estalou B estalou

E para mostrar diferenças, você pode optar por usar ST_SymDifferenceou ST_Difference. Ambos mostram o mesmo resultado para este exemplo.

Mike T
fonte
Boa resposta. Gostaria de saber o que você usou para visualizar os resultados de suas consultas intermediárias. Não parecia imediatamente com qgis, e eu talvez seja algo que renderize um pouco mais rápido?
RoperMaps
1
Eu uso o JTS Testbuilder para visualizar e processar geometrias. É um mecanismo de geometria relacionado ao GEOS e Shapely, mas tem uma GUI baseada em Java.
Mike T
Existe uma maneira de ignorar / pular problemas de 'Interseção sem assentimento entre LINESTRING'? Todos os polígonos de entrada parecem estar bem (verificado com o verificador de geometria QGIS).
44520 Eclipse_by_the_moon #
1
'ST_Boundary (ST_SnapToGrid (geom, 0,001))' em vez de 'ST_Boundary (geom)' resolve o problema.
eclipsed_by_the_moon
6

Eu acho que é um pouco complicado, por causa dos diferentes conjuntos de nós dos dois polígonos (polígono verde A, segmentos diferentes vermelhos do polímero B). A comparação dos segmentos de ambos os polígonos fornece uma pista de quais segmentos do polígono B serão modificados.

Polígono de nós A

poli a

Nós do polígono B de segmentos "diferentes"

seg diff

Infelizmente, isso mostra apenas a diferença na estrutura do segmento, mas espero que seja um ponto de partida e funcione assim:

Após um processo de download e descompactação, importei o conjunto de dados usando o PostgrSQL 9.46, o PostGIS 2.1 no Debian Linux Jessie com os 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

Supondo que os segmentos do polígono A não estejam em B e vice-vera, tento criar a diferença entre os segmentos de ambos os conjuntos de polígonos, negligenciando a associação do segmento aos polígonos de cada grupo (A ou B). Por razões didáticas, formulo o material SQL em várias visualizações.

Correspondendo a este post do GIS-SE , decomponho os dois polígonos em tabelas de segmentos segments_aesegments_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;

O polígono da tabela de segmentos 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

O mesmo procedimento foi aplicado ao 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;

O polígono da tabela de segmentos 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)
...                        

Eu posso criar uma exibição de tabela de diferença chamada segments_diff_{a,b}. A diferença é dada pela não ocorrência de pontos de início ou de fim classificados no conjunto de segmentos A e 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

E o material complementar:

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

Conclusão: Para obter um resultado adequado para os pequenos pequenos segmentos marcados com a seta vermelha, os dois polígonos devem ter a mesma estrutura do nó e uma etapa de interseção no nível do nó (é necessário inserir vértices do polígono A em B). A interseção pode ser feita 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;

Mas com resultados estranhos ...

versão cortada

huckfinn
fonte
Obrigado por seus esforços. Bem, os resultados são estranhos. Gostaria de saber se ST_HausdorffDistance () pode ajudar a responder à pergunta: gis.stackexchange.com/questions/180593/…
Mar Lunar
Hum, st_haudorffdistance fornece uma medida de similaridade, não os segmentos desejados (setas vermelhas apontando para).
Huckfinn 29/02
É apenas uma ideia, ST_HausdorffDistance pode ser usado para comparar as geometrias de ambas as tabelas. Os polígonos não são iguais espacialmente? Tenho que criar linhas. Eu simplesmente não sei como fazer isso.
Mar Lunar
Parece ser uma questão de precisão e topologia ( gis.stackexchange.com/a/182838/26213 e webhelp.esri.com/arcgisdesktop/9.2/… )
huckfinn
1

Observando o exemplo, a alteração implica que os recursos da nova tabela que foram alterados sempre estarão sobrepondo-os à tabela antiga. Portanto, você seria feito com

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

A negação dos toques ocorre porque os recursos também se sobrepõem se apenas suas bordas compartilharem os mesmos locais de vértices.

inclinar
fonte