Encontre polígonos de cabeceira

8

Esta é uma pergunta de acompanhamento para esta pergunta .

Eu tenho uma rede fluvial (multilinha) e alguns polígonos de drenagem (veja a figura abaixo). Meu objetivo é selecionar apenas os polígonos da cabeceira (verde).

insira a descrição da imagem aqui

Com a solução de John, posso extrair facilmente os pontos de partida do rio (estrelas). No entanto, eu posso ter situações (polígono vermelho) em que tenho pontos de partida em um polígono, mas o polígono não é um polígono de cabeceira, porque é levado pelo rio. Eu só quero os polígonos da cabeceira.

Tentei selecioná-los contando o número de interseções entre polígonos e rios (justificativa: um polígono de cabeceira deveria ter apenas 1 interseção com o rio)

SELECT 
    polyg.*
FROM 
    polyg, start_points, stream
WHERE 
    st_contains(polyg.geom, start_points.geom)
    AND ST_Npoints(ST_Intersection(poly.geom, stream.geom)) = 1

, onde poylg são os poylgons, start_points de johns respondem e stream é minha rede fluvial.

No entanto, isso leva uma eternidade e eu não o executei:

"Nested Loop  (cost=0.00..20547115.26 rows=641247 width=3075)"
"  Join Filter: _st_contains(ezg.geom, start_points.geom)"
"  ->  Nested Loop  (cost=0.00..20264906.12 rows=327276 width=3075)"
"        Join Filter: (st_npoints(st_intersection(ezg.geom, rivers.geom)) = 1)"
"        ->  Seq Scan on ezg_2500km2_31467 ezg  (cost=0.00..2161.52 rows=1648 width=3075)"
"              Filter: ((st_area(geom) / 1000000::double precision) < 100::double precision)"
"        ->  Materialize  (cost=0.00..6364.77 rows=39718 width=318)"
"              ->  Seq Scan on stream_typ rivers  (cost=0.00..4498.18 rows=39718 width=318)"
"  ->  Index Scan using idx_river_starts on river_starts start_points  (cost=0.00..0.60 rows=1 width=32)"
"        Index Cond: (ezg.geom && geom)"

Portanto, minha pergunta é: como posso consultar com eficiência polígonos de cabeceira?

Atualização: adicionei alguns dados de amostra à minha caixa de depósito . Os dados são do sudoeste da Alemanha. São dois arquivos de forma - um com fluxos e outro com polígonos.

EDi
fonte
Portanto, para deixar claro, você deseja os polígonos que contêm apenas pontos de início, não os pontos de início. E os pontos de partida são definidos como na sua pergunta anterior (que eu respondi e até onde sei), corretamente?
John Powell
Jupp, apenas os polígonos que contêm pontos de início E não são passados ​​por um rio / são apenas começos do rio. O polígono vermelho acima contém startpoints, mas não é um polígono headwater como o rio flui através dele / não começa dentro do polígono ...
EDi
Portanto, você deseja que o conjunto polygonscontenha apenas os pontos que são fontes de rios (da pergunta anterior) e exclua os locais onde dois rios se encontram. Desculpe, por todas as perguntas, só quero ter certeza.
John Powell
Não, por exemplo, no polígono verde inferior, também dois rios se encontram. Quero excluir aqueles polygonsque têm um rio que passa (o rio entra e sai do polígono) e os mantém com início (e os rios deixam apenas esse polígono).
Edi
1
Eu não conheço nenhum PostGIS, então não posso ajudar com código direto, no entanto, no ArcGIS, eu iria ao longo destas linhas: (1) cruzar entre linhas e polígonos em um arquivo de ponto. (2) excluir (espacialmente) pontos idênticos. (3) adicione um campo numérico ao parâmetro point com o valor 1 para cada ponto. (4) junte espacialmente o polígono nos pontos e use a soma do campo numérico para indicar o tipo de drenagem. Uma soma de 1 significa que é um promontório. Superior a 1 significa que há mais de 1 entrada ou saída.
Mikkel Lydholm Rasmussen

Respostas:

4

Acredito que o esquema geral (parcialmente testado até agora) é:

  1. Encontre os pontos que representam as fontes de fluxo, como nesta resposta .

  2. Interseção com a tabela de polígonos para obter uma contagem dos vértices de origem por polígono.

  3. Use ST_DumpPoints em conjunto com o grupo por geometria para obter uma contagem de cada ponto. A idéia é contar quantos rios se encontram em um determinado ponto.

Um exemplo dessa consulta:

SELECT count(geom), ST_AsText(geom) as wkt
FROM 
   (SELECT (ST_DumpPoints(foo.geom)).geom 
   FROM 
     (SELECT 
        ST_Collect(ST_MakeLine(ST_MakePoint(0,0), ST_MakePoint(10,10)),
                   ST_MakeLine(ST_MakePoint(0,0), ST_MakePoint(20,20))
        ) as geom
     ) foo 
 ) bar 
 GROUP BY geom; 

