将泰森多边形和地图结合起来

Question

我想将 Voronoi 多边形与地图结合起来，以便稍后用于空间分析。我有一些点和形状文件要组合，然后保存为形状文件/空间多边形。为了获得 Voronoi 多边形，我使用此主题中的函数。

以下是我的代码：

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

请注意，此链接需要复制并粘贴到浏览器中才能打开。

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'))))
}

我的代码获取Voronoi多边形：

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

结果：

结果：

我希望在我的地图中将这个沃罗诺伊多边形作为新的空间多边形。

如果能得到答案，我将不胜感激！

编辑：

明确一下，我想实现类似于这样的效果（这些线条应该是从沃罗诺伊多边形创建的）：

- Maciej

---

你想要什么并不清楚。在pzn.coords对象中，你有voronoi多边形作为一个spatialpolygonsdataframe。 - Spacedman

您提供的链接是一个 .rdata 文件。您可以使用 save() 创建一个新的 myvoronietc.rdata 文件，例如 save('pzn.coords','lotsa_other_data',file='myvoronietc.rdata')。 - Carl Witthoft

我想将这两个文件合并，使我的地图内有 Voronoi 多边形，而不是在我的地图边界上有正方形边界的 Voronoi 多边形。 - Maciej

你需要将Voronoi多边形扩展到比你的多边形更大，可能通过向deldir添加虚拟点来实现，然后使用rgeos进行裁剪。我之后可能会提供完整的答案，或者你可以自行尝试。 - Spacedman

只是一条注释（不是答案）：这行代码中有一个拼写错误： rw = as.numeric(t(bbox(pznall)))，请将pznall替换为poly。感谢您的修改， 祝好， Jed - jed.a.long

Answer 1

稍作修改的函数，接受一个额外的空间多边形参数并将其扩展到该框：

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 字符串分配给其中一个对象。

- Spacedman

---

第 7 行出现一个无法识别的 pznall（应该是 poly）。 - jbaums

我同意，rw = as.numeric(t(bbox(pznall))) 应该被替换为 rw = as.numeric(t(bb))（因为 bbox(poly) 已经保存到了 bb）。 - Tamas Ferenci