CODEWITHSUNDEEP
rmapspolygonvoronoispatstat
Bryan
Bryan

Reputation: 1821

Table of shared edges for a Voronoi tesslleation

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

Upvotes: 1

Views: 442

Answers (2)

Bryan
Bryan

Reputation: 1821

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

Upvotes: 2

Ege Rubak
Ege Rubak

Reputation: 4507

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

Upvotes: 0

Related Questions