centroidal voronoi tesselation of a complex polygon in R

tesselate-polygon.r

library('deldir')	# voronoi stuff
library('rgdal')	# postgis and spatial stuff
library('sp') # spatial stuff
library('raster') # more spatial stuff
library('splancs') # csr() for generating random points

# function from http://carsonfarmer.com/2009/09/voronoi-polygons-with-r/
voronoipolygons = function(layer,bounding_box) {
    require(deldir)
    crds = layer@coords
    z = deldir(crds[,1], crds[,2],rw=bounding_box)
    w = tile.list(z)
    polys = vector(mode='list', length=length(w))
    require(sp)
    for (i in seq(along=polys)) {
        pcrds = cbind(w[[i]]$x, w[[i]]$y)
        pcrds = rbind(pcrds, pcrds[1,])
        polys[[i]] = Polygons(list(Polygon(pcrds)), ID=as.character(i))
    }
    SP = SpatialPolygons(polys)
    voronoi = SpatialPolygonsDataFrame(SP, data=data.frame(x=crds[,1], 
        y=crds[,2], row.names=sapply(slot(SP, 'polygons'), 
        function(x) slot(x, 'ID'))))
}

db="PG:dbname='???' user='???' password='???' host='localhost'"

# get the polygon to tesselate
polys = readOGR(db,"cincy_polygon")
Cincy = subset(polys,name=='Cincinnati')

# define a reasonable, simple, bounding box for the polygon to be tesselated
xmin = 697733.99
xmax = 727622.22
ymin = 4325711.21
ymax = 4344302.87
Cbounds = c(xmin,xmax,ymin,ymax)

# random points
x = rnorm(n=550,mean=714794,sd=750)
y = rnorm(n=550,mean=4331915,sd=750)
points = SpatialPoints(cbind(x,y))
# do the initial tesselation
tiles = voronoipolygons(points,Cbounds) * Cincy

for(i in 1:10000){
	print(i)
	# store these for rendering
	old_points = points
	old_tiles = tiles
	# get centroids
	points = SpatialPoints( coordinates(old_tiles) )
	tiles = voronoipolygons(points,Cbounds) * Cincy

# plot iterations?
#	png(paste0('tess_',i,'_.png'),width=1000,height=1000)
#		par(mar=rep(1, 4))
#		# plot the bits from the last iteration
#		plot(old_tiles, axes=FALSE, ann=FALSE,border='red')
#		points(old_points, pch=20, col='red')
#		# plot the new bits
#		plot(tiles,add=TRUE,border='black')
#		points(points,pch=20,col='black')
#		box()
#	dev.off()

}

# output to shapefile
writeOGR(obj=tiles, dsn='shapes', layer='tess10k', driver="ESRI Shapefile")