Nuestros mejores programadores agotaron sus depósitos de café, por su búsqueda diariamente por la solución, hasta que Carlos encontró el resultado en GitLab por lo tanto en este momento la comparte con nosotros.
Solución:
Configuremos algunos datos de muestra:
library(sp)
library(spdep)
example(columbus)
plot(columbus)
Ahora obtendremos algunos puntos. Haga clic 20 veces en el mapa, algunas dentro y otras fuera de los polígonos:
pts = locator(20,type="p")
Convertir a tipo de datos espaciales:
spts = SpatialPoints(pts)
Ahora usa rgeos para calcular distancias punto-polígono, y tomar el mínimo para cada punto:
library(rgeos)
> apply(gDistance(spts, columbus,byid=TRUE),2,min)
1 2 3 4 5 6 7
0.25947439 0.03898401 1.27156515 1.50490316 0.00000000 0.00000000 0.00000000
8 9 10 11 12 13 14
0.00000000 0.00000000 0.31312329 0.26742466 0.53934325 0.07764322 0.03909773
15 16 17 18 19 20
0.11156343 0.29243322 0.08872334 0.00000000 0.00000000 0.00000000
Ahí tienes Esa es la distancia mínima de cada uno de los 20 puntos a cualquiera de los polígonos. Los ceros son para puntos dentro de uno de los polígonos.
Tenga en cuenta que debe usar datos en un sistema de coordenadas espaciales proyectadas y no en latitud y longitud.
Si tiene 17 000 puntos, es posible que desee hacer esto en subconjuntos más pequeños para evitar crear una matriz de 17 000 x 5000 si tiene 5000 polígonos. No dijiste.
Como se mencionó aquí, también se puede usar el geosphere::dist2Line para coordenadas no proyectadas (lat-long). A continuación se muestra un ejemplo:
library(sp)
library(geosphere)
# some country polygons to try on
data(wrld_simpl, package = "maptools")
wrld_subset <- wrld_simpl[[email protected]$ISO2 %in% c("RO","HU","AT","DE","FR"),]
# Generate random points (in and out)
set.seed(2017)
pts <- sp::makegrid(wrld_subset, n = 5)
# compute the shortest distance between points and polygons
# (from ?dist2Line): "returns matrix with distance and lon/lat of the nearest point" &
# "the ID (index) of (one of) the nearest objects"; distance is in meters (default)
dist.mat <- geosphere::dist2Line(p = pts, line = wrld_subset)
# bind results with original points
pts.wit.dist <- cbind(pts, dist.mat)
pts.wit.dist[1:3,]
## x1 x2 distance lon lat ID
## 1 -10.0 40.2 767133.6 -1.7808770 43.35992 2
## 2 -0.2 40.2 282022.2 0.1894124 42.71846 2
## 3 9.6 40.2 134383.0 9.1808320 41.36472 2
Graficar algunos resultados para tener una idea
pts.sp <- sp::SpatialPoints(coords = pts[,c("x1","x2")], # order matters
proj4string = [email protected])
plot(pts.sp, col="red")
plot(wrld_subset, add=TRUE)
# plot arrows to indicate the direction of the great-circle-distance
for (i in 1:nrow(pts.wit.dist))
arrows(x0 = pts.wit.dist[i,1],
y0 = pts.wit.dist[i,2],
x1 = pts.wit.dist[i,4],
y1 = pts.wit.dist[i,5],
length = 0.1,
col = "green")
Si aceptas, tienes la opción de dejar un tutorial acerca de qué le añadirías a esta reseña.