Obliczanie {minimalnej} odległości między wielokątami w R.


9

Obliczyłem powierzchnię rozmieszczenia gatunków (łączenie wielokątów z plików kształtów), ale ponieważ obszar ten może składać się z dość odległych wielokątów, chciałbym obliczyć pewną miarę dyspersji. Do tej pory udało mi się pobrać centroidy każdego wielokąta, obliczyć odległość między nimi i użyć ich do obliczenia współczynnika zmienności, jak w poniższym przykładzie manekina;

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

Wykres trzech wielokątów i ich centroidów

wprowadź opis zdjęcia tutaj

Nie jestem pewien, czy to podejście jest bardzo przydatne, ponieważ w wielu przypadkach niektóre wielokąty (jak niebieski w powyższym przykładzie) są dość duże w porównaniu z resztą, co jeszcze bardziej zwiększa odległość. Np. Centroid Australii ma prawie taką samą odległość do swoich zachodnich granic jak Papau.

Chciałbym uzyskać wkład w alternatywne podejścia. Np. Jak lub przy pomocy jakiej funkcji mogę obliczyć odległość między wielokątami?

Testowałem, aby przekonwertować powyższą ramkę danych SpatialPolygon na PointPatterns (ppp), {spatstat}aby móc uruchomić, nndist() {spatstat}aby obliczyć odległość między wszystkimi punktami. Ale ponieważ mam do czynienia z dość dużymi obszarami (wiele wielokątów i duże), matryca staje się ogromna i nie jestem pewien, jak dalej docierać do minimalnej odległości między wielokątami .

Patrzyłem również na tę funkcję gDistance {rgeos}, ale myślę, że działa ona tylko na rzutowanych danych, co może być dla mnie problemem, ponieważ moje obszary mogą przekraczać kilka EPSG areas. Ten sam problem pojawiłby się dla funkcji crossdist {spatstat}.


1
Czy rozważalibyście użycie postgres/postgisoprócz R? Korzystałem z przepływu pracy, w którym wykonuję większość pracy R, ale przechowuję dane w bazie danych, z której korzystam sqldf. Umożliwia to korzystanie ze wszystkich postgisfunkcji (których odległość między wielokątami jest prosta)
djq

@djq: Dziękujemy za komentowanie. Tak, zdecydowanie postgresspróbuję :) Zacząłem budować bazę danych, ale przestałem, gdy nie wiedziałem (nie szukałem), jak połączyć przepływ pracy / geostaty między bazą danych a R...
JO.

Odpowiedzi:


9

Możesz wykonać tę analizę w pakiecie „spdep”. W przypadku funkcji sąsiada, jeśli użyjesz „longlat = TRUE”, funkcja oblicza wielką odległość koła i zwraca kilometry jako jednostkę odległości. W poniższym przykładzie można przymusić wynikowy obiekt listy odległości („listę odległości”) do matrycy lub ramki danych. Jest to jednak dość wydajne obliczanie statystyk podsumowujących za pomocą 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) } ) )

Dziękujemy za komentowanie i wgląd w spdebpakiet. Tylko dla wyjaśnienia, to podejście daje takie same wyniki, jak w moim przykładzie, prawda?
JO.

Na wypadek, gdybyś nie widział mojego powyższego komentarza
JO.

Chociaż odpowiedź zapewnia użyteczny kod do obliczania odległości między centroidami, nie zajmuje się ona centralnym punktem OP, czyli sposobem znajdowania odległości między dwoma najbliższymi punktami granic wielokąta.
csfowler

Ogromny głupiec i zła forma dla SE, ale nie mogę teraz wykonać całej pracy. Moje własne wyszukiwanie odpowiedzi na to pytanie wydaje się wskazywać, że funkcja gDistance z bibliotek rgeos zrobi to, co zamierzał OP: znajdź najkrótszą odległość między krawędziami. Jeśli w pośpiechu, aby dotrzymać napiętego terminu, źle zinterpretowałem OP lub Jeffrey'a Evansa, szczerze przepraszam.
csfowler
Korzystając z naszej strony potwierdzasz, że przeczytałeś(-aś) i rozumiesz nasze zasady używania plików cookie i zasady ochrony prywatności.
Licensed under cc by-sa 3.0 with attribution required.