How to do Voronoi tessellation of 100 points on a unit square multiple times

Question

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?

Answer 1

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")