Table of shared edges for a Voronoi tesslleation

Question

I'm attempting to create a table of polygon neighbors based on a Voronoi tessellation (aka Dirichlet tessellations or Thiessen polygons) generated by the dirichlet() function of the spatstat library. For example, in the figure below, the upper right and bottom right tiles would each have 2 neighbors, the center right tile 4 neighbors, and the remaining two tiles 3 neighbors each. I want to capture the neighbor pairs in a table, and ideally capture the length of the boundary line that they share: e.g., 'Tile1', 'Tile2', 'shared_edge_length'.

Initially, I attempted to loop through and compare each pair of polygons in the tessellation using intersect.tess(), intersect.own(), as well as polyclip functions, but I'm guessing these don't work because the tiles are by definition not overlapping in area, despite sharing edges. Is there a simple function to accomplish this (the alternative being to loop through the $bdry points perhaps)? It seems like the regeos package has gTouches, but I couldn't find something similar for spatstat.

This is my current non-working approach:

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
} } }

enter image description here

score 2 · Accepted Answer · answered Aug 11 '14 at 19:03

The following function was provided by the package authors. It uses the fact that the deldir() function's dirsgs structure outputs the start/end coordinates of each line in the tessellation along with the point indices. These can be converted to a psp line segment pattern which can easily provide the length of each segment using lengths.psp(). The code below produces a table with one row for each of the 7 edges that can be seen in the plot above.

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)

The output stored in 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

score 0 · Answer 2 · answered Jun 20 '14 at 10:12

0

I just asked Rolf Turner who is author of the package deldir, and he provided the code below for finding the number of neighbours for each tile. I haven't yet looked at your wish to "ideally capture the length of the boundary line that they share".

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

answered Jun 20 '14 at 10:12

Ege Rubak

4,347
1
10
18

After a more detailed reading of the question, I see that this doesn't provide an answer at all. It is just the number of neighbours and not indicating which pairs are neighbours. Hopefully, I will have time to get back to this later. (You could try `delaunay` for finding neighbouring pairs, but I'm absolutely not sure it is useful. Conceptually the Delaunay graph is the relevant thing to look at.). – Ege Rubak Jun 20 '14 at 13:39
thanks for the notes. I've been exploring your suggestion but haven't found a way to employ a Delaunay graph. Is there not a way to do this using pairs of the tiles already created with `dirichlet`? – Bryan Aug 06 '14 at 17:22

Table of shared edges for a Voronoi tesslleation

2 Answers2

Linked