RTree índice espacial no es el resultado en el más rápido de la intersección de la computación

Question

Tengo algo de código que estoy utilizando para determinar que bien formada Polígono/Multipolígonos se cruzan con un número de frijoles LineStrings. A través de las respuestas a esta pregunta el código ha pasado de esto:

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

donde cada intersección se comprueba, a esta:

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

donde el índice espacial se utiliza para reducir el número de intersección de cheques.

Con los shapefiles que tengo (aproximadamente 4000 polígonos, y 4 líneas), el código original se realiza 12936 .intersection() cheques y toma alrededor de 114 segundos. El segundo trozo de código que utiliza el índice espacial realiza sólo 1816 .intersection() los controles, pero también toma aproximadamente 114 segundos.

El código para generar el índice espacial sólo se tarda de 1-2 segundos para correr, por lo que el 1816 comprueba en la segunda pieza de código están teniendo casi la misma cantidad de tiempo para llevar a cabo como la 12936 cheques en el código original (ya que la carga de archivos y convertir a bien formada geometrías es el mismo en ambas piezas de código).

No puedo ver ninguna razón por la que el índice espacial haría el .intersects() de verificación tomar más tiempo, así que estoy en una pérdida de por qué está sucediendo esto.

Solo puedo pensar que estoy utilizando el RTree índice espacial incorrectamente. Los pensamientos?

Answer 1

Mi respuesta se basa esencialmente en la otra respuesta por @gen aquí:

Más Eficiente de unión Espacial en Python sin QGIS, ArcGIS, PostGIS, etc

Él propone la misma solución utilizando dos diferentes métodos, con o sin un índice espacial.

Él (correctamente) declaró:

¿Cuál es la diferencia ?

• Sin el índice, se debe iterar a través de todas las geometrías (polígonos y puntos).
• Con una delimitación espacial del índice (Índice Espacial RTree), se itera sólo a través de las geometrías que tienen una oportunidad para que se cruzan con su geometría actual ('filtro' que puede ahorrar una considerable cantidad de cálculos y el tiempo...)
• pero un Índice Espacial no es una varita mágica. Cuando una gran parte del conjunto de datos tiene que ser recuperada, de un Índice Espacial no puede dar ninguna ventaja de velocidad.

Estas frases son auto-explicativas, pero he preferido citar @gen en lugar de proponer a las mismas conclusiones como la mía (por lo tanto, todos los créditos van a su brillante trabajo!).

Para una mejor comprensión de la Rtree índice espacial, usted puede recuperar algo de información útil de los siguientes enlaces:

• Rtree: los índices Espaciales para Python;
• El uso de Rtree espacial de la indización con OGR.

Otro gran introducción sobre el uso de índices espaciales puede ser este artículo de @Nathan Woodrow.

Answer 2

Sólo para agregar a mgri respuesta.

Es importante entender lo que es un índice espacial (Cómo Implementar Correctamente un Cuadro Delimitador para bien formada & Fiona?). Con mi ejemplo, en Cómo de manera eficiente determinar cuál de miles de polígonos se cruzan con un linestring

Usted puede crear un índice espacial con los polígonos

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

Los límites del índice espacial (polígono de límites en verde)

O con el LineStrings

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

Los límites del índice espacial (LineString encuadernado en rojo)

Ahora, se itera sólo a través de las geometrías que tienen una oportunidad para intersecar con su geometría actual (en amarillo)

Utilizo aquí el LineStrings índice espacial (los resultados son los mismos, pero con su ejemplo de 4000 polígonos y 4 líneas ....).

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

También puede utilizar un generador (example.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)

Usted puede examinar el GeoPandas script sjoin.py para entender el uso de Rtree.

Hay muchas soluciones, pero no olvides que

• Índice espacial no es una varita mágica...

Answer 3

Edit: Para aclarar esta respuesta, me creían equivocadamente que todos la intersección de las pruebas se realizó en aproximadamente la misma cantidad de tiempo. Este no es el caso. La razón por la que no obtenga la velocidad hasta que yo esperaba desde el uso de un índice espacial es que la selección de intersección de las pruebas son los que más tardó en hacer en el primer lugar.

Como gen y mgri ya han dicho, un índice espacial no es una varita mágica. Aunque el índice espacial redujo el número de intersección de las pruebas que necesita para realizarse a partir de 12936 sólo de 1816, el de 1816 pruebas son las pruebas que tuvo la inmensa mayoría de tiempo para calcular en primer lugar.

El índice espacial está siendo utilizado correctamente, pero la suposición de que cada intersección de la prueba toma aproximadamente la misma cantidad de tiempo es lo que es incorrecto. El tiempo requerido por la intersección de la prueba puede variar enormemente (0.05 segundos frente a 0.000007 segundos).