Wie kann ich effizient auf die von QgsSpatialIndex zurückgegebenen Funktionen zugreifen?

Das PyQGIS-Kochbuch erklärt, wie der räumliche Index eingerichtet wird, erklärt jedoch nur die Hälfte seiner Verwendung:

räumlichen Index erstellen - Der folgende Code erstellt einen leeren Index

index = QgsSpatialIndex()

Features zum Index hinzufügen - index nimmt das QgsFeature-Objekt und fügt es der internen Datenstruktur hinzu. Sie können das Objekt manuell erstellen oder eines aus dem vorherigen Aufruf von nextFeature () des Anbieters verwenden.

index.insertFeature(feat)

Sobald der räumliche Index mit einigen Werten gefüllt ist, können Sie einige Abfragen durchführen

# returns array of feature IDs of five nearest features
nearest = index.nearestNeighbor(QgsPoint(25.4, 12.7), 5)

Was ist der effizienteste Schritt, um die tatsächlichen Features zu erhalten, die zu den zurückgegebenen Feature-IDs gehören?

    # assume a list of feature ids returned from index and a QgsVectorLayer 'lyr'
    fids = [1, 2, 4]
    request = QgsFeatureRequest()
    request.setFilterFids(fids)

    features = lyr.getFeatures(request)
    # can now iterate and do fun stuff:
    for feature in features:
        print feature.id(), feature

    1 <qgis._core.QgsFeature object at 0x000000000E987510>
    2 <qgis._core.QgsFeature object at 0x000000000E987400>
    4 <qgis._core.QgsFeature object at 0x000000000E987510>

In einem Blog-Beitrag zu diesem Thema stellt Nathan Woodrow den folgenden Code bereit:

layer = qgis.utils.iface.activeLayer()

# Select all features along with their attributes
allAttrs = layer.pendingAllAttributesList()
layer.select(allAttrs)
# Get all the features to start
allfeatures = {feature.id(): feature for (feature) in layer}

def noindex():
    for feature in allfeatures.values():
        for f in allfeatures.values():
            touches = f.geometry().touches(feature.geometry())
            # It doesn't matter if we don't return anything it's just an example

def withindex():
    # Build the spatial index for faster lookup.
    index = QgsSpatialIndex()
    map(index.insertFeature, allfeatures.values())

    # Loop each feature in the layer again and get only the features that are going to touch.
    for feature in allfeatures.values():
        ids = index.intersects(feature.geometry().boundingBox())
        for id in ids:
            f = allfeatures[id]
            touches = f.geometry().touches(feature.geometry())
            # It doesn't matter if we don't return anything it's just an example

import timeit
print "With Index: %s seconds " % timeit.timeit(withindex,number=1)
print "Without Index: %s seconds " % timeit.timeit(noindex,number=1)

Dadurch wird ein Wörterbuch erstellt, mit dem Sie ein QgsFeature mithilfe seiner FID schnell nachschlagen können.

Ich habe festgestellt, dass dies für sehr große Schichten nicht besonders praktisch ist, da es viel Speicher benötigt. Die alternative Verwendung (zufälliger Zugriff auf das gewünschte Merkmal) layer.getFeatures(QgsFeatureRequest().setFilterFid(fid))scheint jedoch vergleichsweise sehr langsam zu sein. Ich bin mir nicht sicher, warum dies so ist, da der entsprechende Aufruf mit den SWIG OGR-Bindungen layer.GetFeature(fid)viel schneller zu sein scheint.

Für Vergleiche, Blick auf Effizientere Spatial Join in Python ohne QGIS, ArcGIS, PostGIS, etc . Die vorgestellte Lösung verwendet die Python-Module Fiona , Shapely und rtree (Spatial Index).

Mit PyQGIS und demselben Beispiel zwei Ebenen pointund polygon:

1) Ohne räumlichen Index:

polygons = [feature for feature in polygon.getFeatures()]
points = [feature for feature in point.getFeatures()]
for pt in points: 
    point = pt.geometry()
    for pl  in polygons:
        poly = pl.geometry()
        if poly.contains(point):
            print point.asPoint(), poly.asPolygon()
