Calcular {mínima} distancia entre polígonos en R

9

He calculado el área de superficie de las distribuciones de especies (fusionando polígonos de archivos de forma), pero dado que esta área puede estar compuesta de polígonos bastante distantes, me gustaría calcular alguna medida de dispersión. Lo que he hecho hasta ahora es recuperar los centroides de cada polígono, calcular la distancia entre ellos y usarlos para calcular el coeficiente de variación, como en el siguiente ejemplo ficticio;

require(sp)
require(ggplot2)
require(mapdata)
require(gridExtra)
require(scales)
require(rgeos)
require(spatstat)

# Create the coordinates for 3 squares
ls.coords <- list()
ls.coords <- list()
ls.coords[[1]] <- c(15.7, 42.3, # a list of coordinates
                    16.7, 42.3,
                    16.7, 41.6,
                    15.7, 41.6,
                    15.7, 42.3)

ls.coords[[2]] <- ls.coords[[1]]+0.5 # use simple offset

ls.coords[[3]] <- c(13.8, 45.4, # a list of coordinates
                    15.6, 45.4,
                    15.6, 43.7,
                    13.8, 43.7,
                    13.8, 45.4)

# Prepare lists to receive the sp objects and data frames
ls.polys <- list()
ls.sp.polys <- list()

for (ii in seq_along(ls.coords)) {
   crs.args <- "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"
   my.rows <- length(ls.coords[[ii]])/2
   # create matrix of pairs
   my.coords <- matrix(ls.coords[[ii]],nrow = my.rows,ncol = 2,byrow = TRUE)
   # now build sp objects from scratch...
   poly = Polygon(my.coords)
   # layer by layer...
   polys = Polygons(list(poly),1)
   spolys = SpatialPolygons(list(polys))
   # projection is important
   proj4string(spolys) <- crs.args
   # Now save sp objects for later use
   ls.sp.polys[[ii]] <- spolys
   # Then create data frames for ggplot()
   poly.df <- fortify(spolys)
   poly.df$id <- ii
   ls.polys[[ii]] <- poly.df
}

# Convert the list of polygons to a list of owins
w <- lapply(ls.sp.polys, as.owin)
# Calculate the centroids and get the output to a matrix
centroid <- lapply(w, centroid.owin)
centroid <- lapply(centroid, rbind)
centroid <- lapply(centroid, function(x) rbind(unlist(x)))
centroid <- do.call('rbind', centroid)

# Create a new df and use fortify for ggplot
centroid_df <- fortify(as.data.frame(centroid))
# Add a group column
centroid_df$V3 <- rownames(centroid_df)

ggplot(data = italy, aes(x = long, y = lat, group = group)) +
  geom_polygon(fill = "grey50") +
  # Constrain the scale to 'zoom in'
  coord_cartesian(xlim = c(13, 19), ylim = c(41, 46)) +
  geom_polygon(data = ls.polys[[1]], aes(x = long, y = lat, group = group), fill = alpha("red", 0.3)) +
  geom_polygon(data = ls.polys[[2]], aes(x = long, y = lat, group = group), fill = alpha("green", 0.3)) +
  geom_polygon(data = ls.polys[[3]], aes(x = long, y = lat, group = group), fill = alpha("lightblue", 0.8)) + 
  coord_equal() +
  # Plot the centroids
  geom_point(data=centroid_points, aes(x = V1, y = V2, group = V3))

# Calculate the centroid distances using spDists {sp}
centroid_dists <- spDists(x=centroid, y=centroid, longlat=TRUE)

centroid_dists

       [,1]      [,2]     [,3]
[1,]   0.00000  69.16756 313.2383
[2,]  69.16756   0.00000 283.7120
[3,] 313.23834 283.71202   0.0000

# Calculate the coefficient of variation as a measure of polygon dispersion 
cv <- sd(centroid_dist)/mean(centroid_dist)
[1] 0.9835782

Gráfico de los tres polígonos y sus centroides.

ingrese la descripción de la imagen aquí

No estoy seguro de si este enfoque es muy útil ya que en muchos casos, algunos de los polígonos (como el azul en el ejemplo anterior) son bastante grandes en comparación con el resto, lo que aumenta aún más la distancia. Por ejemplo, el centroide de Australia tiene casi la misma distancia a sus residentes occidentales que a Papau.

Lo que me gustaría obtener es alguna información sobre enfoques alternativos. Por ejemplo, ¿cómo o con qué función puedo calcular la distancia entre polígonos?

He probado para convertir el marco de datos SpatialPolygon anterior a PointPatterns (ppp) {spatstat}para poder ejecutar nndist() {spatstat}y calcular la distancia entre todos los puntos. Pero como estoy tratando con áreas bastante grandes (muchos polígonos y grandes), la matriz se vuelve enorme y no estoy seguro de cómo continuar para llegar a la distancia mínima entre los polígonos .

