Jak skutecznie uzyskać dostęp do funkcji zwróconych przez QgsSpatialIndex?

PyQGIS Cookbook wyjaśnia jak skonfigurować indeks przestrzenny, ale to tylko wyjaśnia połowę jego użytkowania:

utwórz indeks przestrzenny - poniższy kod tworzy pusty indeks

index = QgsSpatialIndex()

dodaj funkcje do indeksu - indeks pobiera obiekt QgsFeature i dodaje go do wewnętrznej struktury danych. Możesz utworzyć obiekt ręcznie lub użyć jednego z poprzedniego wywołania nextFeature () dostawcy

index.insertFeature(feat)

po wypełnieniu indeksu przestrzennego pewnymi wartościami można wykonać zapytania

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

Jaki jest najbardziej efektywny krok, aby uzyskać rzeczywiste funkcje należące do zwróconych identyfikatorów funkcji?

    # 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>

W poście na ten temat Nathan Woodrow podaje następujący kod:

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)

Spowoduje to utworzenie słownika umożliwiającego szybkie wyszukiwanie funkcji QgsFeature za pomocą jej FID.

Przekonałem się, że w przypadku bardzo dużych warstw nie jest to szczególnie praktyczne, ponieważ wymaga dużo pamięci. Jednak użycie alternatywne (losowy dostęp do pożądanej funkcji) layer.getFeatures(QgsFeatureRequest().setFilterFid(fid))wydaje się stosunkowo wolne. Nie jestem pewien, dlaczego tak jest, ponieważ równoważne wywołanie za pomocą powiązań SWIG OGR layer.GetFeature(fid)wydaje się znacznie szybsze.

Dla porównania spójrz na Efektywniejsze łączenie przestrzenne w Pythonie bez QGIS, ArcGIS, PostGIS itp . Przedstawione rozwiązanie wykorzystuje moduły Python Fiona , Shapely i rtree (Spatial Index).

Z PyQGIS i tym samym przykładem dwie warstwy pointi polygon:

1) Bez indeksu przestrzennego:

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) Z indeksem przestrzennym PyQGIS R-Tree :

# 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

Co oznaczają te wskaźniki?

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

Te same wnioski, co w przypadku bardziej wydajnego łączenia przestrzennego w Pythonie bez QGIS, ArcGIS, PostGIS itp . :

• Bez i indeksowanie, musisz iterować przez wszystkie geometrie (wielokąty i punkty).
• Z ograniczającym indeksem przestrzennym (QgsSpatialIndex ()), iterujesz tylko przez geometrie, które mają szansę przeciąć się z twoją bieżącą geometrią („filtr”, który może zaoszczędzić znaczną ilość obliczeń i czasu ...).
• W PyQGIS można również używać innych modułów indeksów przestrzennych Pythona ( rtree , Pyrtree lub Quadtree ), jak w przypadku korzystania z indeksu przestrzennego QGIS w celu przyspieszenia kodu (za pomocą QgsSpatialIndex () i rtree )
• ale Indeks przestrzenny nie jest magiczną różdżką. Gdy trzeba pobrać bardzo dużą część zestawu danych, Indeks przestrzenny nie może dać żadnej korzyści w zakresie prędkości.

Inny przykład w GIS se: Jak znaleźć najbliższą linię do punktu w QGIS? [duplikować]

Najwyraźniej jedyną metodą uzyskania dobrej wydajności jest unikanie lub łączenie wywołań z layer.getFeatures (), nawet jeśli filtr jest tak prosty jak fid.

Oto pułapka: wywołanie getFeatures jest drogie. Jeśli wywołasz go na warstwie wektorowej, QGIS będzie musiał skonfigurować nowe połączenie ze składnicą danych (dostawcą warstw), utworzyć zapytanie w celu zwrócenia danych i przeanalizować każdy wynik w miarę jego zwracania od dostawcy. Może to być powolne, szczególnie jeśli pracujesz z jakimś rodzajem zdalnej warstwy, takim jak tabela PostGIS przez połączenie VPN.

źródło: http://nyalldawson.net/2016/10/speeding-up-your-pyqgis-scripts/