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:
- 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.
- Zapytaj ten indeks, aby uzyskać przecinające się identyfikatory.
- 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
.qix
jest MapServer / GDAL / OGR / SpatiaLite QuadTree indexTypeError: 'Polygon' object is not callable
z przykładu aktualizacji, ponieważ zastępujeszshape
funkcję, którą zaimportowałeśfor i, shape in enumerate(gl):
Odpowiedzi:
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.
źródło
Prawie go masz, ale popełniłeś mały błąd. Musisz użyć
intersection
metody na indeksie przestrzennym, zamiast przekazywać indeks dointersection
metody w zbuforowanym punkcie. Po znalezieniu listy obiektów, w których pola graniczne nakładają się, musisz sprawdzić, czy zbuforowany punkt faktycznie przecina geometrie.Jeśli chcesz znaleźć punkty, które znajdują się w minimalnej odległości od twojej klasy lądowej, możesz
distance
zamiast tego skorzystać z tej metody (zamień odpowiednią sekcję z poprzedniej).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:
źródło
Tak, to jest pomysł. Oto fragment tego samouczka na temat używania indeksu przestrzennego r-drzewa w Pythonie , używania foremnej, Fiony i geopandas:
źródło