PostGIS: ST_Equals false quando ST_Intersection = 100% da geometria?

9

Eu tenho 2 conjuntos de dados que consistem em dados de parcelas cadastrais - aproximadamente 125.000 linhas cada. A coluna de geometria é polígonos WKB que representam os limites das parcelas; todos os dados são geometricamente válidos (os polígonos estão fechados etc.).

Alguns dados recentes chegaram em uma projeção diferente dos dados de base usados ​​para um trabalho de comparação - então eu reprojetei o mais recente (a base era 4326; a outra era o WGA94 que foi trazido para o PostGIS como 900914 ... eu o reprojetei para 4326) .

A primeira etapa da análise foi encontrar e armazenar parcelas não correspondentes; parte disso é identificar e armazenar parcelas com geometria idêntica.

Então, eu executei uma consulta muito padrão (o bloco de código abaixo abstrai os detalhes do esquema, etc.):

create table matchdata as
  select  a.*
  from gg2014 a, gg2013 b
  where ST_Equals(a.g1,b.g1)

ZERO results.

"Estranho ..." pensei. "Talvez tenha havido pequenas mudanças de vértices causadas pela reprojeção: isso seria irritante e realmente não deveria acontecer".

Felizmente, existem dados aspaciais abundantes (5 colunas identificadoras) que me permitem estabelecer parcelas que devem ser espacialmente idênticas: aquelas com o mesmo identificador, cuja data de alteração na tabela de 2014 era anterior à data de alteração máxima nos dados de 2013. Isso totalizou 120.086 linhas distintas.

Armazenei os identificadores e geometrias em uma tabela separada ( match_id) e executei a seguinte consulta:

select apid, 
       bpid, 
       ST_Area(ag::geometry) as aa, 
       ST_Area(bg::geometry) as ab,
       ST_Area(ST_Intersection(ag,bg)::geometry)/ST_Area(ag::geometry) as inta,
       ST_Area(ST_Intersection(ag,bg)::geometry)/ST_Area(ag::geometry) as intb
from match_id
order by inta

Os 16 primeiros valores intae intbforam idênticos a zero, os 456 seguintes foram 0,99999999-ish (min 0,99999999999994, máx 0,999999999999999) e as linhas 473 em diante foram 1 - até a linha 120050, quando a área da interseção era maior que a geometria (a maior valor para intae intbera 1.00000000000029, mas ainda assim).

Então, eis o meu enigma: se duas geometrias se cruzam espacialmente entre 99,999999999994% e 100,00000000002929% de suas respectivas áreas, eu gostaria que "ST_Equals" dissesse "Sim ... eu darei a você essa. Fechar o suficiente".

Afinal, é equivalente a ficar em 1 parte em 16 trilhões ... ou seja, como se a dívida nacional dos EUA estivesse abaixo de 93 centavos.

No contexto da circunferência da Terra (a ~ 40.000 km), é como estar em 0,0000000025 km, no máximo (desde que resultem em uma diferença de área tão pequena, qualquer mudança de vértice deve ser ainda menor).

De acordo com o TFD (que eu recomendei), a tolerância ST_Intersects()é de 0,00001m (1mm), de modo que as alterações implícitas nos vértices (que eu confesso que não verifiquei: quero ST_Dump()e faço) pareceriam menores do que a tolerância. (Eu percebo isso ST_Intersects !== ST_Intersection(), mas é a única tolerância mencionada).

Não fui capaz de descobrir a tolerância correspondente à comparação de vértices realizada por ST_Equals()... mas parece realmente estranho que pelo menos 120.000 das minhas linhas devam passar por qualquer avaliação sensata da identidade espacial, mas não o fazem.

(Nota: eu também fiz o mesmo exercício usando ::geography- com resultados que tinham mais variabilidade, mas ainda mais que 110.000 entradas com um bom '1' limpo).

Existe uma maneira de diminuir a tolerância de ST_Equals, que não exija escavação nos interstícios do código? Eu não estou interessado em fazer isso.

