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

9

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?

derNincompoop
la source

Réponses:

6

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:

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

mgri
la source
Je comprends que l'indice spatial fonctionnera mieux lorsqu'il sera en mesure de réduire le moins possible les géométries d'intérêt. C'est pourquoi j'ai comparé le nombre de géométries d'intérêt lors de l'utilisation de la méthode naïve (12936) au nombre de géométries lors de l'utilisation de l'indice spatial (1816). La intersects()méthode prend plus de temps lorsque l'index spatial est utilisé (voir la comparaison temporelle ci-dessus), c'est pourquoi je ne suis pas sûr d'utiliser incorrectement l'index spatial. En lisant la documentation et les articles liés, je pense que je le suis, mais j'espérais que quelqu'un pourrait faire remarquer si je ne l'étais pas.
derNincompoop
Votre déclaration finale était: "Je ne peux que penser que j'utilise mal l'index spatial de Rtree.", Donc je pensais que vous étiez confus quant à l'utilisation de l'index spatial et j'ai souligné sa signification dans ma réponse sujet). Vous essayez d'effectuer une sorte d'analyse statistique, mais le nombre de géométries et de tentatives impliquées ne devrait pas être suffisant pour une meilleure compréhension du problème. Ce comportement peut dépendre du nombre de géométries impliquées (un très petit nombre pour apprécier la puissance de l'indice spatial) ou de votre machine.
mgri
4

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

entrez la description de l'image ici

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)

entrez la description de l'image ici

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)

entrez la description de l'image ici

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

entrez la description de l'image ici

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 ...
gène
la source
Lorsque j'utilise la méthode naïve (où j'effectue un test d'intersection entre chaque combinaison Polygon et LineString), je finis par effectuer 12936 de ces tests. Lorsque j'utilise l'index spatial, je n'ai qu'à effectuer 1816 tests. Je pense que cela signifie que l'indice spatial fournit de la valeur dans ce cas d'utilisation. Cependant, lorsque je chronomètre le code, effectuer les tests 1816 prend COMME LONG comme effectuer les tests 12936. Le code avec l'index spatial ne devrait-il pas être plus rapide puisqu'il y a plus de 11 000 tests de moins en cours?
derNincompoop
J'ai donc examiné cela et constaté que les tests ~ 11000 effectués par le code naïf seulement prenaient moins d'une seconde à effectuer, tandis que les tests 1816 effectués par les deux ensembles de code prenaient 112 secondes. Maintenant, je comprends ce que vous entendez par `` L'index spatial n'est pas une baguette magique '' - même s'il réduisait le nombre de tests nécessaires, ceux qui étaient nécessaires étaient ceux qui avaient le plus contribué à l'époque.
derNincompoop
2

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).

derNincompoop
la source
1
Vous ne pouvez pas prendre en compte la façon dont l'indice spatiale influence la vitesse de l'autre intersection , car il ne fait partie de la complexité des géométries impliquées. Pour votre cas, si les géométries "A" et "B" se coupent en 0,05 seconde, elles se coupent toujours en 0,05 seconde même si vous avez utilisé précédemment un indice spatial (il s'agit évidemment d'un énoncé théorique, car je pense que vous savez que le traitement de quoi que ce soit dans un processeur est lié à de nombreux autres facteurs!).
mgri