Memahami penggunaan indeks spasial dengan RTree?

Saya mengalami kesulitan memahami penggunaan indeks spasial dengan RTree.

Contoh: Saya memiliki 300 titik buffer, dan saya perlu tahu area persimpangan masing-masing buffer dengan bentuk poligon. Shapefile poligon memiliki> 20.000 poligon. Disarankan agar saya menggunakan indeks spasial untuk mempercepat proses.

JADI ... Jika saya membuat indeks spasial untuk bentuk poligon saya, apakah akan "dilampirkan" ke file dengan cara tertentu, atau akankah indeks berdiri sendiri? Yaitu, setelah membuatnya, bisakah saya menjalankan fungsi persimpangan pada file poligon dan mendapatkan hasil yang lebih cepat? Akankah persimpangan "melihat" bahwa ada indeks spasial dan tahu apa yang harus dilakukan? Atau, apakah saya harus menjalankannya di indeks, lalu menghubungkan kembali hasil-hasil itu ke file poligon asli saya melalui FID atau semacamnya?

Dokumentasi RTree tidak banyak membantu saya (mungkin karena saya baru belajar pemrograman). Mereka menunjukkan cara membuat indeks dengan membaca poin yang dibuat secara manual, dan kemudian menanyakannya terhadap poin lain yang dibuat secara manual, yang mengembalikan id yang terkandung dalam jendela. Masuk akal. Tetapi, mereka tidak menjelaskan bagaimana hal itu akan berhubungan kembali dengan beberapa file asli dari mana indeks tersebut berasal.

Saya pikir itu harus seperti ini:

1. Tarik kotak untuk setiap fitur poligon dari shapefile poligon saya dan letakkan di indeks spasial, beri mereka id yang sama dengan id mereka di shapefile.
2. Permintaan indeks itu untuk mendapatkan id yang berpotongan.
3. Kemudian jalankan kembali persimpangan saya hanya pada fitur-fitur di shapefile asli saya yang diidentifikasi dengan menanyakan indeks saya (tidak yakin bagaimana saya akan melakukan bagian terakhir ini).

Apakah saya punya ide yang tepat? Apakah saya kehilangan sesuatu?

Saat ini saya sedang berusaha agar kode ini berfungsi pada satu titik shapefile yang hanya berisi satu fitur titik, dan satu poligon bentuk yang berisi> 20.000 fitur poligon.

Saya mengimpor shapefile menggunakan Fiona, menambahkan indeks spasial menggunakan RTree, dan mencoba melakukan persimpangan menggunakan Shapely.

Kode pengujian saya terlihat seperti ini:

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

Tapi saya terus mendapatkan TypeError: objek 'Polygon' tidak bisa dipanggil

Itulah intinya. R-tree memungkinkan Anda untuk membuat operan pertama yang sangat cepat dan memberi Anda satu set hasil yang akan memiliki "false positive" (kotak pembatas dapat berpotongan ketika geometri tepatnya tidak). Kemudian Anda memeriksa set kandidat (mengambil mereka dari shapefile dengan indeks mereka) dan melakukan tes persimpangan matematis yang tepat menggunakan, misalnya, Shapely. Ini adalah strategi yang sama yang digunakan dalam database spasial seperti PostGIS.

Anda hampir mendapatkannya, tetapi Anda membuat kesalahan kecil. Anda harus menggunakan intersectionmetode pada indeks spasial, daripada meneruskan indeks ke intersectionmetode pada titik buffered. Setelah Anda menemukan daftar fitur di mana kotak pembatas tumpang tindih, maka Anda perlu memeriksa apakah titik buffer Anda benar-benar memotong geometri.

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

Jika Anda tertarik untuk menemukan poin yang berada dalam jarak minimum ke kelas tanah Anda, Anda bisa menggunakan distancemetode ini (menukar bagian yang sesuai dari sebelumnya).

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!')

Jika butuh waktu lama untuk membangun indeks spasial Anda dan Anda akan melakukan ini lebih dari beberapa kali, Anda harus melihat serialisasi indeks ke file. Dokumentasi menjelaskan cara melakukan ini: http://toblerity.org/rtree/tutorial.html#serializing-your-index-to-a-file

Anda juga dapat melihat pemuatan massal kotak pembatas ke rtree menggunakan generator, seperti ini:

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

Ya, itulah idenya. Berikut adalah kutipan dari tutorial ini tentang penggunaan indeks spasial r-tree dalam Python , menggunakan shapely, Fiona, dan geopanda:

R-tree mewakili objek individual dan kotak pembatasnya ("r" adalah untuk "persegi panjang") sebagai level terendah dari indeks spasial. Kemudian agregat objek terdekat dan mewakili mereka dengan kotak pembatas agregat mereka di tingkat indeks berikutnya yang lebih tinggi. Pada tingkat yang lebih tinggi, r-tree mengumpulkan kotak yang mengikat dan mewakilinya dengan kotak pembatas mereka, secara iteratif, sampai semuanya bersarang menjadi satu kotak pembatas tingkat atas. Untuk mencari, r-tree mengambil kotak kueri dan, mulai dari tingkat atas, melihat kotak pembatas mana (jika ada) yang memotongnya. Itu kemudian memperluas setiap kotak lompatan berpotongan dan melihat yang mana dari kotak lompatan anak di dalamnya memotong kotak permintaan. Ini berlangsung secara rekursif sampai semua kotak berpotongan dicari ke tingkat terendah, dan mengembalikan objek yang cocok dari tingkat terendah.