Se não, existe algum argumento que alguém esteja ciente?

Nota: seria bom se o 'kludge' não estivesse fazendo uma comparação bilateral como

where ST_within(g1, ST_Buffer(g2, 0.0000001))
  and ST_within(g2, ST_Buffer(g1, 0.0000001))


   - I've done that: sure, it works... but it's a gigantic documentation PITA).

Posso resolver isso, mas escrever as 20 páginas para documentar a solução alternativa - que só será exibida novamente se tivermos dados duvidosos - é uma PITA que eu preferiria não fazer, já que é provável que seja uma ocorrência única .

(Versões: Postgresql 9.3.5; PostGIS 2.1.3)

GT.
fonte
Apenas um pensamento aqui, mas você tentou canonizar as novas parcelas em uma grade que seja compatível com os dados existentes usando st_snaptogrid?
nickves
Entendo que não quero olhar para o código-fonte, mas sua pergunta me levou a fazê-lo (mesmo que meu C ++ seja uma porcaria), por isso agradeço por isso. Se você estiver interessado, posso postar as seções relevantes, todas em github.com/libgeos .
John Powell
ST_Equalsretorna apenas truequando as geometrias são iguais - tipo de geometria, número de vértices, SRID e valores de vértice (em todas as dimensões, na mesma ordem). Se houver alguma variação, a comparação será interrompida e falseretornada.
Vince
@Vince: como eu o entendo (a partir dos documentos), ST_Equals()ignora a direcionalidade. Entendi que isso significa que, para um polígono 2-D fechado, não faz diferença se os pontos são enumerados no sentido horário versus anti-horário. ST_OrderingEquals()é o teste mais rígido. Dito isso, depois de inspecionar os vértices (usando ST_Dump()e calculando deltas para todos os vértices), fica claro que a incrível resposta de @John Barça está no dinheiro. ST_equals()é contra-indicado, mesmo para dados idênticos conhecidos ex-ante , se uma geometria for reprojetada - a menos que a comparação seja feita com ST_SnapToGrid ().
GT.
Voltando a isso: uma maneira rápida e agradável de obter um teste aceitável de [quase] igualdade espacial é verificar qual proporção de cada geometria faz parte da interseção. É um pouco computacionalmente oneroso; calcular (100*(ST_Area(ST_Intersection(a.g1, b.g1))/ST_Area(a.g1)))::int as int_pcae (100*(ST_Area(ST_Intersection(a.g1, b.g1))/ST_Area(b.g1)))::int as int_pcb(verifique se você JOINinclui ST_Intersects(a.g1,b.g1)). Teste se (int_pca, int_pcb)=(100,100)(ou algum outro conjunto de pontos de corte). Kludgy, mas fará 2,6 milhões de encomendas em ~ 30 min (desde que g1 esteja indexado no GIST).
GT.

Respostas:

20

Meu palpite é que você coordenou transformações introduziram pequenos erros de arredondamento (veja um exemplo abaixo). Como não há como definir a tolerância em ST_Equals, isso faz com que ST_Equals retorne false para algumas geometrias que diferem apenas na enésima casa decimal, pois as geometrias precisam ser idênticas em todos os aspectos - consulte a definição da matriz de interseção nas libgeos . Você pode verificar isso com um exemplo realmente extremo,

SELECT ST_Equals(
      ST_MakePoint(0,0),
      ST_MakePoint(0,0.000000000000000000000000000000000000000000000000000000000001));

que retorna falso .

Se você usa ST_SnapToGrid, pode impor uma precisão específica, por exemplo, a dez casas decimais,

SELECT ST_Equals(
      ST_MakePoint(0,0),
      ST_SnapToGrid(
             ST_MakePoint(0,0.00000000000000000000000000000000000000000000001),
      10));

agora retorna verdadeiro .

Se você fosse correr,

