在保持拓扑完整性的同时缓冲多边形的 Voronoi 方法
Voronoi approach to buffering polygons while preserving topological integrity
据我所知，R 缺少一种方法来以空间独占的方式缓冲多边形，以保留相邻多边形的拓扑结构。所以我正在试验一种生成原始多边形顶点的 voronoi 多边形的方法。除了 voronoi 生成中的明显错误外，结果似乎很有希望。
相当老派的 R，所以更整洁的替代方案可能会更好。此可重现示例使用 US/Canada，但请注意问题是数学几何问题之一，因此海洋边界不相关：
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 但很难让它工作。感谢任何帮助。
这是一个有趣且重要的问题；我认为使用 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"))
罗伯特的轻微改编，供讨论。
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)
据我所知，R 缺少一种方法来以空间独占的方式缓冲多边形，以保留相邻多边形的拓扑结构。所以我正在试验一种生成原始多边形顶点的 voronoi 多边形的方法。除了 voronoi 生成中的明显错误外，结果似乎很有希望。
相当老派的 R，所以更整洁的替代方案可能会更好。此可重现示例使用 US/Canada，但请注意问题是数学几何问题之一，因此海洋边界不相关：
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 但很难让它工作。感谢任何帮助。
这是一个有趣且重要的问题；我认为使用 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"))
罗伯特的轻微改编，供讨论。
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)