O índice espacial RTree não resulta em computação de interseção mais rápida

9

Eu tenho algum código que estou usando para determinar quais polígonos esculturais / MultiPolygons se cruzam com um número de linhas de formas bem torneadas. Através das respostas a esta pergunta, o código passou disso:

import fiona
from shapely.geometry import LineString, Polygon, MultiPolygon, shape

# Open each layer
poly_layer = fiona.open('polygon_layer.shp')
line_layer = fiona.open('line_layer.shp')

# Convert to lists of shapely geometries
the_lines = [shape(line['geometry']) for line in line_layer]
the_polygons = [(poly['properties']['GEOID'], shape(poly['geometry'])) for poly in poly_layer]

# Check for Polygons/MultiPolygons that the LineString intersects with
covered_polygons = {}
for poly_id, poly in the_polygons:
    for line in the_lines:
        if poly.intersects(line):
            covered_polygons[poly_id] = covered_polygons.get(poly_id, 0) + 1

onde todas as interseções possíveis são verificadas, para isso:

import fiona
from shapely.geometry import LineString, Polygon, MultiPolygon, shape
import rtree

# Open each layer
poly_layer = fiona.open('polygon_layer.shp')
line_layer = fiona.open('line_layer.shp')

# Convert to lists of shapely geometries
the_lines = [shape(line['geometry']) for line in line_layer]
the_polygons = [(poly['properties']['GEOID'], shape(poly['geometry'])) for poly in poly_layer]

# Create spatial index
spatial_index = rtree.index.Index()
for idx, poly_tuple in enumerate(the_polygons):
    _, poly = poly_tuple
    spatial_index.insert(idx, poly.bounds)

# Check for Polygons/MultiPolygons that the LineString intersects with
covered_polygons = {}
for line in the_lines:
    for idx in list(spatial_index.intersection(line.bounds)):
        if the_polygons[idx][1].intersects(line):
            covered_polygons[idx] = covered_polygons.get(idx, 0) + 1

onde o índice espacial é usado para reduzir o número de verificações de interseção.

Com os shapefiles que tenho (aproximadamente 4000 polígonos e 4 linhas), o código original executa 12936 .intersection()verificações e leva cerca de 114 segundos para executar. O segundo pedaço de código que usa o índice espacial executa apenas 1816 .intersection()verificações, mas também leva aproximadamente 114 segundos para ser executado.

O código para criar o índice espacial leva apenas 1-2 segundos para ser executado. Portanto, as verificações de 1816 no segundo trecho de código estão demorando quase a mesma quantidade de tempo que as 12936 verificam no código original (desde o carregamento do shapefiles e a conversão para geometrias Shapely são iguais em ambas as partes do código).

Não vejo nenhuma razão para que o índice espacial .intersects()demore mais tempo a verificação; portanto, não sei por que isso está acontecendo.

Só consigo pensar que estou usando o índice espacial RTree incorretamente. Pensamentos?

derNincompoop
fonte

Respostas:

6

Minha resposta é essencialmente baseada em outra resposta da @gene aqui:

Junção espacial mais eficiente em Python sem QGIS, ArcGIS, PostGIS, etc

Ele propôs a mesma solução usando dois métodos diferentes, com ou sem um índice espacial.

Ele (corretamente) afirmou:

Qual é a diferença ?

  • Sem o índice, você deve percorrer todas as geometrias (polígonos e pontos).
  • Com um índice espacial delimitador ( Spatial Index RTree ), você itera apenas pelas geometrias que têm a chance de se cruzar com a geometria atual ('filtro' que pode economizar uma quantidade considerável de cálculos e tempo ...)
  • mas um índice espacial não é uma varinha mágica. Quando uma parte muito grande do conjunto de dados precisa ser recuperada, um Índice Espacial não pode oferecer nenhum benefício de velocidade.

Essas frases são auto-explicativas, mas eu preferi citar @gene em vez de propor as mesmas conclusões que as minhas (então, todos os créditos vão para o seu brilhante trabalho!).

Para uma melhor compreensão do índice espacial da Rtree, você pode recuperar algumas informações úteis seguindo estes links:

Outra ótima introdução sobre o uso de índices espaciais pode ser este artigo de @Nathan Woodrow .

mgri
fonte
Entendo que o índice espacial funcionará melhor quando for capaz de reduzir as geometrias de interesse ao mínimo possível. Por isso, comparei o número de geometrias de interesse ao usar o método ingênuo (12936) com o número de geometrias ao usar o índice espacial (1816). O intersects()método leva mais tempo quando o índice espacial está sendo usado (veja a comparação de tempos acima), e é por isso que não tenho certeza se estou usando o índice espacial incorretamente. Ao ler a documentação e os posts vinculados, acho que sou, mas esperava que alguém pudesse apontar se não o fizesse.
DerNincompoop
Sua declaração final foi: "Só consigo pensar que estou usando o índice espacial de Rtree incorretamente." tópico). Você está tentando executar algum tipo de análise estatística, mas o número de geometrias e tentativas envolvidas não deve ser suficiente para uma melhor compreensão do problema. Esse comportamento pode depender do número de geometrias envolvidas (um número muito pequeno para apreciar a potência do índice espacial) ou da sua máquina.
MGRI
4