CREATE table matchdata AS
SELECT  a.*
FROM gg2014 a, gg2013 b
WHERE ST_Equals(ST_SnapToGrid(a.g1, 5), ST_SnapToGrid(b.g1, 5));

definindo uma tolerância apropriada, suspeito que seus problemas desapareceriam.

Aqui está um link para uma discussão do desenvolvedor do Postgis sobre tolerância, o que sugere que é menos do que trivial de implementar.

Fiz algumas conversões entre a British National Grid (EPSG: 27700) e a lat / lon para ilustrar o ponto sobre arredondar a precisão, tomando um ponto em algum lugar de Londres,

SELECT ST_AsText(ST_Transform(ST_SetSrid(ST_MakePoint(525000, 190000),27700),4326));

retorna POINT(-0.19680497282746 51.5949871603888)

e revertendo isso,

SELECT ST_AsText(ST_Transform(ST_SetSrid(ST_MakePoint(-0.19680497282746, 51.5949871603888),4326),27700));

retorna POINT(525000.000880007 189999.999516211)

que está desativado em menos de um milímetro, mas mais que o suficiente para fazer ST_Equals retornar falso.

John Powell
fonte
A resposta de John Barca está correta - que pequenos erros de arredondamento podem prejudicar o ST_Equals. Minha experiência (desagradável) foi ao trabalhar com dois conjuntos de dados - ambos projetados do EPSG 4326 ao EPSG 3857 - um via ArcCatalog (ArcToolbox -> Ferramentas de gerenciamento de dados -> projeções e transformações) , enquanto o outro via GDAL ogr2ogr.
Ralph Tee
Esta resposta me ajudou muito. Mas notei que os índices geográficos não eram mais usados ​​e as consultas demoravam muito. Minha solução alternativa foi criar tabelas temporárias com as geometrias ajustadas e adicionar índices a elas antes de executar a consulta. Existe uma maneira melhor de acelerar as coisas?
HFS
1
@hfs. Eu acredito que você pode criar um índice funcional usando ST_SnapToGrid. Você está correto que o uso de uma chamada de função dentro de uma operação espacial igual a / cruza / contém etc fará com que o índice não seja usado e criar um índice funcional resolveria isso. Ou você pode atualizar permanentemente seus dados se achar que a precisão é espúria e não precisar usar ST_SnapToGrid na consulta. Depende dos seus dados e caso de uso, é claro.
John Powell
2

Você executou a verificação ST_IsValid nas suas geometrias? Se eles são inválidos, todas as apostas estão desativadas. ST_Intersects e a outra família de funções de relacionamento espacial do GEOS geralmente retornam false porque a área não está bem definida do ponto de vista da matriz de interseção. O motivo de fazer ST_Buffer provavelmente funcionar é porque está convertendo suas geometrias inválidas em válidas. ST_Buffer (..., tinybit) é o que é conhecido como uma ferramenta "homem pobre tenta tornar minha geometria válida".

LR1234567
fonte
A primeira etapa de qualquer novo conjunto de dados é selecionar apenas geometrias válidas usando ST_isValid(g1)- isso foi mencionado (obliquamente) "[a] coluna de geometria são polígonos WKB representando limites de parcelas; todos os dados são geometricamente válidos (os polígonos estão fechados etc.) ".
GT.
0

Minha resposta chega um pouco tarde, mas talvez ajude alguém que tenha o mesmo problema. Da minha experiência, quando duas geometrias que são realmente iguais, mas ST_Equals retorna False, duas coisas podem ajudar:

  1. verifique se as geometrias comparadas são geometrias únicas (sem MultiLinesting, MultiPoin etc.)
  2. tente em ST_Equals(st_astext(a.geom), st_astext(b.geom)) vez deST_Equals(a.geom , b.geom)

O primeiro já é mencionado na documentação . O segundo parece irracional, mas funciona. Não sei, mas acho que tem a ver com o formato binário da geometria padrão do postGIS.

ioanna tsak
fonte