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

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

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

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

Ответы

Ответ 1

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

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 )

enter image description here

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

Ответ 2

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

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

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

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