Bagaimana cara efisien mengakses fitur yang dikembalikan oleh QgsSpatialIndex?

9

The PyQGIS Cookbook menjelaskan cara mengatur indeks spasial tetapi hanya menjelaskan setengah dari penggunaannya:

buat indeks spasial - kode berikut membuat indeks kosong

index = QgsSpatialIndex()

menambahkan fitur ke indeks - indeks mengambil objek QgsFeature dan menambahkannya ke struktur data internal. Anda dapat membuat objek secara manual atau menggunakannya dari panggilan sebelumnya ke penyedia berikutnya ()

index.insertFeature(feat)

setelah indeks spasial diisi dengan beberapa nilai, Anda dapat melakukan beberapa pertanyaan

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

Apa langkah paling efisien untuk mendapatkan fitur aktual milik ID fitur yang dikembalikan?

underdark
sumber

Jawaban:

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>
tuan
sumber
Terima kasih! Snorfalorpagus menyebutkan bahwa setFilterFids akan jauh lebih lambat daripada solusi yang ia posting. Apakah Anda mengonfirmasi hal ini?
underdark
Saya belum menggunakannya pada set hasil yang besar jadi tidak bisa mengkonfirmasi.
gsherman
1
Saya mengkonfirmasi dan, dalam kasus saya, rtree bahkan lebih cepat daripada QgsSpatialIndex () (untuk konstruksi Planar Graphs dari lapisan polylines yang sangat besar, transposisi modul PlanarGraph dengan Shapely in PyQGIS. Tetapi solusi dengan Fiona, Shapely dan rtree masih menjadi tercepat)
gen
1
Saya percaya pertanyaannya adalah tentang mendapatkan fitur aktual dari id fitur yang dikembalikan daripada kecepatan berbagai metode pengindeksan.
gsherman
7

Dalam sebuah posting blog tentang hal itu , Nathan Woodrow memberikan kode berikut:

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)

Ini membuat kamus yang memungkinkan Anda mencari QgsFeature dengan cepat menggunakan FID-nya.

Saya telah menemukan bahwa untuk lapisan yang sangat besar ini tidak terlalu praktis karena memerlukan banyak memori. Namun, penggunaan alternatif (akses acak dari fitur yang diinginkan) layer.getFeatures(QgsFeatureRequest().setFilterFid(fid))tampaknya sangat lambat. Saya tidak yakin mengapa ini terjadi, karena panggilan yang setara menggunakan binding OGR SWIG layer.GetFeature(fid)tampaknya jauh lebih cepat dari ini.

Snorfalorpagus
sumber
1
Menggunakan kamus itu sangat jauh lebih cepat daripada layer.getFeatures(QgsFeatureRequest().setFilterFid(fid)). Saya sedang mengerjakan layer dengan fitur 140k, dan total waktu untuk pencarian 140k berubah dari beberapa menit menjadi beberapa detik.
Håvard Tveite
5

Untuk perbandingan, lihat Penggabungan Tata Ruang yang Lebih Efisien dengan Python tanpa QGIS, ArcGIS, PostGIS, dll . Solusi yang disajikan penggunaan Python modul Fiona , Shapely dan RTree (Indeks Tata Ruang).

Dengan PyQGIS dan contoh yang sama dua layer, pointdan polygon:

masukkan deskripsi gambar di sini

1) Tanpa indeks spasial:

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) Dengan indeks spasial 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

Apa arti indeks ini?

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

Kesimpulan yang sama seperti dalam Tata Ruang yang Lebih Efisien bergabung dengan Python tanpa QGIS, ArcGIS, PostGIS, dll :

  • Tanpa dan indeks, Anda harus mengulangi semua geometri (poligon dan titik).
  • Dengan indeks spasial pembatas (QgsSpatialIndex ()), Anda beralih hanya melalui geometri yang memiliki peluang untuk berpotongan dengan geometri Anda saat ini ('filter' yang dapat menghemat banyak perhitungan dan waktu ...).
  • Anda juga dapat menggunakan modul Python indeks spasial lainnya ( rtree , Pyrtree atau Quadtree ) dengan PyQGIS seperti dalam Menggunakan indeks spasial QGIS untuk mempercepat kode Anda (dengan QgsSpatialIndex () dan rtree )
  • tetapi Indeks Spasial bukan tongkat ajaib. Ketika sebagian besar dataset harus diambil, Indeks Spasial tidak dapat memberikan manfaat kecepatan apa pun.

Contoh lain dalam GIS se: Bagaimana menemukan garis terdekat ke titik di QGIS? [duplikat]

gen
sumber
Terima kasih untuk semua penjelasan tambahannya. Pada dasarnya, solusi Anda menggunakan daftar alih-alih dict seperti yang dilakukan Snorfalorpagus. Jadi sepertinya tidak ada fungsi layer.getFeatures ([ids]) ...
underdark
Tujuan dari penjelasan ini adalah murni geometris dan sangat mudah untuk menambahkan fungsi layer.getFeatures ([ids]) seperti dalam Penggabungan Spasial yang Lebih Efisien dengan Python tanpa QGIS, ArcGIS, PostGIS, dll.
gen
0

Rupanya satu-satunya metode untuk mendapatkan kinerja yang baik adalah dengan menghindari atau menggabungkan panggilan ke layer.getFeatures (), bahkan jika filternya sesederhana fid.

Sekarang, inilah jebakannya: memanggil getFeatures itu mahal. Jika Anda menyebutnya pada layer vektor, QGIS akan diminta untuk mengatur koneksi baru ke penyimpanan data (penyedia lapisan), membuat beberapa permintaan untuk mengembalikan data, dan mengurai setiap hasil saat dikembalikan dari penyedia. Ini bisa lambat, terutama jika Anda bekerja dengan beberapa jenis lapisan jarak jauh, seperti tabel PostGIS melalui koneksi VPN.

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

evod
sumber