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.
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}
.
postgres/postgis
además deR
? He usado un flujo de trabajo donde realizo la mayor parte de mi trabajoR
, pero almaceno los datos en una base de datos a la que accedo usandosqldf
. Esto le permite utilizar todas laspostgis
funciones (de las cuales la distancia entre polígonos es directa)postgres
pero me detuve cuando no sabía (no veía) cómo conectar el flujo de trabajo / geostats entre la base de datos yR
...Respuestas:
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.
fuente
spdeb
paquete. Solo para aclarar, este enfoque produce el mismo resultado que en mi ejemplo, ¿verdad?