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.