r - 在保持拓扑完整性的同时缓冲多边形的 Voronoi 方法

Question

据我了解，R 缺少一种方法来以空间独占的方式缓冲多边形，以保留相邻多边形的拓扑结构。所以我正在试验一种生成原始多边形顶点的 voronoi 多边形的方法。除了 voronoi 生成中的明显错误外，结果似乎很有希望。

相当老派的 R，因此更整洁的替代方案可能会更好。这个可重现的例子使用美国/加拿大，但请注意问题是数学几何问题之一，因此海洋边界不相关:

require(rworldmap)
require(rgeos)
require(dismo)
require(purrr)
require(dplyr)
par(mai = rep(0,4))

p = rworldmap::countriesCoarse[,'ADMIN']
p = p[p$ADMIN %in% c('United States of America', 'Canada'),]
p$ADMIN = as.character(p$ADMIN)
p = rgeos::gBuffer(p, byid=T, width = 0) # precaution to ensure no badly-formed polygon nonsense

# Not critical to the problem, but consider we have points we want to assign to enclosing or nearest polygon
set.seed(42)
pts = data.frame(x = runif(1000, min = p@bbox[1,1], max = p@bbox[1,2]),
                 y = runif(1000, min = p@bbox[2,1], max = p@bbox[2,2]))
coordinates(pts) = pts
pts@proj4string = p@proj4string

# point in polygon classification.
pts$admin = sp::over(pts, p)$ADMIN
pts$admin = replace(pts$admin, is.na(pts$admin), 'unclass')

plot(p)
plot(pts, pch=16, cex=.4, col = c('red','grey','blue')[factor(pts$admin)], add=T)

假设我们想要将灰点分到最近的多边形。我认为最优雅的方法是创建一组新的扩展多边形。这避免了大量的 n 平方最近邻计算。接下来我们尝试对原始多边形顶点进行 voronoi 分割:

vertices1 = map_df(p@polygons, ~ map2_df(.x@Polygons, rep(.x@ID, length(.x@Polygons)),
                               ~ as.data.frame(..1@coords) %>% `names<-`(c('x','y')) %>% mutate(id = ..2)))
print(head(vertices1))
#>           x        y     id
#> 1 -56.13404 50.68701 Canada
#> 2 -56.79588 49.81231 Canada
#> 3 -56.14311 50.15012 Canada
#> 4 -55.47149 49.93582 Canada
#> 5 -55.82240 49.58713 Canada
#> 6 -54.93514 49.31301 Canada
coordinates(vertices1) = vertices1[,1:2]

# voronois
vor1 = dismo::voronoi(vertices1)

# visualise
plot(p)
plot(vertices1, add=T, pch=16, cex=.5, col = c('red','blue')[factor(vertices1$id)])
plot(vor1, add=T, border='#00000010', col = c('#FF000040','#0000FF40')[factor(vor1$id)])

这里有很多错误。可能是由于不同的多边形共享一些顶点。让我们尝试使用小的负缓冲区来帮助算法:

p_buff2 = rgeos::gBuffer(p, byid=T, width = -.00002) # order of 1 metre

vertices2 = map_df(p_buff2@polygons, ~ map2_df(.x@Polygons, rep(.x@ID, length(.x@Polygons)), 
                                     ~ as.data.frame(..1@coords) %>% `names<-`(c('x','y')) %>% mutate(id = ..2)))
coordinates(vertices2) = vertices2[,1:2]

vor2 = dismo::voronoi(vertices2)

plot(p_buff2)
plot(vertices2, add=T, pch=16, cex=.4, col = c('red','blue')[factor(vertices2$id)])
plot(vor2, add=T, border='#00000010', col = c('#FF000040','#0000FF40')[factor(vor2$id)])

一些改进 - 几乎验证了我认为的方法。但是我们仍然有一些错误，例如不列颠哥伦比亚的蓝色大块和阿拉斯加复活节边境地区的粉红色细带。最后，我使用更大的缓冲区进行绘图，以帮助显示各个顶点发生的情况(单击以获得更大的分辨率):

p_buff3 = rgeos::gBuffer(p, byid=T, width = -.5, ) # order of 30kms I think

vertices3 = map_df(p_buff3@polygons, ~ map2_df(.x@Polygons, rep(.x@ID, length(.x@Polygons)), 
                                     ~ as.data.frame(..1@coords) %>% `names<-`(c('x','y')) %>% mutate(id = ..2)))
