Блог

Я пытаюсь создать таблицу соседей полигонов на основе тесселяции Вороного (также известной как тесселяции Дирихле или полигоны Тиссена), сгенерированную dirichlet() функцией spatstat библиотеки. Например, на рисунке ниже у верхнего правого и нижнего правого тайлов будет по 2 соседа, у центрального правого тайла 4 соседа, а у оставшихся двух тайлов по 3 соседа у каждого. Я хочу захватить соседние пары в таблице и в идеале захватить длину граничной линии, которую они разделяют: например, ‘Tile1’, ‘Tile2’, ‘shared_edge_length’.

Первоначально я попытался перебрать и сравнить каждую пару полигонов в тесселяции, используя intersect.tess() , intersect.own() , а также polyclip функции, но я предполагаю, что это не сработает, потому что плитки по определению не перекрываются по площади, несмотря на общие ребра. Существует ли простая функция для выполнения этого (возможно, альтернативой является перебор $bdry точек)? Похоже regeos , что в пакете есть gTouches , но я не смог найти что-то подобное spatstat .

Это мой текущий нерабочий подход:

 library(spatstat)
points <- ppp(x=c(-77.308703, -77.256582, -77.290600,  -77.135668, -77.097144),
         y=c(39.288603, 39.147019, 39.372818, 39.401898, 39.689203),
         window=owin(xrange=c(-77.7,-77), yrange=c(39.1, 39.7)))
vt <- dirichlet(points) # Dirichlet tesselation
plot(vt)

tilesA <- tiles(vt)
n_tiles <- length(tilesA)
boundary_calcs <- data.frame('area1_id'=numeric(), 'area2_id'=numeric(), 'neighbor'=logical()) # Store boundary pairs
for (i in 1:n_tiles) {
  for (j in 1:n_tiles) {
    intersection <- intersect.owin(tilesA[[i]], tilesA[[j]], fatal=FALSE) # does not work
    if (!is.empty(intersection)) {
      boundary_calcs[nrow(boundary_calcs) 1, ] <- c(i, j, TRUE) # add to data table as new row
} } }
  

Следующая функция была предоставлена авторами пакета. Он использует тот факт, что структура deldir() функции dirsgs выводит начальные / конечные координаты каждой строки в тесселяции вместе с индексами точек. Они могут быть преобразованы в шаблон psp линейного сегмента, который может легко обеспечить длину каждого используемого сегмента lengths.psp() . Приведенный ниже код создает таблицу с одной строкой для каждого из 7 ребер, которые можно увидеть на графике выше.

 library(spatstat)
library(deldir)
points <- ppp(x=c(-77.308703, -77.256582, -77.290600,  -77.135668, -77.097144),
         y=c(39.288603, 39.147019, 39.372818, 39.401898, 39.689203),
         window=owin(xrange=c(-77.7,-77), yrange=c(39.1, 39.7)))

sharededge <- function(X) {
   verifyclass(X, "ppp")
   Y <- X[as.rectangle(X)]
   dX <- deldir(Y)
   DS <- dX$dirsgs
   xyxy <- DS[,1:4]
   names(xyxy) <- c("x0","y0","x1","y1")
   sX <- as.psp(xyxy,window=dX$rw)
   marks(sX) <- 1:nobjects(sX)
   sX <- sX[as.owin(X)]
   tX <- tapply(lengths.psp(sX), marks(sX), sum)
   jj <- as.integer(names(tX))
   ans <- data.frame(ind1=DS[jj,5], 
                     ind2=DS[jj,6], 
                     leng=as.numeric(tX))
   return(ans)
}

shared_edge_lengths <- sharededge(points)
  

Выходные данные, хранящиеся в shared_edge_lengths :

   ind1 ind2       leng
1    2    1 0.17387212
2    3    1 0.13444458
3    4    1 0.05791519
4    4    2 0.10039321
5    4    3 0.25842530
6    5    3 0.09818828
7    5    4 0.17162429
  

Я только что спросил Рольфа Тернера, который является автором пакета deldir , и он предоставил приведенный ниже код для определения количества соседей для каждой плитки. Я еще не рассмотрел ваше желание «идеально зафиксировать длину линии границы, которую они разделяют».

 x <- c(-77.308703, -77.256582, -77.290600,  -77.135668, -77.097144)
y <- c(39.288603, 39.147019, 39.372818, 39.401898, 39.689203)
rw <- c(-77.7,-77,39.1,39.7)
require(deldir)
dxy <- deldir(x,y,rw=rw)

dxy$summary$n.tside