Zrozumienie zastosowania indeksów przestrzennych z RTree?

Mam problem ze zrozumieniem użycia indeksów przestrzennych w RTree.

Przykład: Mam 300 zbuforowanych punktów i muszę znać obszar przecięcia każdego bufora za pomocą wielokąta. Plik kształtu wielokąta ma> 20 000 wielokątów. Zasugerowano użycie wskaźników przestrzennych w celu przyspieszenia procesu.

SO ... Jeśli utworzę indeks przestrzenny dla mojego pliku kształtu wielokąta, czy zostanie on w jakiś sposób „dołączony” do pliku, czy też indeks będzie samodzielny? Czy po utworzeniu go mogę po prostu uruchomić funkcję przecięcia na pliku wielokąta i uzyskać szybsze wyniki? Czy skrzyżowanie „zobaczy”, że istnieją wskaźniki przestrzenne i będą wiedziały, co robić? Czy też muszę uruchomić go w indeksie, a następnie powiązać te wyniki z powrotem z moim oryginalnym plikiem wielokąta za pomocą FID lub innych podobnych?

Dokumentacja RTree niewiele mi pomaga (prawdopodobnie dlatego, że dopiero się uczę programowania). Pokazują, jak utworzyć indeks, wczytując ręcznie utworzone punkty, a następnie sprawdzając je względem innych ręcznie utworzonych punktów, które zwracają identyfikatory zawarte w oknie. Ma sens. Ale nie wyjaśniają, w jaki sposób odnosi się to do jakiegoś oryginalnego pliku, z którego pochodziłby indeks.

Myślę, że musi to wyglądać mniej więcej tak:

1. Wyciągnij boksy dla każdego elementu wielokąta z mojego pliku kształtu wielokąta i umieść je w indeksie przestrzennym, dając im identyfikator taki sam jak ich identyfikator w pliku kształtu.
2. Zapytaj ten indeks, aby uzyskać przecinające się identyfikatory.
3. Następnie ponownie uruchom moje skrzyżowanie tylko na funkcjach w moim oryginalnym pliku kształtu, które zostały zidentyfikowane przez zapytanie mojego indeksu (nie jestem pewien, jak to zrobię w tej ostatniej części).

Czy mam odpowiedni pomysł? Czy coś mi brakuje?

W tej chwili staram się, aby ten kod działał na jednym pliku kształtu punktowego, który zawiera tylko jeden element punktowy i jeden plik kształtu wielokąta, który zawiera> 20 000 elementów wielokąta.

Importuję pliki kształtów za pomocą Fiony, dodam indeks przestrzenny za pomocą RTree i próbuję wykonać skrzyżowanie za pomocą Shapely.

Mój kod testowy wygląda następująco:

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

Ale ciągle pojawia się błąd TypeError: Obiekt „Polygon” nie jest wywoływalny

To jest sedno tego. R-drzewo pozwala na wykonanie bardzo szybkiego pierwszego przejścia i daje zestaw wyników, które będą miały „fałszywie dodatnie” (obwiednie mogą się przecinać, gdy geometria dokładnie tego nie robi). Następnie przeglądasz zestaw kandydatów (pobierając ich z pliku kształtu według ich indeksu) i wykonujesz matematycznie dokładny test przecięcia, używając np. Shapely. Jest to ta sama strategia, która jest stosowana w przestrzennych bazach danych, takich jak PostGIS.

Prawie go masz, ale popełniłeś mały błąd. Musisz użyć intersectionmetody na indeksie przestrzennym, zamiast przekazywać indeks do intersectionmetody w zbuforowanym punkcie. Po znalezieniu listy obiektów, w których pola graniczne nakładają się, musisz sprawdzić, czy zbuforowany punkt faktycznie przecina geometrie.

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

Jeśli chcesz znaleźć punkty, które znajdują się w minimalnej odległości od twojej klasy lądowej, możesz distancezamiast tego skorzystać z tej metody (zamień odpowiednią sekcję z poprzedniej).

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

Jeśli zbudowanie indeksu przestrzennego zajmuje dużo czasu, a zamierzasz to zrobić więcej niż kilka razy, powinieneś rozważyć szeregowanie indeksu do pliku. Dokumentacja opisuje, jak to zrobić: http://toblerity.org/rtree/tutorial.html#serializing-your-index-to-a-file

Można również spojrzeć na ładowanie zbiorcze obwiedni do drzewa rtree za pomocą generatora, takiego jak ten:

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

Tak, to jest pomysł. Oto fragment tego samouczka na temat używania indeksu przestrzennego r-drzewa w Pythonie , używania foremnej, Fiony i geopandas:

Drzewo r reprezentuje pojedyncze obiekty i ich obwiednie („r” oznacza „prostokąt”) jako najniższy poziom indeksu przestrzennego. Następnie agreguje pobliskie obiekty i reprezentuje je za pomocą pola agregacji na kolejnym wyższym poziomie indeksu. Na jeszcze wyższych poziomach r-drzewo agreguje obwiednie i reprezentuje je po swojej obwiedni, iteracyjnie, aż wszystko zostanie zagnieżdżone w jednym obwiedni najwyższego poziomu. Aby wyszukać, r-drzewo pobiera pole zapytania i, zaczynając od najwyższego poziomu, widzi, które (jeśli w ogóle) pola ograniczające go przecinają. Następnie rozszerza każde przecinające się pole ograniczające i widzi, które z podrzędnych pól ograniczających w nim przecinają pole zapytania. Odbywa się to rekurencyjnie, aż wszystkie przecinające się pola zostaną przeszukane w dół do najniższego poziomu i zwróci pasujące obiekty z najniższego poziomu.