(184127,122472) [[(183372,123361), (184078,123130), (184516,122631),   (184516,122265), (183676,122144), (183067,122570), (183128,123105), (183372,123361)]]
(183457,122850) [[(183372,123361), (184078,123130), (184516,122631), (184516,122265), (183676,122144), (183067,122570), (183128,123105), (183372,123361)]]
(184723,124043) [[(184200,124737), (185368,124372), (185466,124055), (185515,123714), (184955,123580), (184675,123471), (184139,123787), (184200,124737)]]
(182179,124067) [[(182520,125175), (183348,124286), (182605,123714), (182252,123544), (181753,123799), (181740,124627), (182520,125175)]]

2) Mit dem räumlichen Index R-Tree PyQGIS:

# build the spatial index with all the polygons and not only a bounding box
index = QgsSpatialIndex()
for poly in polygons:
     index.insertFeature(poly)

# intersections with the index 
# indices of the index for the intersections
for pt in points:
    point = pt.geometry()
    for id in index.intersects(point.boundingBox()):
    print id
0
0
1
2

Was bedeuten diese Indizes?

for i, pt in enumerate(points):
     point = pt.geometry()
     for id in index.intersects(point.boundingBox()):
        print "Point ", i, points[i].geometry().asPoint(), "is in Polygon ", id, polygons[id].geometry().asPolygon()
Point  1 (184127,122472) is in Polygon  0 [[(182520,125175), (183348,124286), (182605,123714), (182252,123544), (181753,123799), (181740,124627), (182520,125175)]]
Point  2 (183457,122850) is in Polygon  0 [[(182520,125175), (183348,124286), (182605,123714), (182252,123544), (181753,123799), (181740,124627), (182520,125175)]]
Point  4 (184723,124043) is in Polygon  1 [[(182520,125175), (183348,124286), (182605,123714), (182252,123544), (181753,123799), (181740,124627), (182520,125175)]]
Point  6 (182179,124067) is in Polygon  2 [[(182520,125175), (183348,124286), (182605,123714), (182252,123544), (181753,123799), (181740,124627), (182520,125175)]]

Gleiche Schlussfolgerungen wie bei Effizienterer räumlicher Join in Python ohne QGIS, ArcGIS, PostGIS usw . :

• Ohne und Index müssen Sie alle Geometrien (Polygone und Punkte) durchlaufen.
• Mit einem begrenzten räumlichen Index (QgsSpatialIndex ()) iterieren Sie nur durch die Geometrien, die sich mit Ihrer aktuellen Geometrie überschneiden können ('Filter', der eine beträchtliche Menge an Berechnungen und Zeit sparen kann ...).
• Sie können mit PyQGIS auch andere Python-Module mit räumlichem Index ( rtree , Pyrtree oder Quadtree ) verwenden, wie unter Verwenden eines räumlichen QGIS-Index, um Ihren Code zu beschleunigen (mit QgsSpatialIndex () und rtree ).
• Ein räumlicher Index ist jedoch kein Zauberstab. Wenn ein sehr großer Teil des Datensatzes abgerufen werden muss, kann ein räumlicher Index keinen Geschwindigkeitsvorteil bieten.

Anderes Beispiel in GIS se: Wie finde ich die nächstgelegene Linie zu einem Punkt in QGIS? [Duplikat]

Anscheinend besteht die einzige Methode, um eine gute Leistung zu erzielen, darin, Aufrufe von layer.getFeatures () zu vermeiden oder zu bündeln, selbst wenn der Filter so einfach wie ein FID ist.

Hier ist die Falle: Das Aufrufen von getFeatures ist teuer. Wenn Sie es auf einer Vektorebene aufrufen, muss QGIS eine neue Verbindung zum Datenspeicher (dem Ebenenanbieter) herstellen, eine Abfrage erstellen, um Daten zurückzugeben, und jedes Ergebnis analysieren, sobald es vom Anbieter zurückgegeben wird. Dies kann langsam sein, insbesondere wenn Sie mit einer Art Remote-Schicht arbeiten, z. B. einer PostGIS-Tabelle über eine VPN-Verbindung.

Quelle: http://nyalldawson.net/2016/10/speeding-up-your-pyqgis-scripts/