L'indice spatial RTree n'entraîne pas un calcul d'intersection plus rapide

J'ai du code que j'utilise pour déterminer quels polygones / multi-polygones Shapely se croisent avec un certain nombre de chaînes de lignes Shapely. Grâce aux réponses à cette question, le code est passé de ceci:

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

où chaque intersection possible est vérifiée, à ceci:

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

où l'indice spatial est utilisé pour réduire le nombre de vérifications d'intersection.

Avec les fichiers de formes que j'ai (environ 4000 polygones et 4 lignes), le code d'origine effectue 12936 .intersection()vérifications et prend environ 114 secondes pour s'exécuter. Le deuxième morceau de code qui utilise l'index spatial effectue seulement 1816 .intersection()vérifications, mais il faut également environ 114 secondes pour s'exécuter.

Le code pour construire l'index spatial ne prend que 1 à 2 secondes à exécuter, donc les vérifications 1816 dans le deuxième morceau de code prennent à peu près le même temps à effectuer que les vérifications 12936 dans le code d'origine (depuis le chargement de les fichiers de formes et la conversion en géométries Shapely sont les mêmes dans les deux morceaux de code).

Je ne vois aucune raison pour laquelle l'index spatial rendrait la .intersects()vérification plus longue, donc je ne comprends pas pourquoi cela se produit.

Je peux seulement penser que j'utilise incorrectement l'index spatial RTree. Pensées?

Ma réponse est essentiellement basée sur une autre réponse de @gene ici:

Jointure spatiale plus efficace en Python sans QGIS, ArcGIS, PostGIS, etc.

Il a proposé la même solution en utilisant deux méthodes différentes, avec ou sans indice spatial.

Il a (à juste titre) déclaré:

Quelle est la différence ?

• Sans l'index, vous devez parcourir toutes les géométries (polygones et points).
• Avec un index spatial englobant ( Spatial Index RTree ), vous itérez uniquement à travers les géométries qui ont une chance de se croiser avec votre géométrie actuelle (`` filtre '' qui peut économiser une quantité considérable de calculs et de temps ...)
• mais un indice spatial n'est pas une baguette magique. Lorsqu'une très grande partie de l'ensemble de données doit être récupérée, un index spatial ne peut pas bénéficier de la vitesse.

Ces phrases sont explicites, mais j'ai préféré citer @gene au lieu de proposer les mêmes conclusions que les miennes (donc, tous les crédits vont à son brillant travail!).

Pour une meilleure compréhension de l'index spatial Rtree, vous pouvez récupérer quelques informations utiles en suivant ces liens:

• Rtree: Indexation spatiale pour Python ;
• Utilisation de l'indexation spatiale Rtree avec OGR .

Une autre excellente introduction à l'utilisation des index spatiaux peut être cet article de @Nathan Woodrow .

Juste pour ajouter à la réponse mgri.

Il est important de comprendre ce qu'est un index spatial ( Comment implémenter correctement un cadre de délimitation pour Shapely & Fiona? ). Avec mon exemple dans Comment déterminer efficacement lequel des milliers de polygones se croisent avec une chaîne de lignes

Vous pouvez créer un index spatial avec les polygones

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

Limites de l'index spatial (limites des polygones en vert)

Ou avec les LineStrings

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

Limites de l'index spatial (LineString lié en rouge)

Maintenant, vous parcourez uniquement les géométries qui ont une chance de se croiser avec votre géométrie actuelle (en jaune)

J'utilise ici l'index spatial LineStrings (les résultats sont les mêmes mais avec votre exemple de 4000 polygones et 4 lignes ....).

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

Vous pouvez également utiliser un générateur ( exemple.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)

Vous pouvez examiner le script GeoPandas sjoin.py pour comprendre l'utilisation de Rtree .

Il existe de nombreuses solutions mais n'oubliez pas que

• L'index spatial n'est pas une baguette magique ...

Edit: Pour clarifier cette réponse, j'ai cru à tort que tous les tests d'intersection ont pris environ le même temps. Ce n'est pas le cas. La raison pour laquelle je n'ai pas obtenu la vitesse à laquelle j'attendais en utilisant un indice spatial est que la sélection des tests d'intersection est celle qui a pris le plus de temps à faire en premier lieu.

Comme l'ont déjà dit gène et mgri, un indice spatial n'est pas une baguette magique. Bien que l'indice spatial ait réduit le nombre de tests d'intersection qui devaient être effectués de 12936 à seulement 1816, les tests de 1816 sont ceux qui ont pris la grande majorité du temps à calculer en premier lieu.

L'indice spatial est utilisé correctement, mais l'hypothèse selon laquelle chaque test d'intersection prend à peu près le même temps est ce qui est incorrect. Le temps requis par le test d'intersection peut varier considérablement (0,05 seconde contre 0,000007 seconde).