Zrozumienie zastosowania indeksów przestrzennych z RTree?

13

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

CSB
źródło
1
Indeks przestrzenny jest integralny i przejrzysty dla zestawu danych (zawarty, nie jest pojedynczym bytem z perspektywy użytkownika) Oprogramowanie, które wykonuje skrzyżowania, jest świadome i wykorzysta wskaźniki przestrzenne, aby utworzyć krótką listę do wykonania prawdziwego skrzyżowania, szybko informując oprogramowanie, którego funkcje należy rozważyć w celu dokładniejszej inspekcji i które wyraźnie nie są przecinające się. Sposób jego utworzenia zależy od oprogramowania i typu danych ... proszę podać więcej informacji na temat tych punktów, aby uzyskać bardziej szczegółową pomoc. W przypadku pliku kształtu jest to plik .shx.
Michael Stimson,
4
.shx nie jest indeksem przestrzennym. Jest to po prostu plik przesunięcia dynamicznego dostępu do rekordu o zmiennej szerokości. .sbn / .sbx to para indeksów przestrzennych Shapefile typu ArcGIS, chociaż ich specyfikacja nie została wydana.
Vince
1
Również .qixjest MapServer / GDAL / OGR / SpatiaLite QuadTree index
Mike T
Twój pomysł jest idealnie odpowiedni dla Spatialite, który nie ma prawdziwego indeksu przestrzennego. Większość innych formatów, jeśli w ogóle obsługują indeksy przestrzenne, robi to w przejrzysty sposób.
user30184,
2
Ciągle korzystasz TypeError: 'Polygon' object is not callablez przykładu aktualizacji, ponieważ zastępujesz shapefunkcję, którą zaimportowałeś for i, shape in enumerate(gl):
foremnie

Odpowiedzi:

12

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.

sgillies
źródło
1
Fajna gra słów (GiST)! GiST jest zwykle opisywany jako wariant B-drzewa, ale Postgresql ma implementację Gi-R dla drzewa. Chociaż wiki niekoniecznie jest najlepszym odniesieniem do cytowania, ma ładny schemat wyjaśniający wyszukiwania w obwiedni.
MappaGnosis,
Warto nauczyć się ręcznego korzystania z indeksu R-drzewa, jak w krokach 2 i 3. Ten blog o OGC GeoPackage, który obsługuje również R-drzewa, ponieważ osobne tabele bazy danych pokazują niektóre SQL i zrzuty ekranu openjump.blogspot.fi / 2014/02 /… .
user30184,
9

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))
Snorfalorpagus
źródło
2

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.

eos
źródło