A maneira mais rápida de juntar espacialmente um CSV de ponto com um polígono Shapefile

19

Eu tenho um arquivo CSV de 1 bilhão de pontos e um shapefile com cerca de 5.000 polígonos. Qual seria a maneira mais rápida de unir espacialmente pontos e polígonos? Para cada ponto, preciso obter o ID do polígono. (Os polígonos não se sobrepõem.)

Normalmente, eu carregava os dois conjuntos de dados no PostGIS. Existe uma maneira mais rápida de realizar o trabalho?

Estou procurando uma solução de código aberto.

open-source-gis optimization point-in-polygon underdark
fonte

Respostas:

16

Se "mais rápido" inclui o montante de seu tempo que é gasto, a solução vai depender do que o software que você está confortável com e pode usar de forma expedita. As observações a seguir, conseqüentemente, se concentram em idéias para alcançar os tempos de computação mais rápidos possíveis .

Se você usar um programa fixo, quase certamente o melhor que poderá fazer é pré-processar os polígonos para configurar uma estrutura de dados point-in-polygon, como uma árvore KD ou quadtree, cujo desempenho normalmente será O (log (V ) * (N + V)) onde V é o número total de vértices nos polígonos e N é o número de pontos, porque a estrutura de dados exigirá pelo menos O (log (V) * V) esforço para criar e, em seguida, deve ser sondado para cada ponto a um custo por ponto O (log (V)).

Você pode obter um desempenho substancialmente melhor primeiro agrupando os polígonos, explorando a suposição de que não há sobreposições. Cada célula da grade está inteiramente no interior de um polígono (incluindo o interior do "polígono universal"); nesse caso, rotule a célula com o ID do polígono ou então contenha uma ou mais arestas de polígono. O custo dessa rasterização, igual ao número de células de grade referenciadas durante a rasterização de todas as arestas, é O (V / c), em que c é o tamanho de uma célula, mas a constante implícita na notação O-grande é pequena.

(Uma vantagem dessa abordagem é que você pode explorar rotinas gráficas padrão. Por exemplo, se você tiver um sistema que (a) desenhe os polígonos em uma tela virtual usando (b) uma cor distinta para cada polígono e (c) permita para ler a cor de qualquer pixel que você queira endereçar, você o criou.)

Com essa grade instalada, faça uma pré-seleção dos pontos computando a célula que contém cada ponto (uma operação O (1) exigindo apenas alguns relógios). A menos que os pontos estejam agrupados em torno dos limites do polígono, isso normalmente deixará apenas cerca de O (c) pontos com resultados ambíguos. O custo total de construção da grade e pré-triagem é, portanto, O (V / c + 1 / c ^ 2) + O (N). É necessário usar outro método (como qualquer um dos recomendados até agora) para processar os pontos restantes (ou seja, aqueles que estão próximos aos limites dos polígonos), a um custo de O (log (V) * N * c) .

À medida que c diminui, cada vez menos pontos estarão na mesma célula da grade com uma aresta e, portanto, cada vez menos exigirá o processamento subsequente de O (log (V)). Agindo contra isso é a necessidade de armazenar células da grade O (1 / c ^ 2) e gastar tempo O (V / c + 1 / c ^ 2) rasterizando os polígonos. Portanto, haverá um tamanho de grade ideal c. Utilizando-o, o custo computacional total é O (log (V) * N), mas a constante implícita é tipicamente muito menor do que o uso de procedimentos em lata, devido à velocidade O (N) da pré-triagem.

Há 20 anos, testei essa abordagem (usando pontos espaçados uniformemente em toda a Inglaterra e no exterior e explorando uma grade relativamente bruta de cerca de 400 mil células oferecidas pelos buffers de vídeo da época) e obtive duas ordens de aceleração de magnitude em comparação com o melhor algoritmo publicado que pude encontrar. Mesmo quando os polígonos são pequenos e simples (como triângulos), você tem praticamente a garantia de uma ordem de magnitude acelerada.

Na minha experiência, o cálculo foi tão rápido que toda a operação foi limitada pelas velocidades de E / S de dados, não pela CPU. Antecipando que a E / S pode ser o gargalo, você obteria os resultados mais rápidos armazenando os pontos no formato mais compactado possível para minimizar os tempos de leitura dos dados. Pense também em como os resultados devem ser armazenados, para que você possa limitar as gravações em disco.

whuber
fonte
6
Muito bom tempo gasto na realização da solução versus tempo de computação. Levar muito tempo para chegar a uma solução ótima só é benéfico se você realizar essas economias por meio da otimização (especialmente do ponto de vista do empregador).
Sasa Ivetic
5

