Voronoi approach to buffering polygons while preserving topological integrity

Question

As I understand it R lacks a methods to buffer polygons in a spatially exclusive way that preserves the topology of adjacent polygons. So I'm experimenting with an approach that generates voronoi polygons of the original polygon vertices. Results seem quite promising except for apparent errors in the voronoi generation.

Fairly old school R, so it's possible a tidier alternative may work better. This reproducible example uses US/Canada, but note the problem is one of mathematical geometry so marine boundaries are not relevant:

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)

Let's say we want to bin the grey points to nearest polygon. I think the most elegant approach would be to create a new expanded set of polygons. This avoids lots of n-squared nearest neighbour calculations. Next we try a voronoi tesselation of the original polygon vertices:

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)])

Lots of errors in here. Maybe due to different polygons sharing some vertices. Let's try small negative buffer to help the algorithm:

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)])

Some improvements - almost validating the approach I think. But again we still have some errors, e.g. blue chunk of British Colombia and a thin pink strip of easter border area in Alaska. Lastly I plot with a bigger buffer to help show what is happening with individual vertices (click for bigger resolution):

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)])

Is anyone able to shed light on the problem, or possible suggest an alternative voronoi method that works? I've tried ggvoronoi but struggled to get that working. Any assistance appreciated.

@robert-hijmans If not already seen you might be interested in this new package https://twitter.com/pjs_228/status/1307940206054694912 — geotheory, Sep 21 '20 at 09:28

score 2 · Accepted Answer · edited Jul 06 '20 at 13:08

That is an interesting, and important, problem; and I think it is a good idea to use voronoi. The apparent errors arise from the distribution of the vertices. For example, the border between Canada and the USA hardly has vertices in the west. This leads to undesired results, but they are not wrong. A step in the right direction might be to add vertices, using 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"]])

Dissolve the polygons by country and remove holes

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)

Now use this to cut the buffer by country

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")

And now for something more focused. The below does essentially the same as the above, but focuses just on the regions that matter (coastal borders) making it computationally less intensive --- although not so much for this example with a rather large buffer.

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"))

Really helpful, thanks Robert. I take your point about not being errors. Do you know why the slivers you identify come about? They seem to exist only inside the original polygons. Running `spatialEco::remove.holes(v)` does a good job of removing these particular slivers. Is this something you'd recommend as a general solution? — geotheory, Jul 05 '20 at 08:49
Great. I have expanded the code to take if further. spatialEco::remove.holes is fine, but I show an alternative that might be good enough for this case.. — Robert Hijmans, Jul 05 '20 at 18:31
Thanks again. For me now the challenge here is how to decide geosphere::makePoly's interval parameter. Too big and you see the artifacts I first found. Too small and you hit your compute limit at the voronoi step. There's not necessarily a sweet spot! I've renamed the question to be more useful. Feel free to edit. — geotheory, Jul 05 '20 at 20:42
The new method is a big improvement computationally. Also if this does fully avoid the sliver holes (?) then you don't have a problem of how to handle enclaved polygons. I've made suggestion for a more general solution of _n_ polygons. What do you think? — geotheory, Jul 06 '20 at 13:04
Hmm looks like a bug in the original version (pre my edit): set `bsz <- 1` and various areas of buffered zone gets misclassified. — geotheory, Jul 06 '20 at 14:30
I've made a couple of suggested changes in [this version](https://stackoverflow.com/a/62773258/1156245). What's your view? — geotheory, Jul 07 '20 at 10:34

score 1 · Answer 2 · answered Jul 07 '20 at 10:31

Slight adaptation from Robert's, for discussion.

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)

Alas testing on `w <- w[w$ADMIN %in% c('United States of America', 'Canada', 'Mexico'),]` fails at voronoi with deldir error `Cannot find an enclosing triangle`. Any idea what causes this? — geotheory, Jul 07 '20 at 10:48

Voronoi approach to buffering polygons while preserving topological integrity

2 Answers2