Combine Voronoi polygons and maps
Asked Answered
C

1

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

Cannonry answered 28/8, 2012 at 9:34 Comment(5)
Its not clear what you want. You have the voronoi polygons as a spatialpolygonsdataframe in the pzn.coords object.Hertahertberg
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') .Loco
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.Cannonry
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.Hertahertberg
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, JedSharpeared
H
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.

Hertahertberg answered 28/8, 2012 at 13:2 Comment(2)
You have a stray pznall on line 7 (should be poly)Utility
I concur, rw = as.numeric(t(bbox(pznall))) should be replaced with rw = as.numeric(t(bb)) (as bbox(poly) was already saved to bb).Flyer

© 2022 - 2024 — McMap. All rights reserved.