看到一段生成在R语言中生成泰森多边形的代码,感觉很有用。这段代码的最初来源是Carson Farmer的个人网站,网址是:https://carsonfarmer.com/2009/09/voronoi-polygons-with-r/。
voronoipolygons = function(layer) {
require(deldir)
require(sp)
crds = layer@coords
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'))))
}
这个函数只有一个参数layer
,用起来很方便。不过它是基于sp结构的对象运行的,如果使用的矢量数据是sf结构的,需要进行相关转换。但是直接使用as
函数进行转换会报错:
library(sf)
library(sp)
usa <- albersusa::counties_sf(proj = "laea")
usa.sp <- as(usa, "Spatial")
v <- voronoipolygons(usa.sp) # error
经过摸索,发现下面的代码是可行的:
library(sf)
library(sp)
usa <- albersusa::counties_sf(proj = "laea")
# to sp
coord <- st_coordinates(st_centroid(usa))
usa.data <- st_drop_geometry(usa)
usa.sp <- SpatialPointsDataFrame(coords = coord, data = usa.data)
# 生成泰森多边形
v <- voronoipolygons(usa.sp)
# 可视化
plot(v)
网址https://www.it1352.com/790552.html也提供了两个使用示例。
往期推荐阅读: