Пространственное сопоставление больших наборов данных

Question

У меня есть набор данных, содержащий около 100000 точек и еще один набор данных с примерно 3000 полигонами. Для каждой из точек мне нужно найти ближайший многоугольник (пространственное совпадение). Точки внутри полигона должны совпадать с этим полигоном.

Вычисление расстояний между парами возможно, но занимает немного больше времени, чем необходимо. Есть ли R-пакет, который будет использовать пространственный индекс для такого рода проблем соответствия?

Я знаю о пакете sp и функции over, но в документации ничего не говорится об индексах.

Answer 1

Вы можете попробовать и использовать функцию gDistance в пакете rgeos. В качестве примера рассмотрим приведенный ниже пример, который я переработал из этого old thread. Надеюсь, поможет.

require(rgeos) 
require(sp) 

# Make some polygons 
grd <- GridTopology(c(1,1), c(1,1), c(10,10)) 
polys <- as.SpatialPolygons.GridTopology(grd) 

# Make some points and label with letter ID 
set.seed(1091) 
pts = matrix(runif(20 , 1 , 10) , ncol = 2) 
sp_pts <- SpatialPoints(pts) 
row.names(pts) <- letters[1:10] 

# Plot 
plot(polys) 
text(pts , labels = row.names(pts) , col = 2 , cex = 2) 
text(coordinates(polys) , labels = row.names(polys) , col = "#313131" , cex = 0.75) 

# Find which polygon each point is nearest 
cbind(row.names(pts) , apply(gDistance(sp_pts , polys , byid = TRUE) , 2 , which.min)) 
# [,1] [,2] 
#1 "a" "86" 
#2 "b" "54" 
#3 "c" "12" 
#4 "d" "13" 
#5 "e" "78" 
#6 "f" "25" 
#7 "g" "36" 
#8 "h" "62" 
#9 "i" "40" 
#10 "j" "55" 

Answer 2

Я ничего об R не знаю, но я буду предлагать один из возможных вариантов решения с использованием PostGIS. Вы можете загружать данные в PostGIS и обрабатывать их быстрее, чем вы можете использовать только R.

Учитывая две таблицы planet_osm_point (80k строк) и planet_osm_polygon (30k строк), следующий запрос выполняется в 30s вокруг

create table knn as 
select 
    pt.osm_id point_osm_id, 
    poly.osm_id poly_osm_id 
from planet_osm_point pt, planet_osm_polygon poly 
where poly.osm_id = (
    select p2.osm_id 
    from planet_osm_polygon p2 
    order by pt.way <-> p2.way limit 1 
); 

Результатом является приближение на основе расстояния между точкой и Centre- точка ограничивающей рамки многоугольника (а не центральная точка самого полигона). С немного большей работой этот запрос может быть адаптирован для получения ближайшего многоугольника на основе центральной точки самого многоугольника, хотя он не будет выполняться так быстро.