0

Briefly, I am going to do the Voronoi tessellation of 100 points, and create different sets of 100 points 1000 times and do the tessellation for each set of the points.

I created points as,

x=matrix(runif(100*1000),nrow=100,ncol=1000)
y=matrix(runif(100*1000),nrow=100,ncol=1000)

Basic code to perform Voronoi tessellation using spatstat package is

dir = ppp(x=x, y=y, window = square(c(0,1)))
tess = dirichlet(dir)
plot(tess, main = "")
points(dir, pch=19, cex = 0.5)

However, I need to do the Voronoi tessellation for 1000 samples and trying to create a loop. I want to choose each column of x and y and end up with 1000 'dir'. Then do the tessellation 'tess' for 1000 'dir'. Also I need to calculate the areas of voronoi cells using the function area=tile.areas(tess)

The loop that I created was

for (i in 1000) {
  dir[i] = ppp(x=x[,i], y=y[,i], window = square(c(0,1)))
}

but all I get is errors and warnings. Do you have any idea how to do it?

Fatih
  • 5
  • 1

1 Answers1

0

You need to store the output in an object, in this case, let's put it in a list called dirList. Also, you have to specify a sequence to iterate over. for (i in 100) does nothing, it has to be for (i in 1:100)

library(deldir)
library(spatstat)

x <- matrix(runif(100 * 1000), nrow = 100, ncol = 1000)
y <- matrix(runif(100 * 1000), nrow = 100, ncol = 1000)

dirList <- list()

for (i in 1:1000){
  dirList[[i]] <- ppp(x = x[ , i], y = y[ , i], window = square(c(0, 1)))
}

Then to plot, you access the object with [[]]

tess = dirichlet(dirList[[1]])
plot(tess, main = "")

To the second part of your question, using your object tess for one set of points:

Credit to @DJack for the general process

library(deldir)
library(spatstat)
library(maptools)
library(rgeos)

x <- matrix(runif(100 * 1000), nrow = 100, ncol = 1000)
y <- matrix(runif(100 * 1000), nrow = 100, ncol = 1000)

x <- x[ , 1]
y <- y[ , 1]

dir = ppp(x=x, y=y, window = square(c(0,1)))
tess = dirichlet(dir)

t1 <- as(tess, "SpatialPolygons")
t2 <- unionSpatialPolygons(t1, ID=rep(1,length(t1)))
vor_border <- which(gRelate(t1,t2, byid=TRUE) %in% c("2FF11F212"))

par(mfcol=c(1,3))
plot(t1)
plot(t1[vor_border,], col="red")
# you could also invert the mask to get just the inner polygons
plot(t1[-vor_border,], col="lightblue")

enter image description here

Mako212
  • 6,787
  • 1
  • 18
  • 37
  • Thank you for the answer. I am very new to loops and the Voronoi. Is there any simpler way to detect the cells which intersect with the boundary? I am increasing the boundaries of the unit square by 0.1 and re-calculating the areas of the cells. Then I check which cells' areas increased, so that is automatically the edge cell. I used tripack as well but the function `voronoi.findrejectsites` but it requires a triangulation object and takes the boundary as convex hull of the points, not a specified boundary such as unit square. – Fatih Sep 21 '18 at 20:07
  • @Fatih [this link](https://gis.stackexchange.com/questions/76060/simple-way-to-detect-select-border-polygons-of-a-shapefile) seems to outline a good process, although you may need to convert your tessellation to a `SpatialPolygons` object to do it. – Mako212 Sep 21 '18 at 20:36
  • Here's actual code: https://stackoverflow.com/questions/19784793/looking-for-a-function-allowing-to-select-identify-polygon-that-share-a-line-seg – Mako212 Sep 21 '18 at 21:32
  • Thank you Mako212 for your answer. I got the idea but the code `t1 = as(tess, "SpatialPolygons")` always giving me an error saying `Error in as(tess, "SpatialPolygons") : no method or default for coercing “tess” to “SpatialPolygons”` and I looked for ways to deal with it but I couldn't – Fatih Sep 22 '18 at 15:11
  • @Fatih try again now. I missed a library `maptools` (many apologies), and also had a name conflict that might have caused an issue. Let me know if it works now! – Mako212 Sep 22 '18 at 20:08
  • Mako212, it works perfectly now. Thank you so much. I will go ahead and do further analysis. – Fatih Sep 22 '18 at 21:02
  • @Fatih Great, don't forget to click the checkmark on the answer indicating that it answered your question! – Mako212 Sep 23 '18 at 07:46