Corrispondenza spaziale di dataset grandi

Question

Ho un set di dati con circa 100000 punti e un altro set di dati con circa 3000 poligoni. Per ciascuno dei punti ho bisogno di trovare il poligono più vicino (corrispondenza spaziale). I punti all'interno di un poligono devono corrispondere a quel poligono.

Calcolare le distanze su tutte le coppie è possibile, ma richiede un po 'più del necessario. Esiste un pacchetto R che farà uso di un indice spaziale per questo tipo di problema di abbinamento?

Sono a conoscenza del pacchetto sp e della funzione over, ma la documentazione non indica nulla sugli indici.

Answer 1

Si può provare e utilizzare la funzione gDistance nel pacchetto rgeos per questo. Ad esempio, guarda l'esempio seguente, che ho rielaborato da questo old thread. Spero che sia d'aiuto.

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

Io non so nulla di R ma io offrirò una possibile soluzione usando PostGIS. Potresti essere in grado di caricare i dati in PostGIS ed elaborarli più velocemente di quanto tu possa utilizzare la sola R.

Date due tabelle planet_osm_point (80k righe) e planet_osm_polygon (30k righe), la seguente query esegue in circa 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 
); 

Il risultato è un'approssimazione basata sulla distanza tra il punto e sul centro punto del riquadro di delimitazione del poligono (non il punto centrale del poligono stesso). Con un po 'più di lavoro questa query può essere adattata per ottenere il poligono più vicino basato sul punto centrale del poligono stesso sebbene non venga eseguito altrettanto rapidamente.