Finding neighbours from voronoi polygons within a set map area

Question

After a rather unsuccessfully written question, I hope this one is more clear and direct, and any help at all with it is much appreciated.

I want to create voronoi/thiessen polygons around a set of points within a given 'map' in order to determine which points are neighbours (share a boundary line) with each other within this given area.

library(sp); library(rgeos); library(deldir)

Considering the situation where I have 14 locations I am interested in:

x<-c(0.9,1.7,2.4,2.9,4.83, 0.73, 2.31, 3.69, 4.23, 2.86, 1.91, 4.32, 4.60, 1.82)
y<-c(1.9,0.9,2.8,1.9,1.81, 1.66, 4.54, 5.66, 1.99, 4.03, 4.32, 5.98, 5.56, 3.41)
crds<-cbind(x,y)

within the given polygon (map)

x.p<-c(0.1,0.1,3.5,3,5,1,6,6,0.1)
y.p<-c(0.1,5,4.8,1,5,5.5,6.5,1,0.1)
poly<-cbind(x.p,y.p)

and this polygon (map) has a set hole in it of:

x.h<-c(1,1.1,1.5,2.1,1.9,2.3,3,1)
y.h<-c(1,2.9,3.1,3,2.8,2.2,1.5,1)
hole<-cbind(x.h,y.h)

I now need to know the first order neighbours between each of the points of interest, but where the voronoi polygons around the point of interest cannot extend out of the map boundaries, or into/through the hole.

deldir(crds[,1],crds[,2]) 

simply gives the voronoi polygons (and first order neighbours) without these constraints.

For illustrative purposes only and to further explain what I mean with an example, if we were to plot the polygons in a usual way:

voronoipolygons <- function(crds) {
z <- deldir(crds[,1], crds[,2],rw=c(0,7,0,7))
w <- tile.list(z)
polys <- vector(mode='list', length=length(w))
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=(1:nrow(crds))[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'))))
return(voronoi)
}
SP<-voronoipolygons(crds[,1:2])
plot(SP)

and then if I plot my constraints over this to demonstrate what I mean

ybox<-xbox<-c(0,7)

polypath(c(poly[,1], NA, c(xbox, rev(xbox))), c(poly[,2], NA, rep(ybox, each=2)), col="light blue", rule="evenodd")
polygon(hole[,1],hole[,2],col="light blue")
text(crds[,1],crds[,2],1:nrow(crds))

I am wanting to have it so when I use the deldir command (or something similar?), '3' would not be classed as neighbours with '1' or '2', nor would '10' and '8' be classed as neighbours with each other etc.

Edit:

I have since found this. It sounds similar to what I need (with the 'extent' option) but was hoping to carry it out without the use of ArcGIS software.

Answer 1

The following seems to help for the simpler question (without the hole).

https://github.com/cran/deldir/blob/master/inst/code.discarded/triang.list.R.save

In particular, the 'tlist' object defined by line 24.

You may not want to include as neighbours those whose tiles only meet a long way outside the plotting area. In that case, it'd be easier to take as neighbours each pair of points with an entry in the 'ind1' and 'ind2' columns of the 'dirsgs' component of the 'deldir' output. Again, this is for the simpler version of the question (no holes).

Neal