Indeks spasial RTree tidak menghasilkan perhitungan persimpangan yang lebih cepat

Saya punya beberapa kode yang saya gunakan untuk menentukan Shapely Polygon / MultiPolygons berpotongan dengan sejumlah Shapely LineStrings. Melalui jawaban atas pertanyaan ini , kodenya telah berubah dari ini:

import fiona
from shapely.geometry import LineString, Polygon, MultiPolygon, shape

# Open each layer
poly_layer = fiona.open('polygon_layer.shp')
line_layer = fiona.open('line_layer.shp')

# Convert to lists of shapely geometries
the_lines = [shape(line['geometry']) for line in line_layer]
the_polygons = [(poly['properties']['GEOID'], shape(poly['geometry'])) for poly in poly_layer]

# Check for Polygons/MultiPolygons that the LineString intersects with
covered_polygons = {}
for poly_id, poly in the_polygons:
    for line in the_lines:
        if poly.intersects(line):
            covered_polygons[poly_id] = covered_polygons.get(poly_id, 0) + 1

di mana setiap persimpangan yang mungkin diperiksa, untuk ini:

import fiona
from shapely.geometry import LineString, Polygon, MultiPolygon, shape
import rtree

# Open each layer
poly_layer = fiona.open('polygon_layer.shp')
line_layer = fiona.open('line_layer.shp')

# Convert to lists of shapely geometries
the_lines = [shape(line['geometry']) for line in line_layer]
the_polygons = [(poly['properties']['GEOID'], shape(poly['geometry'])) for poly in poly_layer]

# Create spatial index
spatial_index = rtree.index.Index()
for idx, poly_tuple in enumerate(the_polygons):
    _, poly = poly_tuple
    spatial_index.insert(idx, poly.bounds)

# Check for Polygons/MultiPolygons that the LineString intersects with
covered_polygons = {}
for line in the_lines:
    for idx in list(spatial_index.intersection(line.bounds)):
        if the_polygons[idx][1].intersects(line):
            covered_polygons[idx] = covered_polygons.get(idx, 0) + 1

di mana indeks spasial digunakan untuk mengurangi jumlah pemeriksaan persimpangan.

Dengan shapefile yang saya miliki (sekitar 4000 poligon, dan 4 baris), kode asli melakukan 12936 .intersection()pemeriksaan dan membutuhkan waktu sekitar 114 detik untuk dijalankan. Sepotong kode kedua yang menggunakan indeks spasial hanya melakukan 1816 .intersection()pemeriksaan tetapi juga membutuhkan waktu sekitar 114 detik untuk dijalankan.

Kode untuk membangun indeks spasial hanya membutuhkan waktu 1-2 detik untuk dijalankan, sehingga pemeriksaan 1816 pada potongan kode kedua menghabiskan waktu yang hampir sama untuk melakukan sebagaimana 12936 memeriksa dalam kode asli (sejak pemuatan shapefile dan mengkonversi ke geometri Shapely adalah sama di kedua bagian kode).

Saya tidak dapat melihat alasan mengapa indeks spasial akan membuat .intersects()cek lebih lama, jadi saya bingung mengapa ini terjadi.

Saya hanya bisa berpikir bahwa saya menggunakan indeks spasial RTree secara tidak benar. Pikiran?

Jawaban saya pada dasarnya didasarkan pada jawaban lain oleh @ gen di sini:

Penggabungan Spasial yang Lebih Efisien dengan Python tanpa QGIS, ArcGIS, PostGIS, dll

Dia mengusulkan solusi yang sama menggunakan dua metode berbeda, dengan atau tanpa indeks spasial.

Dia (dengan benar) menyatakan:

Apa bedanya ?

• Tanpa indeks, Anda harus mengulangi semua geometri (poligon dan titik).
• Dengan indeks spasial pembatas ( Indeks Tata Ruang RTree ), Anda beralih hanya melalui geometri yang memiliki peluang untuk berpotongan dengan geometri Anda saat ini ('filter' yang dapat menghemat banyak perhitungan dan waktu ...)
• tetapi Indeks Spasial bukan tongkat ajaib. Ketika sebagian besar dataset harus diambil, Indeks Spasial tidak dapat memberikan manfaat kecepatan apa pun.