coordinates(vertices3) = vertices3[,1:2]

vor3 = dismo::voronoi(vertices3)

plot(p_buff3)
plot(vertices3, add=T, pch=16, cex=.4, col = c('red','blue')[factor(vertices3$id)])
plot(vor3, add=T, border='#00000010', col = c('#FF000040','#0000FF40')[factor(vor3$id)])

有没有人能够阐明这个问题，或者可能建议一种可行的替代 voronoi 方法？我试过 ggvoronoi 但很难让它工作。任何帮助表示赞赏。

Answer 1

这是一个有趣且重要的问题；我认为使用 voronoi 是个好主意。明显的错误来自顶点的分布。例如，加拿大和美国之间的边界几乎没有西边的顶点。这会导致不希望的结果，但它们并没有错。朝正确方向迈出的一步可能是使用 geosphere::makePoly 添加顶点

library(dismo)
library(geosphere)
library(rworldmap)
library(rgeos)

w <- rworldmap::countriesCoarse[,'ADMIN']
w <- w[w$ADMIN %in% c('United States of America', 'Canada'),]
p <- geosphere::makePoly(w, 25000)
p$ADMIN = as.character(p$ADMIN)

p <- buffer(p, width = 0, dissolve=FALSE)
p_buff <- buffer(p, width = -.00002, dissolve=FALSE) # order of 1 metre

g <- geom(p_buff)
g <- unique(g)

vor <- dismo::voronoi(g[,c("x", "y")])

plot(p_buff)
points(g[,c("x", "y")], pch=16, cex=.4, col= c('red','blue')[g[,"object"]])
plot(vor, add=T, border='#00000010', col = c('#FF000040','#0000FF40')[g[,"object"]])

按国家分解多边形并去除孔

v <- aggregate(vor, list(g[,"object"]), FUN=length)   
gg <- data.frame(geom(v))
v <- as(gg[gg$hole==0, ], "SpatialPolygons")

lines(v, col="yellow", lwd=4)

现在用这个按国家削减缓冲区

pp <- buffer(p, width = 10)
buf <- v * (pp - p)   # intersect(v, erase(pp, p))
buf <- SpatialPolygonsDataFrame(buf, data=data.frame(p), match.ID = FALSE)
x <- bind(p, buf)
z <- aggregate(x, "ADMIN")

lines(z, lwd=2, col="dark green")

现在是更专注的事情。下面的内容与上面的内容基本相同，但只关注重要的区域(沿海边界)，使其计算强度较低 --- 尽管对于这个具有相当大缓冲区的示例来说并不是那么多。

library(dismo)
library(rworldmap)
library(rgeos)

w <- rworldmap::countriesCoarse[,'ADMIN']
w <- w[w$ADMIN %in% c('United States of America', 'Canada', 'Mexico'),]
p <- geosphere::makePoly(w, 25000)
p$ADMIN = as.character(p$ADMIN)
p <- buffer(p, width = 0, dissolve=FALSE)
#p <- buffer(p, width = -.00002, dissolve=FALSE) # order of 1 metre

bsz <- 10
mbuf <- buffer(p, width = bsz, dissolve=FALSE)
# e <- mbuf[1,] * mbuf[2,]

# -----------
# general solution for e?

poly_combs = expand.grid(p1 = seq_along(mbuf), p2 = seq_along(mbuf))
poly_combs = poly_combs[poly_combs$p1 < poly_combs$p2,]

# pairwise overlaps
e_pw = plyr::compact(lapply(1:nrow(poly_combs), FUN = function(i){
  pair = poly_combs[i,]
  pairing = suppressWarnings(mbuf[pair$p1,] * mbuf[pair$p2,])
  return(pairing)
}))

e = e_pw[[1]]
for(i in 2:length(e_pw)) e = e + e_pw[[i]]
# -----------

f <- e - p
b <- buffer(f, bsz)
# bp is the area that matters
bp <- b * p

g <- data.frame(geom(bp))
# getting rid of duplicated and shared vertices
g <- aggregate(g[,1,drop=FALSE], g[,5:6], min)  
v <- dismo::voronoi(g[,c("x", "y")], extent(p)+ 2 * bsz)
v <- aggregate(v, list(g[,"object"]), FUN=length)   

v <- v- p
buf1 <- buffer(p, width = bsz, dissolve=TRUE)
v <- v * buf1
v@data <- p@data

plot(v, col=c("red", "blue", "green"))