Qual a melhor forma de corrigir um problema de interseção sem aceno no PostGIS?

38

Estou usando uma PL/Rfunção e PostGISpara gerar polígonos voronoi em torno de um conjunto de pontos. A função que estou usando é definida aqui . Quando uso essa função em um conjunto de dados específico, recebo a seguinte mensagem de erro:

Error : ERROR:  R interpreter expression evaluation error
DETAIL:  Error in pg.spi.exec(sprintf("SELECT %3$s AS id,   
st_intersection('SRID='||st_srid(%2$s)||';%4$s'::text,'%5$s') 
AS polygon FROM %1$s WHERE st_intersects(%2$s::text,'SRID='||st_srid(%2$s)||';%4$s');",  
:error in SQL statement : Error performing intersection: TopologyException: found non-noded 
intersection between LINESTRING (571304 310990, 568465 264611) and LINESTRING (568465 
264611, 594406 286813) at 568465.05533706467 264610.82749605528
CONTEXT:  In R support function pg.spi.exec In PL/R function r_voronoi

Ao examinar esta parte da mensagem de erro:

Error performing intersection: TopologyException: found non-noded intersection between
LINESTRING (571304 310990, 568465 264611) and LINESTRING (568465 264611, 594406 286813) 
at 568465.05533706467 264610.82749605528

É assim que o problema listado acima se parece:

insira a descrição da imagem aqui

Inicialmente, pensei que essa mensagem pudesse ser causada pela existência de pontos idênticos e tentei resolver isso usando a st_translate()função, usada da seguinte maneira:

ST_Translate(geom, random()*20, random()*20) as geom 

Isso resolve o problema, mas minha preocupação é que agora estou traduzindo todos os pontos até ~ 20m na ​​direção x / y. Também não sei dizer qual é a quantidade necessária de tradução. Por exemplo, neste conjunto de dados por tentativa e erro a 20m * random numberestá ok, mas como posso saber se isso precisa ser maior?

Com base na imagem acima, acho que o problema é que o ponto está cruzando com a linha enquanto o algoritmo está tentando cruzar o ponto com um polígono. Não tenho certeza do que devo fazer para garantir que o ponto esteja dentro de um polígono, em vez de cruzar com uma linha. O erro está ocorrendo nesta linha:

"SELECT 
  %3$s AS id, 
  st_intersection(''SRID=''||st_srid(%2$s)||'';%4$s''::text,''%5$s'') AS polygon 
FROM 
  %1$s 
WHERE 
  st_intersects(%2$s::text,''SRID=''||st_srid(%2$s)||'';%4$s'');"

Eu li essa pergunta anterior: o que é um "cruzamento sem aceno"? para tentar entender melhor esse problema e gostaria de receber conselhos sobre a melhor maneira de resolvê-lo.

djq
fonte
Se suas entradas não são válidas para começar, execute ST_MakeValid () nelas. Se eles são válidos, adicionar entropia, como você está fazendo, é o próximo truque disponível, e talvez o último truque por enquanto.
Paul Ramsey
Sim, estou usando WHERE ST_IsValid(p.geom)para filtrar os pontos inicialmente.
DJQ

Respostas:

30

Na minha experiência, esse problema quase sempre é causado por:

  1. Alta precisão em suas coordenadas (43.231499999999996), combinada com
  2. Linhas quase coincidentes, mas não idênticas

A abordagem "cutucada" das ST_Buffersoluções permite que você se safe da segunda posição, mas tudo o que você pode fazer para resolver essas causas subjacentes, como encaixar sua geometria em uma grade 1e-6, facilitará sua vida. As geometrias com buffer geralmente são boas para cálculos intermediários, como a área de sobreposição, mas você deve ter cuidado ao retê-las, pois elas podem piorar seus problemas próximos, mas não muito, a longo prazo.

O recurso de tratamento de exceções do PostgreSQL permite escrever funções de invólucro para lidar com esses casos especiais, armazenando em buffer somente quando necessário. Aqui está um exemplo para ST_Intersection; Eu uso uma função semelhante para ST_Difference. Você precisará decidir se o buffer e o retorno potencial de um polígono vazio são aceitáveis ​​em sua situação.

CREATE OR REPLACE FUNCTION safe_isect(geom_a geometry, geom_b geometry)
RETURNS geometry AS
$$
BEGIN
    RETURN ST_Intersection(geom_a, geom_b);
    EXCEPTION
        WHEN OTHERS THEN
            BEGIN
                RETURN ST_Intersection(ST_Buffer(geom_a, 0.0000001), ST_Buffer(geom_b, 0.0000001));
                EXCEPTION
                    WHEN OTHERS THEN
                        RETURN ST_GeomFromText('POLYGON EMPTY');
    END;
END
$$
LANGUAGE 'plpgsql' STABLE STRICT;

Outro benefício dessa abordagem é que você pode identificar as geometrias que realmente estão causando seus problemas; basta adicionar algumas RAISE NOTICEinstruções no EXCEPTIONbloco para gerar o WKT ou outra coisa que o ajudará a rastrear o problema.

dbaston
fonte
Essa é uma ideia inteligente. Eu sempre descobri que os problemas de interseção vêm de cadeias de linhas que aparecem durante combinações de uniões, diferenças, buffers, etc., que podem ser corrigidos armazenando tudo em buffer ou descartando tudo e selecionando somente polígonos / mutlipolígonos. Esta é uma abordagem interessante.
John Powell
Você mencionou encaixar a geometria na grade 1e-6, mas estou sentado aqui imaginando se encaixar na potência 2 seria melhor. PostGIS (e GEOS) usando números de ponto flutuante; portanto, ajustar para uma potência de 10 pode não realmente truncar muito as coordenadas, pois o número pode não ter uma representação binária de comprimento finito. Mas se você digitar 2 ^ -16, acredito que seria truncado qualquer parte fracionária para apenas 2 bytes. Ou estou pensando errado?
jpmc26 17/01
12

Através de muitas tentativas e erros, finalmente percebi que o non-noded intersectionresultado era de um problema de auto-interseção. Eu encontrei um tópico que sugeria o uso ST_buffer(geom, 0)pode ser usado para corrigir o problema (embora o torne muito mais lento no geral). Eu tentei usar ST_MakeValid()e, quando aplicado diretamente à geometria, antes de qualquer outra função. Isso parece resolver o problema de maneira robusta.

ipoint <- pg.spi.exec(
        sprintf(
            "SELECT 
                    %3$s AS id, 
                    st_intersection(ST_MakeValid(''SRID=''||st_srid(%2$s)||'';%4$s''::text), ST_MakeValid(''%5$s'', 0)) AS polygon 
            FROM %1$s 
            WHERE 
                ST_Intersects(ST_MakeValid(%2$s::text),ST_MakeValid(''SRID=''||st_srid(%2$s)||'';%4$s''));",
            arg1,
            arg2,
            arg3,
            curpoly,
            buffer_set$ewkb[1]
        )
    )

Marquei isso como a resposta, pois parece ser a única abordagem que resolve meu problema.

djq
fonte
11

Corri para o mesmo problema (Postgres 9.1.4, PostGIS 2.1.1), e a única coisa que funcionou para mim foi envolver a geometria com um buffer muito pequeno.

SELECT ST_Intersection(
    (SELECT geom FROM table1), ST_Union(ST_Buffer(geom, 0.0000001))
) FROM table2

ST_MakeValidnão funcionou para mim, nem a combinação de ST_Nodee ST_Dump. O buffer parece não resultar em nenhuma degradação no desempenho, mas se eu o diminuísse ainda recebia um erro de interseção sem aceno.

Feio, mas funciona.

Atualizar:

A estratégia ST_Buffer parece funcionar bem, mas encontrei um problema em que ela produzia erros ao converter a geometria para a geografia. Por exemplo, se um ponto está originalmente em -90.0 e é armazenado em buffer por 0.0000001, agora está em -90.0000001, que é uma geografia inválida.

Isso significava que, embora tenha ST_IsValid(geom)sido t, ST_Area(geom::geography)retornou NaNpara muitos recursos.

Para evitar o problema de interseção sem aceno, mantendo a geografia válida, acabei usando ST_SnapToGridassim

SELECT ST_Union(ST_MakeValid(ST_SnapToGrid(geom, 0.0001))) AS geom, common_id
    FROM table
    GROUP BY common_id;
jczaplew
fonte
6

No postgis, o ST_Node deve quebrar uma série de linhas nos cruzamentos, o que deve resolver o problema de cruzamento sem aceno. Quebrar isso em ST_Dump gera a matriz composta das linhas quebradas.

Ligeiramente relacionado, há uma apresentação impressionante PostGIS: Dicas para Usuários Avançados, que descreve claramente esse tipo de problemas e soluções.

WolfOdrade
fonte
Essa é uma ótima apresentação (obrigado @PaulRamsey). Como devo usar ST_Nodee ST_Dump? Imagino que seria necessário para usá-los perto esta parte da função, mas não estou certo: st_intersection(''SRID=''||st_srid(%2$s)||'';%4$s''::text,''%5$s'')em
DJQ
Hmmm, eu não notei que as duas linhas tinham uma coordenada idêntica, o que deveria ser bom. Se você traçar essas coordenadas, o ponto de interseção fica a cerca de 18 cm da interseção. Não é realmente uma solução, apenas uma observação.
WolfOdrade
Não está totalmente claro como eu uso st_nodeaqui - posso usá-lo antes st_intersection?
DJQ
1
A apresentação não está mais disponível. Estou preso com o mesmo problema, quando se tenta ST_Clip (rast, polígono)
Jackie
1
@ Jackie: Corrigi o link para a apresentação na resposta: PostGIS: Dicas para usuários avançados .
Pete
1

Na minha experiência, resolvi meu non-noded intersectionerro usando a função St_SnapToGrid , que resolveu o problema de ter uma alta precisão nas coordenadas do vértice dos polígonos.

SELECT dissolve.machine, dissolve.geom FROM (
        SELECT machine, (ST_Dump(ST_Union(ST_MakeValid(ST_SnapToGrid(geom,0.000001))))).geom 
        FROM cutover_automatique
        GROUP BY machine) as dissolve
WHERE ST_isvalid(dissolve.geom)='t' AND ST_GeometryType(dissolve.geom) = 'ST_Polygon';
CptGasse
fonte