Cálculo da distância {mínima} entre polígonos em R

9

Eu calculei a área de superfície das distribuições de espécies (mesclando polígonos de shapefiles), mas como essa área pode ser composta por polígonos bastante distanciados, eu gostaria de calcular alguma medida de dispersão. O que fiz até agora é recuperar os centróides de cada polígono, calcular a distância entre eles e usá-los para calcular o coeficiente de variação, como no exemplo fictício abaixo;

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

Trama dos três polígonos e seus centróides

insira a descrição da imagem aqui

Não tenho certeza se essa abordagem é muito útil, pois em muitos casos, alguns dos polígonos (como o azul no exemplo acima) são bastante grandes em comparação com o restante, aumentando assim a distância ainda mais. Por exemplo, o centróide da Austrália tem quase a mesma distância das fronteiras ocidentais que de Papau.

O que eu gostaria de receber é uma contribuição sobre abordagens alternativas. Por exemplo, como ou com que função posso calcular a distância entre polígonos?

Eu testei para converter o quadro de dados SpatialPolygon acima em PointPatterns (ppp) {spatstat}para poder executar o nndist() {spatstat}cálculo da distância entre todos os pontos. Mas, como estou lidando com áreas muito grandes (muitos polígonos e áreas grandes), a matriz fica enorme e não sei como continuar a chegar à distância mínima entre os polígonos .

Também examinei a função gDistance {rgeos}, mas acho que ela só funciona com dados projetados, o que pode ser um problema para mim, pois minhas áreas podem atravessar várias EPSG areas. O mesmo problema surgiria para a função crossdist {spatstat}.

jO.
fonte
11
Você consideraria usar postgres/postgisalém de R? Eu usei um fluxo de trabalho no qual realizo a maior parte do meu trabalho R, mas armazeno os dados em um banco de dados ao qual estou acessando sqldf. Isso permite que você use todas as postgisfunções (das quais a distância entre polígonos é direta)
djq
@djq: Obrigado por comentar. Sim, eu definitivamente daria uma chance :) Comecei a construir um banco de dados, postgresmas parei quando não sabia (não parecia) como conectar o fluxo de trabalho / geoestados entre o banco de dados e R...
jO.

Respostas:

9

Você pode fazer essa análise no pacote "spdep". Nas funções vizinhas relevantes, se você usar "longlat = TRUE", a função calculará a grande distância do círculo e retornará quilômetros como a unidade de distância. No exemplo abaixo, você pode coagir o objeto resultante da lista de distâncias ("dist.list") a uma matriz ou data.frame, no entanto, é bastante eficiente calcular estatísticas de resumo 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
fonte
Obrigado por comentar e informações sobre o spdebpacote. Apenas para esclarecimento, essa abordagem produz a mesma saída que no meu exemplo, certo?
jO.
Caso você não tenha visto meu comentário acima
jO.
Embora a resposta forneça um código útil para calcular distâncias entre centróides, ela não lida com o ponto central do OP, que é como encontrar a distância entre os dois pontos mais próximos das bordas dos polígonos.
csfowler
Um policial enorme e péssimo para o SE, mas não posso fazer o trabalho completo agora. Minha própria busca por uma resposta para essa pergunta parece indicar que a função gDistance da biblioteca rgeos fará o que o OP pretendia: encontre a menor distância entre as bordas. Se, na pressa de cumprir um prazo apertado, interpretei mal o OP ou Jeffrey Evans, minhas sinceras desculpas.
csfowler