如何使用 R 包从点创建泰森多边形？答案

【问题标题】：How to create thiessen polygons from points using R packages?如何使用 R 包从点创建泰森多边形？
【发布时间】：2012-03-13 07:17:54
【问题描述】：

我有多组积分（不同年份~20）

我想使用 r 个空间包为每组点生成泰森多边形。

我知道这可以使用 GIS 来完成，但是因为我想要一个批处理，所以 R 中的东西会是

有帮助。

【问题讨论】：

install.packages("sos"); library("sos"); findFn("thiessen")
我现在也在使用 ArcGIS ...

标签： r geometry polygon spatial

【解决方案1】：

您尚未允许我们访问您的数据，但这里有一个代表世界城市的点的示例，使用了 Carson Farmer 在his blog 上描述的方法。希望它能让你开始......

# Carson's Voronoi polygons function
voronoipolygons <- function(x) {
  require(deldir)
  require(sp)
  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))
  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'))))
}

示例 1：输入为 SpatialPointsDataFrame：

# Read in a point shapefile to be converted to a Voronoi diagram
library(rgdal)
dsn <- system.file("vectors", package = "rgdal")[1]
cities <- readOGR(dsn=dsn, layer="cities")

v <- voronoipolygons(cities)

plot(v)

示例 2：输入是 x、y 坐标的向量：

dat <- data.frame(x=runif(100), y=runif(100))
v2 <- voronoipolygons(dat)
plot(v2)

【讨论】：

我已经修改了函数，使其接受坐标向量（按 x、y 的顺序）以及 SpatialPointsDataFrame。
有什么办法可以为不规则多边形边界做到这一点？你可以在这里找到我的问题stackoverflow.com/questions/29746150/…
谢谢，这正是我要找的！ ;)

【解决方案2】：

原理与jbaums相同，但代码更简单：

library(dismo)
library(rgdal)
cities <- shapefile(file.path(system.file("vectors", package = "rgdal")[1], "cities"))

v <- voronoi(cities)
plot(v)

【讨论】：

【解决方案3】：

注意Voronoi cells are also known as Thiessen polygons。

或者，可以使用Simple Features for R，具体来说，是sf::st_voronoi() 函数。以下是受this help page启发的示例：

library(sf)
#> Linking to GEOS 3.8.0, GDAL 3.0.4, PROJ 6.3.1

# generate some random points
set.seed(2020-05-27)
n <- 50
points <- runif(n) %>% 
  matrix(ncol = 2) %>% 
  st_multipoint()

# Voronoi tesselation
voronoi_grid <- st_voronoi(points)

plot(voronoi_grid, col = NA)
plot(points, add = TRUE, col = "blue", pch = 16)

^{由reprex package (v0.3.0) 于 2020 年 5 月 27 日创建}

由于您提到您有多组点，一组一年，您可以使用列表并对其进行迭代，例如：

library(sf)
#> Linking to GEOS 3.8.0, GDAL 3.0.4, PROJ 6.3.1

# generate a list of length 20, each element containing with random points 
set.seed(2020-05-27)
yrs <- 20
n <- 50
points_lst <- vector(mode = "list", length = yrs)
for (i in 1:yrs) {
  points_lst[[i]] <- runif(n) %>% 
    matrix(ncol = 2) %>% 
    st_multipoint()
}

# Voronoi tesselation for each element of the list
voronoi_grids_lst <- lapply(points_lst, st_voronoi)

# plot 1st element
plot(voronoi_grids_lst[[1]], col = NA)

^{由reprex package (v0.3.0) 于 2020 年 5 月 27 日创建}

【讨论】：