Usando ST_Difference para remover recursos sobrepostos?

11

Estou tentando usar ST_Difference para criar um conjunto de polígonos (processing.trimmedparcelsnew) que não contêm nenhuma área coberta por outro conjunto de polígonos (test.single_geometry_1) usando o PostGis 2.1 (e o Postgres SQL 9.3). Aqui está a minha consulta:

CREATE TABLE processing.trimmedparcelsnew AS
SELECT
    orig.id, ST_Difference(orig.geom, cont.geom) AS difference
FROM 
    test.single_geometry_1 cont,
    test.multi_geometry_1 orig;

Mas os polígonos resultantes não foram aparados; eles parecem ter sido divididos onde se cruzam com a outra camada. Tentei apenas executar o select sem colocar o resultado em uma tabela e tudo o mais em que consigo pensar, mas não consigo fazer essa função funcionar.

Anexei uma foto do resultado

insira a descrição da imagem aqui


Após os comentários, tentei adicionar uma cláusula WHERE. Quero que as parcelas que não têm interseções e as áreas de interseção das outras parcelas sejam removidas (a camada test.single_geometry representa a contaminação que eu quero remover das minhas parcelas). Eu tentei uma interseção, mas é claro que eu realmente quero as não interseções, então agora estou tentando uma desconexão. Também tentei adicionar o orig à minha tabela, mas a documentação para ST_Difference ( http://postgis.net/docs/ST_Difference.html ) diz que ele retorna a geometria exata que eu preciso (uma geometria que representa essa parte da geometria A que não se cruza com a geometria B); portanto, estou confuso sobre o motivo pelo qual desejaria o polígono original na minha tabela. Enfim, aqui está o meu código modificado:

CREATE TABLE processing.trimmedparcelsnew AS
SELECT
    orig.id, ST_Difference(orig.geom, cont.geom) AS difference, orig.geom AS geom
FROM 
    test.single_geometry_1 cont,
    test.multi_geometry_1 orig
WHERE ST_Disjoint(orig.geom, cont.geom);

Seguindo a resposta do dbaston, tentei agora:

CREATE TABLE processing.parcels_trimmed AS
SELECT id, COALESCE(ST_Difference(geom, (SELECT ST_Union(b.geom) 
                                         FROM test.single_geometry_1 b
                                         WHERE ST_Intersects(a.geom, b.geom)
                                         AND a.id != b.id)), a.geom)
FROM test.multi_geometry_1 a;

O resultado disso é apenas uma cópia do test.multi_geometry_1. Embora agora a divisão não esteja mais ocorrendo.

Eu tentei a versão anterior, mas novamente é só pegar uma cópia do test.multi_geometry_1:

CREATE TABLE processing.parcels_trimmed_no_coalesce AS
SELECT id, COALESCE(ST_Difference(geom, (SELECT ST_Union(b.geom) 
                                         FROM test.single_geometry_1 b
                                         WHERE ST_Intersects(a.geom, b.geom)
                                         AND a.id != b.id)), a.geom)
FROM test.multi_geometry_1 a;

Estou começando a me perguntar se há algo mais que estou fazendo errado? A declaração de processo é:

DROP TABLE IF EXISTS processing.parcels_trimmed_no_coalesce;

E estou executando as consultas na janela de consulta SQL do PostgreSQL e no Openjump.

A declaração que eu uso para ver a tabela é:

SELECT * FROM processing.parcels_trimmed_no_coalesce;

No interesse da simplificação, reduzi agora essa consulta para apenas:

SELECT id, COALESCE(ST_Difference(geom, (SELECT ST_Union(b.geom) 
                                         FROM test.geometriestocutagainst b
                                         WHERE ST_Intersects(a.geom, b.geom)
                                         AND a.id != b.id)), a.geom)
FROM test.geometriestocut a;

Isso ainda resulta apenas nos polígonos originais (teste.geometriestocut) quando o resultado desejado é o original aparado em comparação com teste.geometriestocut.

Mart
fonte
Você não especificou uma WHEREcláusula, portanto, pode haver expansão polinomial na tabela resultante. Quantas linhas estão trimmedparcelsnew?
Vince
Se você quiser apenas a diferença de onde eles se cruzam, tente adicionar WHERE ST_Intersects (orig.geom, cont.geom). Caso contrário, a diferença de dois polígonos que não se cruzam é ​​o polígono original.
John Powell
Há 24 linhas novas em parcelas aparadas. Quero a diferença mesmo quando elas não se cruzam. Por isso, corrijo que preciso usar o orig.geom na tabela em vez da diferença?
24715 Mart
Um teste disjuntos deve produzir expansão polinomial - cada recurso em cont a aparecer uma vez para cada recurso no orig que não se sobrepõem, ea diferença não alterar uma geometria de entrada
Vince
Ok, obrigado pelo esclarecimento, mas ainda não tenho certeza porque o código original não está funcionando. Se ST_Difference (orig.geom, cont.geom) retorna as geometrias em a que não se cruzam com b, por que a tabela contém as geometrias divididas em vez das geometrias em que não se cruzam b.
Mart

Respostas:

14

Uma união automática permite que você opere no relacionamento entre pares de dois recursos. Mas não acho que você esteja interessado em pares: para cada recurso, você deseja operar o relacionamento entre esse recurso e todos os outros recursos do seu conjunto de dados. Você pode fazer isso com uma expressão de subconsulta:

CREATE TABLE parcels_trimmed AS
SELECT id, ST_Difference(geom, (SELECT ST_Union(b.geom) 
                                FROM parcels b
                                WHERE ST_Intersects(a.geom, b.geom)
                                  AND a.id != b.id))
FROM parcels a;

Você pode ver algo estranho nos resultados, no entanto. As encomendas que não possuem sobreposições estão sendo totalmente descartadas! Isso ST_Unionocorre porque o agregado em um conjunto de registros vazio será NULLe ST_Difference(geom, NULL)é NULL. Para suavizar isso, você precisa encerrar sua ST_Differencechamada em COALESCE:

CREATE TABLE parcels_trimmed AS
SELECT id, COALESCE(ST_Difference(geom, (SELECT ST_Union(b.geom) 
                                         FROM parcels b
                                         WHERE ST_Intersects(a.geom, b.geom)
                                         AND a.id != b.id)), a.geom)
FROM parcels a;

Isso significa que, se o resultado de ST_Differencefor NULL, a expressão coalescente será avaliada para a geometria original.

A consulta acima removerá totalmente as áreas sobrepostas do seu domínio. Se preferir escolher um vencedor, você pode fazer a.id < b.id, ou algum outro critério, em vez de a.id != b.id.

dbaston
fonte
Obrigado pela resposta, infelizmente estou tendo problemas para fazer isso funcionar para mim. Em vez disso, acabe com o polígono original (a). Vou editar minha pergunta com mais informações. Obrigado novamente.
Mart
2

Eu tive o mesmo problema que você. Não sei se você já encontrou uma solução para o seu problema, mas modifiquei a resposta aceita acima e consegui o que queria.

CREATE TABLE parcels_trimmed AS
SELECT id, COALESCE(ST_Difference(geom, (SELECT ST_Collect(b.geom) 
                                         FROM parcels b
                                         WHERE ST_Intersects(a.geom, b.geom)
                                         )), a.geom)
FROM parcels a;
Dag
fonte
1

Eu uso ST_DifferenceAgg () dos PostGIS Addons . Você precisa mesclar as duas tabelas, ter um identificador exclusivo e um índice na coluna de geometria. Aqui está um pequeno exemplo:

WITH overlappingtable AS (
  SELECT 1 id, ST_GeomFromText('POLYGON((0 1, 3 2, 3 0, 0 1), (1.5 1.333, 2 1.333, 2 0.666, 1.5 0.666, 1.5 1.333))') geom
  UNION ALL
  SELECT 2 id, ST_GeomFromText('POLYGON((1 1, 3.8 2, 4 0, 1 1))')
  UNION ALL
  SELECT 3 id, ST_GeomFromText('POLYGON((2 1, 4.6 2, 5 0, 2 1))')
  UNION ALL
  SELECT 4 id, ST_GeomFromText('POLYGON((3 1, 5.4 2, 6 0, 3 1))')
  UNION ALL
  SELECT 5 id, ST_GeomFromText('POLYGON((3 1, 5.4 2, 6 0, 3 1))')
)
SELECT a.id, ST_DifferenceAgg(a.geom, b.geom) geom
FROM overlappingtable a,
     overlappingtable b
WHERE a.id = b.id OR -- Make sure to pass at least once the polygon with itself
      ((ST_Contains(a.geom, b.geom) OR -- Select all the containing, contained and overlapping polygons
        ST_Contains(b.geom, a.geom) OR
        ST_Overlaps(a.geom, b.geom)) AND
       (ST_Area(a.geom) < ST_Area(b.geom) OR -- Make sure bigger polygons are removed from smaller ones
        (ST_Area(a.geom) = ST_Area(b.geom) AND -- If areas are equal, arbitrarily remove one from the other but in a determined order so it's not done twice.
         a.id < b.id)))
GROUP BY a.id
HAVING ST_Area(ST_DifferenceAgg(a.geom, b.geom)) > 0 AND NOT ST_IsEmpty(ST_DifferenceAgg(a.geom, b.geom));

Isso mesclará as partes sobrepostas com o maior polígono sobreposto. Se você deseja manter a parte sobreposta separada, veja o exemplo ST_splitAgg ().

Pierre Racine
fonte