Da minha parte, provavelmente carregaria dados CSV em um arquivo shp e, em seguida, escreveria um script python usando shapefile e shapely para obter o ID do polígono e atualizar o valor do campo.

Não sei se geotools e JTS são mais rápidos que shapefile / shapely ... Não tenho tempo para testá-lo!

edit : A propósito, a conversão de csv para o formato shapefile provavelmente não é necessária, pois os valores podem ser facilmente formatados para serem testados com objetos espaciais do shapefile do polígono.

simo
fonte
4
Eu carregava os dados diretamente usando um leitor csv e preenchia um índice espacial Rtree . A combinação de Rtree e Shapely tem um desempenho impressionante (muito melhor que o PostGIS; não posso comparar com o JTS porque não conheço Java).
Mike T
2
Boa ideia, desde que você não precise armazenar todos os pontos 1b na memória de uma só vez. Com um mínimo de 16 bytes por ponto (X / Y), você está analisando dados de 16 GB. Se a Rtree criar o índice no armazenamento local, definitivamente melhorará o desempenho. A importação de pontos 1b para um único arquivo de forma também não funcionará. Os shapefiles das especificações do OGR são limitados a 8 GB (recomenda-se 4 GB). Uma forma de ponto único usa 20 bytes.
Sasa Ivetic
4

Acabei convertendo os polígonos em uma varredura e amostrando-os nas posições dos pontos. Como meus polígonos não se sobrepunham e a alta precisão não era necessária (polígonos representavam classes de uso da terra e suas fronteiras eram consideradas bastante incertas de qualquer maneira), essa era a solução mais eficiente em termos de tempo que eu consegui encontrar.

underdark
fonte
3

Eu rapidamente escrever um pequeno programa Java baseado no leitor shapefile de GeoTools ea operação contém de JTS . Eu não sei o quão rápido pode ser ...

julien
fonte
1
Se você possui os dados no PostGIS, o GeoTools pode usar índices de essência etc.
Ian Turton
3

Use Spatialite .

Faça o download da GUI. Você pode abrir o Shapefile e o CSV como tabelas virtuais. Isso significa que na verdade você não os importa para o banco de dados, mas eles aparecem como tabelas e você pode ingressar rapidamente e consultá-los da maneira que desejar.

Sean
fonte
3

Você pode fazer isso rapidamente usando OGR em C / C ++ / Python (o Python deve ser o mais lento dos 3). Faça um loop em todos os polígonos e defina um filtro nos pontos, faça um loop nos pontos filtrados e você saberá que cada um dos pontos em que fará o loop pertencerá ao polígono atual. Aqui está um exemplo de código em python usando OGR que percorrerá os polígonos e os pontos de filtro de acordo. O código C / C ++ será muito parecido com isso, e eu imagino que você obterá um aumento significativo de velocidade em comparação com python. Você precisará adicionar algumas linhas de código para atualizar o CSV à medida que avança:

from osgeo import ogr
from osgeo.gdalconst import *

inPolyDS = ogr.Open("winnipeg.shp", GA_ReadOnly)
inPolyLayer = inPolyDS.GetLayer(0)
inPointDS = ogr.Open("busstops.vrt", GA_ReadOnly)   
inPointLayer = inPointDS.GetLayerByName("busstops")

inPolyFeat = inPolyLayer.GetNextFeature()
while inPolyFeat is not None:
  inPtFeat = inPointLayer.GetNextFeature()
  while inPtFeat is not None:
    ptGeom = inPtFeat.GetGeometryRef()
    # Do work here...

    inPtFeat = inPointLayer.GetNextFeature()

  inPolyFeat = inPolyLayer.GetNextFeature()

Arquivo VRT (busstops.vrt):

<OGRVRTDataSource>
  <OGRVRTLayer name="busstops">
    <SrcDataSource>busstops.csv</SrcDataSource>
    <GeometryType>wkbPoint</GeometryType>
    <LayerSRS>WGS84</LayerSRS>
    <GeometryField encoding="PointFromColumns" x="X" y="Y" reportSrcColumn="FALSE" />
  </OGRVRTLayer>
</OGRVRTDataSource>

Arquivo CSV (busstops.csv):

FID,X,Y,stop_name
1,-97.1394781371062,49.8712241633646,Southbound Osborne at Mulvey

Arquivo CSVT (busstops.csvt, o OGR precisa dele para identificar os tipos de coluna, caso contrário, ele não executará o filtro espacial):

Integer,Real,Real,String
Sasa Ivetic
fonte
2
Isso não passa por 1 bilhão de pontos 5000 vezes (uma vez para cada polígono)?
Underdark
Um índice espacial é um absoluto mosto . Eu mencionei Rtree antes, e vou mencionar novamente!
Mike T