通过 Voronoi 方法对多边形进行缓冲，同时保持拓扑完整性

Question

通过 Voronoi 方法对多边形进行缓冲，同时保持拓扑完整性

3

据我所知，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但是遇到了困难。任何帮助都将不胜感激。

- geotheory

如果你还没看过的话，你可能会对这个新软件包感兴趣 https://twitter.com/pjs_228/status/1307940206054694912 - geotheory

2个回答

1

轻微修改自Robert的代码，用于讨论。

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

w <- rworldmap::countriesCoarse[,'ADMIN']
# w <- w[w$ADMIN %in% c('United States of America', 'Canada'),]
w <- w[w$ADMIN %in% c('Guyana', 'Suriname','French Guiana'),]
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 <- .5

# outward buffer
mbuf = buffer(p, width = bsz, dissolve=F)

# overlay between two country buffers
# e <- mbuf[1,] * mbuf[2,]
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]]

# contested buffer zones - overlap minus original polys
f <- e - p
f@data = data.frame(id = seq_along(f))

# buffer the contested zones
b <- buffer(f, bsz)

# bp is the area that matters
bp <- b * p

# vertices
bp = buffer(bp, width = -0.00002, dissolve=F)
g0 <- data.frame(data.frame(geom(bp)))
# getting rid of duplicated and shared vertices
# g <- aggregate(g0[,'object', drop=FALSE], g0[,c('x','y')], min)
g = unique(g0)

v0 <- dismo::voronoi(g[,c("x", "y")], extend(extent(p), 2 * bsz))
v0$id = g$object
v <- raster::aggregate(v0, list(g[,"object"]), FUN=length)
v@proj4string = p@proj4string
v = v * f
v@data = data.frame(ADMIN = p$ADMIN[v$Group.1])

# full buffer
fb = raster::bind(mbuf - p - f, v, p)
fb = raster::aggregate(fb, list(fb$ADMIN), FUN = function(x)x[1])[,'ADMIN']
fb@proj4string = p@proj4string

#----------------------------------

par(mai=c(0,0,0,0))
plot(p, border='grey')
plot(mbuf, add=T, border='pink')
plot(e, add=T, col='#00000010', border=NA)
plot(f, add=T, border='purple', lwd=1.5)
plot(b, add=T, border='red')
plot(bp, add=T, col='#ffff0040', border=NA)
# plot(v, add=T, col=c("#ff770020", "#0077ff20"), border=c("#ff7700", "#0077ff"))
plot(fb, add=T, col=c("#ff000020", "#00ff0020", "#0000ff20"), border=NA)

- geotheory

哎呀，对 w <- w[w$ADMIN %in% c('United States of America', 'Canada', 'Mexico'),] 进行测试时，在 Voronoi 中出现了 deldir 错误 Cannot find an enclosing triangle。你有什么想法是什么原因导致的吗？ - geotheory

网页内容由stack overflow 提供, 点击上面的

可以查看英文原文，
原文链接

- Robert Hijmans · Accepted Answer

这是一个有趣且重要的问题；我认为使用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"))