14

I would like to combine Voronoi polygons with map, in order to use this later for spatial analysis. I have number of points and shapefile that i want to combine and then save as shapefile/spatial polygons. To get voronoi polygons i use function from this topic.

My code is as follows:

coords<-data.frame(LONG=c(16.9252,16.9363,16.9408,16.8720,16.9167,16.9461,16.9093,16.9457,16.9171,16.8506,16.9471,16.8723,16.9444,16.9212,16.8809,16.9191,16.8968,16.8719,16.9669,16.8845),
LAT=c(52.4064,52.4266,52.3836,52.3959,52.4496,52.3924,52.4012,52.3924,52.3777,52.4368,52.4574,52.3945,52.4572,52.3962,52.3816,52.3809,52.3956,52.3761,52.4236,52.4539))

My map is available here: https://docs.google.com/file/d/0B-ZJyVlQBsqlSURiN284dF9YNUk/edit

library(rgdal)
voronoipolygons <- function(x) {
  require(deldir)
  if (.hasSlot(x, 'coords')) {
    crds <- x@coords  
  } else crds <- x
  z <- deldir(crds[,1], crds[,2])
  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'))))
}

And my code to get voronoipolygons:

pzn.coords<-voronoipolygons(coords)
plot(pznall)
plot(pzn.coords,add=T)
points(coords$LONG,coords$LAT)

result: enter image description here

I wan to have this voronoi polygon inside my map as new spatialpolygon.

I would be grateful for anwsers!

Edit:

To be clear, i want to achieve something like this (this lines should be created from voronoi polygons):

enter image description here

Community
  • 1
  • 1
Maciej
  • 3,255
  • 1
  • 28
  • 43
  • 3
    Its not clear what you want. You have the voronoi polygons as a spatialpolygonsdataframe in the pzn.coords object. – Spacedman Aug 28 '12 at 10:30
  • The link you provide is an `.rdata` file. You can create a new `myvoronietc.rdata` file using `save()` , e.g. `save('pzn.coords','lotsa_other_data',file='myvoronietc.rdata')` . – Carl Witthoft Aug 28 '12 at 12:02
  • I want to combine both files, to have voronoi polygons inside my map, not to have this square boundaries of voronoi polygons but in my map boundaries. – Maciej Aug 28 '12 at 12:35
  • 1
    You have to extend the voronoi polygons to be bigger than your polygon, possibly by adding dummy points to deldir, and then use rgeos to clip. I'll do a full answer later maybe, or try it yourself. – Spacedman Aug 28 '12 at 12:41
  • Just a comment (Not an answer) There is a typo in the line: rw = as.numeric(t(bbox(pznall))) replace pznall with poly Thanks for the modification, Cheers, Jed – jed.a.long May 23 '13 at 20:23

1 Answers1

12

Slightly modified function, takes an additional spatial polygons argument and extends to that box:

voronoipolygons <- function(x,poly) {
  require(deldir)
  if (.hasSlot(x, 'coords')) {
    crds <- x@coords  
  } else crds <- x
  bb = bbox(poly)
  rw = as.numeric(t(bb))
  z <- deldir(crds[,1], crds[,2],rw=rw)
  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'))))

  return(voronoi)

}

Then do:

pzn.coords<-voronoipolygons(coords,pznall)
library(rgeos)
gg = gIntersection(pznall,pzn.coords,byid=TRUE)
plot(gg)

Note that gg is a SpatialPolygons object, and you might get a warning about mismatched proj4 strings. You may need to assign the proj4 string to one or other of the objects.

Spacedman
  • 92,590
  • 12
  • 140
  • 224