que retorna:

count  |  wkt      
-------+--------------
 2     | POINT(0 0)
 1     | POINT(10 10)
 1     | POINT(20 20)
  1. Execute uma interseção 3contra a tabela de polígonos, para obter uma contagem (soma dos vértices) de junções de rio por polígono.

  2. Una os polígonos a partir 2de então 4, rejeitando aqueles em que a contagem (soma dos vértices) dos pontos em uma junção é maior que a soma das fontes do rio, obtida pela soma das fontes por polígono das etapas 1 e 2. Se essa condição for válida, significa que pelo menos um dos rios que se encontra em uma junção teve origem fora do polígono em questão.

Todos eles podem ser combinados em uma seqüência ampla de CTEs, a menos que algumas tabelas tenham sido criadas a partir das etapas que envolvem pontos (e indexadas).

Não tenho idéia de qual será o tempo de execução disso em um conjunto de dados completo, tendo testado apenas parte disso em um subconjunto, mas com um índice espacial na tabela de polígonos, haverá alguma assistência - obviamente não é possível aplique um índice aos pontos que emergem de ST_DumpPoints, para que uma verificação completa seja necessária lá, embora eles devam estar na memória até então.

Isso não está sendo publicado como uma resposta completa , mas como um trabalho em andamento e como uma chance de encontrar falhas na lógica. Consultas de trabalho em breve.

EDIT 1

Essa é a consulta que eu criei, que parece funcionar em um pequeno subconjunto de seus dados, mas é executada por horas no conjunto de dados completo.

CREATE TABLE good_polys as  
   WITH 
     rivers as 
       (SELECT (ST_DUMP(ST_LineMerge(geom))).geom as geom FROM streams),
     start_points as
       (SELECT ST_StartPoint(geom) as geom FROM rivers),
     end_points as 
        (SELECT ST_EndPoint(geom) as geom FROM rivers),
     junctions as 
        (SELECT (ST_DumpPoints(geom)).geom 
        FROM (SELECT geom FROM streams) s),
     source_polygons as 
        (SELECT 
            count(rivers.geom) as source_count, 
            polygons.geom, 
            polygons.gid 
         FROM rivers, polygons
         WHERE st_intersects(polygons.geom, rivers.geom) 
         GROUP BY polygons.geom, polygons.gid),
     junction_polygons as 
        (SELECT 
            count(junctions.geom) as junction_count, 
            polygons.geom, 
            polygons.gid 
         FROM junctions, polygons
         WHERE st_intersects(polygons.geom, junctions.geom) 
         GROUP BY polygons.geom, polygons.gid)
    SELECT 
       jp.gid 
    FROM 
       junction_polygons jp, source_polygons sp 
    WHERE ST_Equals(jp.geom, sp.geom) 
    AND junction_count <= source_count;

EDIT 2 . Embora isso pareça produzir respostas corretas em um pequeno subconjunto, o tempo de execução no conjunto de dados completo é horrível , presumivelmente porque a consulta final está fazendo comparações n ^ 2 e não está usando um índice espacial. A solução provável seria quebrar a consulta e criar tabelas a partir dos pontos iniciais e das consultas dos polígonos, que podem ser indexadas espacialmente antes da etapa final.

John Powell
fonte
A consulta está atualmente em execução na minha área de trabalho. Não tenho idéia de quanto tempo levará ou se estará correto, embora pareça razoável a partir de uma pequena amostra de seus dados. Você tem idéia de quantos polígonos atendem aos seus critérios?
19415 John Powell
Vou executar a consulta em um servidor. Eu acho que apenas uma pequena parte dos polígonos irá satisfazer os critérios de selecção ...
EDi
Foi o que encontrei em um subconjunto. Vou postar a minha consulta assim que terminar #
John Powell
Simplificação amanhã.
John Powell
Desculpe, muito ocupado hoje. Acho que a resposta é executar a consulta de origem e as consultas de junções de rio primeiro, cruzar com a tabela de polígonos para obter contagens por polígono, salvá-las como tabelas e indexá-las. Em seguida, execute a etapa final, onde as geometrias são iguais e comparando a contagem de pontos das duas tabelas. Espero que isso use um índice, em vez de fazer comparações n² como presentes. Voltará mais tarde.
19415 John Powell
0

No pseudo-código, isso deve funcionar:

  • selecionar todos os polígonos
  • (EXTERIOR CHEIO?), Junte-se a pontos no polígono que cruza pontos
  • (COMPLETO EXTERIOR?) Une linhas em que o polígono cruza linhas
  • foram line.riverid não é igual a point.riverid
  • agrupar por polygonid
  • count (pointid)> 0

Não tenho muita certeza de como criar a consulta e não posso testá-la sem um banco de dados para testar. É uma pergunta bem louca, eu acho. Mas deve funcionar!

Alex Leith
fonte