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

9

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?

Unterdunkel
quelle

Antworten:

12
    # 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>
gsherman
quelle
Vielen Dank! Snorfalorpagus erwähnte, dass setFilterFids erheblich langsamer sein würde als die von ihm veröffentlichte Lösung. Bestätigen Sie das?
Underdark
Ich habe es nicht für große Ergebnismengen verwendet und kann es daher nicht bestätigen.
Gsherman
1
Ich bestätige und in meinem Fall ist rtree sogar noch schneller als QgsSpatialIndex () (für die Erstellung von Planar-Graphen aus sehr großen Polylinienschichten, Transposition des Moduls PlanarGraph mit Shapely in PyQGIS. Aber die Lösung mit Fiona, Shapely und rtree ist immer noch die am schnellsten)
Gen
1
Ich glaube, die Frage bezieht sich eher auf das Abrufen der tatsächlichen Features aus den zurückgegebenen Feature-IDs als auf die Geschwindigkeit verschiedener Indizierungsmethoden.
Gsherman
7

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.

Snorfalorpagus
quelle
1
Die Verwendung eines Wörterbuchs war sehr viel schneller als layer.getFeatures(QgsFeatureRequest().setFilterFid(fid)). Ich habe an einer Ebene mit 140.000 Funktionen gearbeitet, und die Gesamtzeit für die 140.000 Suchvorgänge ging von vielen Minuten auf Sekunden.
Håvard Tveite
5

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:

Geben Sie hier die Bildbeschreibung ein

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]

Gen
quelle
Vielen Dank für die zusätzliche Erklärung. Grundsätzlich verwendet Ihre Lösung eine Liste anstelle eines Diktats wie Snorfalorpagus. Es scheint also wirklich keine Funktion auf layer.getFeatures ([ids]) zu geben ...
underdark
Der Zweck dieser Erklärung ist rein geometrisch und es ist sehr einfach, eine layer.getFeatures ([ids]) - Funktion hinzuzufügen, wie in Effizientere räumliche Verknüpfung in Python ohne QGIS, ArcGIS, PostGIS usw.
Gen
0

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/

evod
quelle