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 但很难让它发挥作用。任何帮助表示赞赏。

4

2 回答 2

2

这是一个有趣且重要的问题。我认为使用 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"))
于 2020-07-05T02:49:18.510 回答
1

稍微改编自罗伯特的，供讨论。

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)

在此处输入图像描述

于 2020-07-07T10:31:54.343 回答