Kalimat-kalimat ini cukup jelas, tetapi saya lebih suka mengutip @gene daripada mengusulkan kesimpulan yang sama seperti milik saya (jadi, semua kreditnya adalah hasil kerjanya yang brilian!).

Untuk pemahaman yang lebih baik tentang indeks spasial Rtree, Anda dapat mengambil beberapa informasi berguna berikut tautan ini:

• Rtree: Pengindeksan spasial untuk Python ;
• Menggunakan Rtree indexing spasial dengan OGR .

Pengantar hebat lainnya tentang penggunaan indeks spasial mungkin artikel ini oleh @Nathan Woodrow .

Tambahkan saja ke mgri jawab.

Penting untuk memahami apa itu indeks spasial ( Bagaimana Cara Menerapkan Kotak Bounding untuk Shapely & Fiona? ). Dengan contoh saya di Bagaimana cara efisien menentukan mana dari ribuan poligon bersinggungan dengan linestring

Anda dapat membuat indeks spasial dengan poligon

idx = index.Index()
for feat in poly_layer:
    geom = shape(feat['geometry'])
    id = int(feat['id'])
    idx.insert(id, geom.bounds,feat)

Batas indeks spasial (batas poligon berwarna hijau)

Atau dengan LineStrings

  idx = index.Index()
  for feat in line_layer:
      geom = shape(feat['geometry'])
      id = int(feat['id'])
      idx.insert(id, geom.bounds,feat)

Batas indeks spasial (LineString terikat warna merah)

Sekarang, Anda beralih hanya melalui geometri yang memiliki kesempatan untuk berpotongan dengan geometri Anda saat ini (berwarna kuning)

Saya menggunakan di sini indeks spasial LineStrings (hasilnya sama tetapi dengan contoh Anda 4000 poligon dan 4 baris ....).

for feat1 in poly_layer:
    geom1 = shape(feat1['geometry'])
    for id in idx.intersection(geom1.bounds):
        feat2 = line_layer[id]
        geom2 = shape(feat2['geometry'])
        if geom2.intersects(geom1):
            print 'Polygon {} intersects line {}'.format(feat1['id'], feat2['id'])

  Polygon 2 intersects line 0
  Polygon 3 intersects line 0
  Polygon 6 intersects line 0
  Polygon 9 intersects line 0

Anda juga dapat menggunakan generator ( example.py )

def build_ind():
     with fiona.open('polygon_layer.shp') as shapes:
         for s in shapes:
             geom = shape(s['geometry'])
             id = int(s['id'])
             yield (id, geom.bounds, s)

 p= index.Property()
 tree = index.Index(build_ind(), properties=p)
 # first line of line_layer
 first = shape(line_layer.next()['geometry'])
 # intersection of first with polygons
 tuple(tree.intersection(first.bounds))
 (6, 2, 3, 9)

Anda dapat memeriksa skrip GeoPandas sjoin.py untuk memahami penggunaan Rtree .

Ada banyak solusi tetapi jangan lupakan itu

• Indeks Spasial bukan tongkat sihir ...

Sunting: Untuk memperjelas jawaban ini, saya salah meyakini bahwa semua tes persimpangan memakan waktu yang kira-kira sama. Ini bukan kasusnya. Alasan saya tidak mendapatkan kecepatan yang saya harapkan dari menggunakan indeks spasial adalah bahwa pemilihan tes persimpangan adalah yang paling lama dilakukan di tempat pertama.

Seperti yang dikatakan oleh gen dan mgri, indeks spasial bukanlah tongkat ajaib. Meskipun indeks spasial mengurangi jumlah tes persimpangan yang perlu dilakukan dari 12936 menjadi hanya 1816, tes 1816 adalah tes yang mengambil sebagian besar waktu untuk menghitung di tempat pertama.

Indeks spasial digunakan dengan benar, tetapi asumsi bahwa setiap tes persimpangan memakan waktu yang kira-kira sama adalah yang salah. Waktu yang dibutuhkan oleh uji persimpangan dapat sangat bervariasi (0,05 detik versus 0,000007 detik).