Comprendre l'utilisation des index spatiaux avec RTree?

J'ai du mal à comprendre l'utilisation des index spatiaux avec RTree.

Exemple: J'ai 300 points tamponnés et j'ai besoin de connaître la zone d'intersection de chaque tampon avec un fichier de formes polygonal. Le fichier de formes des polygones contient> 20 000 polygones. Il a été suggéré d'utiliser des indices spatiaux pour accélérer le processus.

SO ... Si je crée un index spatial pour mon fichier de formes de polygones, sera-t-il "attaché" au fichier d'une manière ou d'une autre, ou l'index sera-t-il autonome? Autrement dit, après l'avoir créé, puis-je simplement exécuter ma fonction d'intersection sur le fichier polygonal et obtenir des résultats plus rapides? L'intersection "verra-t-elle" qu'il existe des indices spatiaux et saura-t-elle quoi faire? Ou, dois-je l'exécuter sur l'index, puis relier ces résultats à mon fichier polygonal d'origine via des FID ou quelque chose du genre?

La documentation RTree ne m'aide pas beaucoup (probablement parce que j'apprends juste la programmation). Ils montrent comment créer un index en lisant des points créés manuellement, puis en l'interrogeant sur d'autres points créés manuellement, ce qui renvoie les identifiants contenus dans la fenêtre. Logique. Mais, ils n'expliquent pas comment cela se rapporterait à un fichier d'origine dont proviendrait l'index.

Je pense que ça doit aller quelque chose comme ça:

1. Tirez des bboxes pour chaque entité polygonale de mon fichier de formes de polygone et placez-les dans un index spatial, en leur donnant un identifiant identique à leur identifiant dans le fichier de formes.
2. Recherchez cet index pour obtenir les identifiants qui se croisent.
3. Ensuite, réexécutez mon intersection sur uniquement les entités de mon fichier de formes d'origine qui ont été identifiées en interrogeant mon index (je ne sais pas comment je ferais cette dernière partie).

Ai-je la bonne idée? Suis-je en train de manquer quelque chose?

En ce moment, j'essaie de faire fonctionner ce code sur un fichier de formes ponctuelles qui ne contient qu'une seule entité ponctuelle et un fichier de formes polygonales qui contient> 20 000 entités surfaciques.

J'importe les fichiers de formes à l'aide de Fiona, j'ajoute l'index spatial à l'aide de RTree et j'essaie de faire l'intersection à l'aide de Shapely.

Mon code de test ressemble à ceci:

#point shapefile representing location of desired focal statistic
traps = fiona.open('single_pt_speed_test.shp', 'r') 

#polygon shapefile representing land cover of interest 
gl = MultiPolygon([shape(pol['geometry']) for pol in fiona.open('class3_aa.shp', 'r')]) 

#search area
areaKM2 = 20

#create empty spatial index
idx = index.Index()

#set initial search radius for buffer
areaM2 = areaKM2 * 1000000
r = (math.sqrt(areaM2/math.pi))

#create spatial index from gl
for i, shape in enumerate(gl):
    idx.insert(i, shape.bounds)

#query index for ids that intersect with buffer (will eventually have multiple points)
for point in traps:
        pt_buffer = shape(point['geometry']).buffer(r)
        intersect_ids = pt_buffer.intersection(idx)

Mais je reçois toujours TypeError: l'objet 'Polygon' n'est pas appelable

Voilà l'essentiel. L'arbre R vous permet de faire un premier passage très rapide et vous donne un ensemble de résultats qui auront des "faux positifs" (les boîtes englobantes peuvent se croiser lorsque les géométries ne le sont pas précisément). Ensuite, vous passez en revue l'ensemble des candidats (en les récupérant du fichier de formes par leur index) et effectuez un test d'intersection mathématiquement précis en utilisant, par exemple, Shapely. Il s'agit de la même stratégie que celle utilisée dans les bases de données spatiales comme PostGIS.