También he analizado la función gDistance {rgeos}, pero creo que solo funciona en datos proyectados, lo que puede ser un problema para mí, ya que mis áreas pueden cruzar varios EPSG areas. El mismo problema surgiría para la función crossdist {spatstat}.

jO.
fuente
1
¿Consideraría usar postgres/postgisademás de R? He usado un flujo de trabajo donde realizo la mayor parte de mi trabajo R, pero almaceno los datos en una base de datos a la que accedo usando sqldf. Esto le permite utilizar todas las postgisfunciones (de las cuales la distancia entre polígonos es directa)
djq
@djq: Gracias por comentar. Sí, definitivamente lo probaría :) Empecé a construir una base de datos, postgrespero me detuve cuando no sabía (no veía) cómo conectar el flujo de trabajo / geostats entre la base de datos y R...
jO.

Respuestas:

9

Puede hacer este análisis en el paquete "spdep". En las funciones vecinas relevantes, si usa "longlat = TRUE", la función calcula la distancia del gran círculo y devuelve kilómetros como la unidad de distancia. En el ejemplo a continuación, podría forzar el objeto de la lista de distancia resultante ("dist.list") a una matriz o data.frame, sin embargo, es bastante eficiente calcular estadísticas de resumen usando lapply.

require(sp)
require(spdep)

# Create SpatialPolygonsDataFrame for 3 squares
poly1 <- Polygons(list(Polygon(matrix(c(15.7,42.3,16.7,42.3,16.7,41.6,15.7,41.6,15.7,42.3), 
                   nrow=5, ncol=2, byrow=TRUE))),"1")     
poly2 <- Polygons(list(Polygon(matrix(c(15.7,42.3,16.7,42.3,16.7,41.6,15.7,41.6,15.7,42.3)+0.5, 
                   nrow=5, ncol=2, byrow=TRUE))),"2")     
poly3 <- Polygons(list(Polygon(matrix(c(13.8, 45.4, 15.6, 45.4,15.6, 43.7,13.8, 43.7,13.8, 45.4), 
                   nrow=5, ncol=2, byrow=TRUE))),"3")                      
spolys = SpatialPolygons(list(poly1,poly2,poly3),1:3)
 spolys <- SpatialPolygonsDataFrame(spolys, data.frame(ID=sapply(slot(spolys, "polygons"), 
                                    function(x) slot(x, "ID"))) )   
   proj4string(spolys) <- "+proj=longlat +datum=WGS84 +no_defs +ellps=WGS84 +towgs84=0,0,0"

# Centroid coordinates (not used but provided for example) 
coords <- coordinates(spolys)

# Create K Nearest Neighbor list
skNN.nb <- knn2nb(knearneigh(coordinates(spolys), longlat=TRUE), 
                  row.names=spolys@data$ID)

# Calculate maximum distance for all linkages 
maxDist <- max(unlist(nbdists(skNN.nb, coordinates(spolys), longlat=TRUE)))

# Create spdep distance object
sDist <- dnearneigh(coordinates(spolys), 0, maxDist^2, row.names=spolys@data$ID)
  summary(sDist, coordinates(spolys), longlat=TRUE)

# Plot neighbor linkages                  
plot(spolys, border="grey") 
  plot(sDist, coordinates(spolys), add=TRUE)  

# Create neighbor distance list 
( dist.list <- nbdists(sDist, coordinates(spolys), longlat=TRUE) )

# Minimum distance 
( dist.min <- lapply(dist.list, FUN=min) )

# Distance coefficient of variation    
( dist.cv <- lapply(dist.list, FUN=function(x) { sd(x) / mean(x) } ) )
Jeffrey Evans
fuente
Gracias por comentar y la información sobre el spdebpaquete. Solo para aclarar, este enfoque produce el mismo resultado que en mi ejemplo, ¿verdad?
jO.
En caso de que no hayas visto mi comentario anterior
jO.
Aunque la respuesta proporciona un código útil para calcular distancias entre los centroides, no trata con el punto central del OP, que es cómo encontrar la distancia entre los dos puntos más cercanos de los bordes del polígono.
csfowler
Una gran policía y mala forma para SE, pero no puedo hacer todo el trabajo en este momento. Mi propia búsqueda de una respuesta a esta pregunta parece indicar que la función gDistance de la biblioteca rgeos hará lo que pretendía el OP: encontrar la distancia más corta entre los bordes. Si, en mi apuro por cumplir un plazo ajustado, he malinterpretado a OP o Jeffrey Evans mis sinceras disculpas.
csfowler