问题描述
我想将Voronoi多边形与地图结合起来,以便以后将其用于空间分析.我有一些要合并的点和shapefile,然后另存为shapefile/空间多边形.要获取voronoi多边形,我可以使用此主题.
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.
我的代码如下:
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))
我的地图在这里可用: https://docs.google.com /file/d/0B-ZJyVlQBsqlSURiN284dF9YNUk/edit
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'))))
}
我的代码得到voronoipolygons:
And my code to get voronoipolygons:
pzn.coords<-voronoipolygons(coords)
plot(pznall)
plot(pzn.coords,add=T)
points(coords$LONG,coords$LAT)
结果:
我想在地图中将此voronoi多边形作为新的spacepolygon.
I wan to have this voronoi polygon inside my map as new spatialpolygon.
我将感谢安瑟斯!
需要明确的是,我想实现以下目标(此行应从voronoi多边形创建):
推荐答案
稍作修改的函数,采用附加的空间多边形参数并扩展到该框:
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)
}
然后做:
pzn.coords<-voronoipolygons(coords,pznall)
library(rgeos)
gg = gIntersection(pznall,pzn.coords,byid=TRUE)
plot(gg)
请注意,gg
是SpatialPolygons对象,您可能会收到有关proj4字符串不匹配的警告.您可能需要将proj4字符串分配给一个或其他对象.
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.
这篇关于合并Voronoi多边形和地图的文章就介绍到这了,希望我们推荐的答案对大家有所帮助,也希望大家多多支持!