r - 包含在地理边界中的 Voronoi 图多边形

Question

我正在尝试在固定地理区域内为一组点创建 Voronoi 多边形(又名 Dirichlet 镶嵌或泰森多边形)。但是，我无法在 R 中找到一种将多边形绑定(bind)在 map 边界内的方法。我的主要目标是获得准确的面积计算(而不仅仅是生成视觉图)。例如，以下内容直观地传达了我想要实现的目标:

library(maps)
library(deldir)
data(countyMapEnv)
counties <- map('county', c('maryland,carroll','maryland,frederick', 'maryland,montgomery', 'maryland,howard'), interior=FALSE)
x <- c(-77.208703, -77.456582, -77.090600,  -77.035668, -77.197144)
y <- c(39.188603, 39.347019, 39.672818, 39.501898, 39.389203)
points(x,y)
vt <- deldir(x, y, rw=counties$range)
plot(vt, wlines="tess", lty="solid", add=TRUE)

产生以下内容:



从概念上讲，我想将 counties 与 vt 相交，它应该提供一组以县边界为界的多边形，并为每个多边形提供准确的面积计算。现在， vt$summary 为每个多边形提供了面积计算，但是对于除了一个内部多边形之外的所有多边形，它们显然被夸大了，并且 deldir() 似乎只接受其 rw 参数的矩形包围。我对 R 的地理空间功能不熟悉，因此对我上面概述的其他方法持开放态度。

Answer 1

您应该能够为此使用 spatstat 函数 dirichlet。

第一个任务是将县转换为 owin 类的 spatstat 观察窗口(代码部分基于@jbaums 的答案):

library(maps)
library(maptools)
library(spatstat)
library(rgeos)

counties <- map('county', c('maryland,carroll', 'maryland,frederick', 
                            'maryland,montgomery', 'maryland,howard'), 
                fill=TRUE, plot=FALSE)
# fill=TRUE is necessary for converting this map object to SpatialPolygons
countries <- gUnaryUnion(map2SpatialPolygons(counties, IDs=counties$names,
                                 proj4string=CRS("+proj=longlat +datum=WGS84")))
W <- as(countries, "owin")

然后你只需将五个点放入 ppp 格式，进行 Dirichlet tesslation 并计算面积:

X <- ppp(x=c(-77.208703, -77.456582, -77.090600,  -77.035668, -77.197144),
         y=c(39.188603, 39.347019, 39.672818, 39.501898, 39.389203), window = W)

y <- dirichlet(X) # Dirichlet tesselation
plot(y) # Plot tesselation
plot(X, add = TRUE) # Add points
tile.areas(y) #Areas