我有多组点数据(不同的年份约20组)
我想使用R空间包为每组点生成泰森多边形。
我知道这可以在GIS中完成,但由于我需要批量处理,所以用R完成会更方便。
您没有给我们访问您的数据的权限,但是这里有一个关于代表世界城市的点的例子,使用了Carson Farmer在他的博客中描述的一种方法。希望这能帮助您入门...
# 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)
与jbaums所示的原理相同,但代码更简单:
library(dismo)
library(rgdal)
cities <- shapefile(file.path(system.file("vectors", package = "rgdal")[1], "cities"))
v <- voronoi(cities)
plot(v)
sf::st_voronoi()
函数。以下是从此帮助页面启发的示例: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包(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)
2020年5月27日使用reprex package (v0.3.0)创建
install.packages("sos"); library("sos");
,搜索"thiessen"函数:findFn("thiessen")
。 - Ben Bolker