Vous l'avez presque compris, mais vous avez fait une petite erreur. Vous devez utiliser la intersectionméthode sur l'index spatial, plutôt que de passer l'index à la intersectionméthode sur le point tamponné. Une fois que vous avez trouvé une liste d'entités où les zones de délimitation se chevauchent, vous devez vérifier si votre point tamponné intersecte réellement les géométries.

import fiona
from shapely.geometry import mapping
import rtree
import math

areaM2 = areaKM2 * 1000000
r = (math.sqrt(areaM2/math.pi))

# open both layers
with fiona.open('single_pt_speed_test.shp', 'r') as layer_pnt:
    with fiona.open('class3_aa.shp', 'r') as layer_land:

        # create an empty spatial index object
        index = rtree.index.Index()

        # populate the spatial index
        for fid, feature in layer_land.items():
            geometry = shape(feature['geometry'])
            idx.insert(fid, geometry.bounds)

        for feature in layer_pnt:
            # buffer the point
            geometry = shape(feature['geometry'])
            geometry_buffered = geometry.buffer(r)

            # get list of fids where bounding boxes intersect
            fids = [int(i) for i in index.intersection(geometry_buffered.bounds)]

            # access the features that those fids reference
            for fid in fids:
                feature_land = layer_land[fid]
                geometry_land = shape(feature_land['geometry'])

                # check the geometries intersect, not just their bboxs
                if geometry.intersects(geometry_land):
                    print('Found an intersection!')  # do something useful here

Si vous êtes intéressé à trouver des points qui se trouvent à une distance minimale de votre classe de terrain, vous pouvez utiliser la distanceméthode à la place (remplacer la section appropriée de la précédente).

for feature in layer_pnt:
    geometry = shape(feature['geometry'])

    # expand bounds by r in all directions
    bounds = [a+b*r for a,b in zip(geometry.bounds, [-1, -1, 1, 1])]

    # get list of fids where bounding boxes intersect
    fids = [int(i) for i in index.intersection(geometry_buffered.bounds)]

    for fid in fids:
        feature_land = layer_land[fid]
        geometry_land = shape(feature_land['geometry'])

        # check the geometries are within r metres
        if geometry.distance(geometry_land) <= r:
            print('Found a match!')

Si la construction de votre index spatial prend beaucoup de temps et que vous allez le faire plusieurs fois, vous devriez envisager de sérialiser l'index dans un fichier. La documentation décrit comment procéder: http://toblerity.org/rtree/tutorial.html#serializing-your-index-to-a-file

Vous pouvez également regarder le chargement en masse des boîtes englobantes dans le rtree à l'aide d'un générateur, comme ceci:

def gen(collection):
    for fid, feature in collection.items():
        geometry = shape(feature['geometry'])
        yield((fid, geometry.bounds, None))
index = rtree.index.Index(gen(layer_land))

Oui, c'est l'idée. Voici un extrait de ce tutoriel sur l'utilisation d'un index spatial r-tree en Python , en utilisant shapely, Fiona et geopandas:

Un arbre r représente les objets individuels et leurs boîtes englobantes (le «r» signifie «rectangle») comme le niveau le plus bas de l'index spatial. Il agrège ensuite les objets voisins et les représente avec leur boîte englobante agrégée au niveau supérieur suivant de l'index. À des niveaux encore plus élevés, l'arborescence r agrège les boîtes englobantes et les représente par leur boîte englobante, de manière itérative, jusqu'à ce que tout soit imbriqué dans une boîte englobante de niveau supérieur. Pour effectuer une recherche, l'arborescence r prend une zone de requête et, en commençant au niveau supérieur, voit quelles (le cas échéant) zones de délimitation l'intersectent. Il développe ensuite chaque zone de délimitation entrecroisée et voit lesquelles des zones de délimitation enfant à l'intérieur se croisent dans la zone de requête. Cela se poursuit récursivement jusqu'à ce que toutes les cases qui se croisent soient recherchées jusqu'au niveau le plus bas et renvoie les objets correspondants du niveau le plus bas.