Menghitung jarak {minimal} antara poligon di R

9

Saya telah menghitung luas permukaan distribusi spesies (menggabungkan poligon dari shapefile), tetapi karena area ini dapat terdiri dari poligon yang cukup jauh, saya ingin menghitung beberapa ukuran dispersi. Apa yang telah saya lakukan sejauh ini adalah untuk mengambil centroid dari setiap poligon, menghitung jarak antara mereka dan menggunakannya untuk menghitung koefisien variasi, seperti pada contoh dummy di bawah ini;

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

Plot ketiga poligon dan centroidnya

masukkan deskripsi gambar di sini

Saya tidak yakin apakah pendekatan ini sangat berguna karena dalam banyak kasus, beberapa poligon (seperti yang biru dalam contoh di atas) cukup besar dibandingkan dengan yang lain, sehingga meningkatkan jarak lebih jauh. Misalnya centroid Australia memiliki jarak yang hampir sama dengan asrama sebelah barat seperti ke Papau.

Yang ingin saya dapatkan adalah beberapa masukan tentang pendekatan alternatif. Misalnya bagaimana atau dengan fungsi apa saya dapat menghitung jarak antar poligon?

Saya telah menguji untuk mengubah kerangka data SpatialPolygon di atas ke PointPatterns (ppp) {spatstat}untuk dapat berjalan nndist() {spatstat}untuk menghitung jarak antara semua titik. Tetapi karena saya berurusan dengan area yang cukup besar (banyak poligon dan yang besar), matriks menjadi besar dan saya tidak yakin bagaimana untuk terus sampai pada jarak minimal antara poligon .

Saya juga telah melihat fungsinya gDistance {rgeos}, tetapi saya pikir itu hanya bekerja pada data yang diproyeksikan yang dapat menjadi masalah bagi saya karena area saya dapat melewati beberapa EPSG areas. Masalah yang sama akan muncul untuk fungsi tersebut crossdist {spatstat}.

JO
sumber
1
Apakah Anda akan mempertimbangkan untuk menggunakan postgres/postgisselain R? Saya telah menggunakan alur kerja di mana saya melakukan sebagian besar pekerjaan saya R, tetapi menyimpan data dalam database yang saya akses menggunakan sqldf. Ini memungkinkan Anda untuk menggunakan semua postgisfungsi (yang jarak antar
poligonnya
@djq: Terima kasih telah berkomentar. Ya, saya pasti akan mencobanya :) Saya mulai membangun database postgrestetapi berhenti ketika saya tidak tahu (tidak melihat) bagaimana menghubungkan alur kerja / geostat antara database dan R...
jO.

Jawaban:

9

Anda dapat melakukan analisis ini dalam paket "spdep". Dalam fungsi tetangga yang relevan, jika Anda menggunakan "longlat = TRUE", fungsi menghitung jarak lingkaran yang hebat dan mengembalikan kilometer sebagai satuan jarak. Dalam contoh di bawah ini, Anda dapat memaksa objek daftar jarak yang dihasilkan ("dist.list") ke sebuah matriks atau data.frame, namun cukup efisien menghitung statistik ringkasan menggunakan 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
sumber
Terima kasih atas komentar dan wawasan tentang spdebpaket ini. Hanya untuk klarifikasi, pendekatan ini menghasilkan output yang sama seperti pada contoh saya, bukan?
jO.
Kalau-kalau Anda tidak melihat komentar saya di atas
jO.
Meskipun responsnya menyediakan kode yang berguna untuk menghitung jarak antara centroid, ia tidak berurusan dengan titik pusat OP yaitu bagaimana menemukan jarak antara dua titik terdekat dari batas poligon.
csfowler
Polisi yang hebat dan bentuk yang buruk untuk SE, tapi saya tidak bisa menyelesaikannya sekarang. Pencarian saya sendiri untuk jawaban atas pertanyaan ini tampaknya menunjukkan bahwa fungsi gDistance dari perpustakaan rgeos akan melakukan apa yang dimaksudkan OP: menemukan jarak terpendek antara tepi. Jika, dengan tergesa-gesa untuk memenuhi tenggat waktu yang ketat, saya telah salah menafsirkan OP atau Jeffrey Evans, permintaan maaf tulus saya.
csfowler