3

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.

Community
  • 1
  • 1
user3122539
  • 41
  • 5
  • Also, if anyone has any ideas of how to do this even without the hole in the middle of the polygon then that'd be great too. Thanks – user3122539 Jan 08 '14 at 22:12
  • further, In a previous related [post](http://stackoverflow.com/questions/12156475/combine-voronoi-polygons-and-maps/) I think @Spacedman suggests that something like this can be achieved with dummy points, but I'm unsure how this could be done. Good Luck to anyone giving that a try – user3122539 Jan 12 '14 at 13:40

1 Answers1

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

Neal Alexander
  • 11
  • 2
  • 1
    Welcome at StackOverflow! Please consider including the essential parts of the linked resource into your answer. Link only answers are considered bad answers and will most likely get deleted fast. – VMai Jul 30 '14 at 13:18