Apenas para adicionar à resposta mgri.

É importante entender o que é um índice espacial ( como implementar corretamente uma caixa delimitadora para Shapely & Fiona? ). Com o meu exemplo em Como determinar com eficiência quais milhares de polígonos se cruzam com uma cadeia de linhas

insira a descrição da imagem aqui

Você pode criar um índice espacial com os polígonos

idx = index.Index()
for feat in poly_layer:
    geom = shape(feat['geometry'])
    id = int(feat['id'])
    idx.insert(id, geom.bounds,feat)

Limites do índice espacial (limites do polígono em verde)

insira a descrição da imagem aqui

Ou com as LineStrings

  idx = index.Index()
  for feat in line_layer:
      geom = shape(feat['geometry'])
      id = int(feat['id'])
      idx.insert(id, geom.bounds,feat)

Limites do índice espacial (LineString vinculado em vermelho)

insira a descrição da imagem aqui

Agora, você itera apenas pelas geometrias que têm a chance de se cruzarem com a geometria atual (em amarelo)

insira a descrição da imagem aqui

Eu uso aqui o índice espacial LineStrings (os resultados são os mesmos, mas com o seu exemplo de 4000 polígonos e 4 linhas ...).

for feat1 in poly_layer:
    geom1 = shape(feat1['geometry'])
    for id in idx.intersection(geom1.bounds):
        feat2 = line_layer[id]
        geom2 = shape(feat2['geometry'])
        if geom2.intersects(geom1):
            print 'Polygon {} intersects line {}'.format(feat1['id'], feat2['id'])

  Polygon 2 intersects line 0
  Polygon 3 intersects line 0
  Polygon 6 intersects line 0
  Polygon 9 intersects line 0

Você também pode usar um gerador ( exemplo.py )

def build_ind():
     with fiona.open('polygon_layer.shp') as shapes:
         for s in shapes:
             geom = shape(s['geometry'])
             id = int(s['id'])
             yield (id, geom.bounds, s)

 p= index.Property()
 tree = index.Index(build_ind(), properties=p)
 # first line of line_layer
 first = shape(line_layer.next()['geometry'])
 # intersection of first with polygons
 tuple(tree.intersection(first.bounds))
 (6, 2, 3, 9)

Você pode examinar o script GeoPandas sjoin.py para entender o uso do Rtree .

Existem muitas soluções, mas não esqueça que

  • O índice espacial não é uma varinha mágica ...
gene
fonte
Quando uso o método ingênuo (onde realizo um teste de interseção entre todas as combinações Polygon e LineString), acabo executando 12936 desses testes. Quando uso o índice espacial, só preciso realizar 1816 testes. Acredito que isso significa que o índice espacial fornece valor nesse caso de uso. No entanto, quando cronometro o código, a execução dos testes de 1816 leva tanto tempo quanto a realização dos testes de 12936. O código com o índice espacial não deveria ser mais rápido, pois há mais de 11000 testes a serem realizados?
DerNincompoop 8/02
Então, analisei isso e descobri que os ~ 11000 testes realizados apenas pelo código ingênuo levam menos de 1 segundo para serem executados, enquanto os testes de 1816 realizados pelos dois conjuntos de códigos levam 112 segundos para serem executados. Agora entendo o que você quer dizer com 'Índice espacial não é uma varinha mágica' - mesmo que diminuísse o número de testes necessários, os necessários foram os que mais contribuíram para a época.
DerNincompoop
2

Editar: para esclarecer esta resposta, eu acreditava incorretamente que todos os testes de interseção levavam aproximadamente a mesma quantidade de tempo. Este não é o caso. O motivo pelo qual não obtive a velocidade que eu esperava usando um índice espacial é que a seleção dos testes de interseção é a que demorou mais tempo para fazer.

Como gene e mgri já disseram, um índice espacial não é uma varinha mágica. Embora o índice espacial reduzisse o número de testes de interseção que precisavam ser executados de 12936 a apenas 1816, os testes de 1816 são o teste que levou a maior parte do tempo para calcular em primeiro lugar.

O índice espacial está sendo usado corretamente, mas a suposição de que cada teste de interseção leva aproximadamente a mesma quantidade de tempo é o que está incorreto. O tempo requerido pelo teste de interseção pode variar bastante (0,05 segundos versus 0,000007 segundos).

derNincompoop
fonte
1
Você não pode levar em consideração como o índice espacial influencia a velocidade da interseção posterior, porque pertence apenas à complexidade das geometrias envolvidas. No seu caso, se as geometrias "A" e "B" se cruzarem em 0,05 segundos, elas ainda se cruzarão em 0,05 segundos, mesmo que você tenha usado anteriormente um índice espacial (essa é obviamente uma afirmação teórica, porque acho que você sabe que o processamento de qualquer coisa dentro de um processador está relacionada a muitos outros fatores!).
MGRI