У меня проблемы с пониманием использования пространственных индексов с RTree.
Пример: у меня есть 300 буферизованных точек, и мне нужно знать область пересечения каждого буфера с помощью полигонального шейп-файла. Шейп-файл полигонов имеет> 20000 полигонов. Было предложено использовать пространственные индексы для ускорения процесса.
ТАК ... Если я создам пространственный индекс для своего многоугольного шейп-файла, будет ли он каким-то образом "прикреплен" к файлу или индекс будет автономным? То есть после его создания я могу просто запустить функцию пересечения для файла полигонов и получить более быстрые результаты? Будет ли пересечение «видеть», что существуют пространственные индексы, и знать, что делать? Или мне нужно запустить его в индексе, а затем связать эти результаты с моим исходным файлом многоугольника через FID или что-то подобное?
Документация RTree мне не очень помогает (вероятно, потому что я только учусь программированию). Они показывают, как создать индекс, читая созданные вручную точки, а затем запрашивая его относительно других созданных вручную точек, который возвращает идентификаторы, содержащиеся в окне. Имеет смысл. Но они не объясняют, как это может относиться к какому-то оригинальному файлу, из которого пришел бы индекс.
Я думаю, что это должно пойти примерно так:
- Вытащите bbox для каждого полигонального объекта из моего шейп-файла полигонов и поместите их в пространственный индекс, предоставив им идентификатор, совпадающий с их идентификатором в шейп-файле.
- Запросите этот индекс, чтобы получить пересекающиеся идентификаторы.
- Затем повторно запустите мое пересечение только тех объектов в моем исходном шейп-файле, которые были идентифицированы путем запроса моего индекса (не уверен, как бы я выполнил эту последнюю часть).
У меня есть правильная идея? Я что-то упустил?
Сейчас я пытаюсь заставить этот код работать с одним точечным шейп-файлом, который содержит только один точечный объект, и одним многоугольным шейп-файлом, который содержит> 20000 многоугольных объектов.
Я импортирую шейп-файлы, используя Fiona, добавляю пространственный индекс, используя RTree, и пытаюсь сделать пересечение, используя Shapely.
Мой тестовый код выглядит так:
#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)
Но я продолжаю получать TypeError: объект 'Polygon' не вызывается
.qix
является картсервер / GDAL / OGR / SpatiaLite квадрантов индексTypeError: 'Polygon' object is not callable
пример обновления, потому что перезаписываетеshape
функцию, которую вы импортировали из shapely, объектом Polygon, который вы создаете с помощью этой строки:for i, shape in enumerate(gl):
Ответы:
Это суть этого. R-дерево позволяет сделать очень быстрый первый проход и дает набор результатов, которые будут иметь «ложные срабатывания» (ограничивающие рамки могут пересекаться, когда геометрия точно не пересекается). Затем вы просматриваете набор кандидатов (выбирая их из шейп-файла по их индексу) и выполняете математически точный тест пересечения с использованием, например, Shapely. Это та же самая стратегия, которая используется в пространственных базах данных, таких как PostGIS.
источник
Вы почти получили это, но сделали небольшую ошибку. Вам нужно использовать
intersection
метод для пространственного индекса, а не передавать индексintersection
методу в буферизованной точке. Как только вы нашли список объектов, в которых ограничивающие прямоугольники перекрываются, вам нужно проверить, действительно ли ваша буферизованная точка пересекает геометрию.Если вы заинтересованы в поиске точек, которые находятся на минимальном расстоянии от вашего класса земли, вы можете использовать
distance
метод вместо этого (поменяйте местами соответствующий раздел из предыдущего).Если создание пространственного индекса занимает много времени, и вы собираетесь делать это несколько раз, вам следует изучить сериализацию индекса в файл. Документация описывает, как это сделать: http://toblerity.org/rtree/tutorial.html#serializing-your-index-to-a-file
Вы также можете посмотреть на массовую загрузку ограничивающих рамок в дерево с помощью генератора, например:
источник
Да, это идея. Вот выдержка из этого руководства по использованию пространственного индекса r-дерева в Python , с использованием shapely, Fiona и геопанд:
источник