Search code examples
rbufferterra

Negative buffer in Terra has unexpected results for complex multi-polygons


I have been playing around with the buffer function in terra this week. I was applying a negative buffer to get an area within my study area. My initial input polygons created multiple odd results and crashed R several times. I found the following solutions but I suspect that they are work-arounds and there might be better ways to solve them.

  • Issue no. 1 : aggregated polygons create an odd polygon that is outside the original polygons

  • Solution : run disagg on original file

  • Issue no. 2 : polygons with holes (or inner rings) have a buffer that is created outside the holes/rings which can go beyond the outer polygon

  • Solution : run fillHoles on original file

  • Issue no. 3 : R crashed very often (although irregularly)

  • Solution : use unprojected coordinate system instead of original projected crs of input file

I updated terra last week and am running version 1.7-55 on a Windows 11. I have a 13th Gen Intel(R) Core(TM) i7-1355U 1.70 GHz and 16GB RAM.

library(terra)

issue 1 : aggregated polygons create an odd polygon that is outside the original polygons

#three-polygon spatVector which generates odd inset buffer
p1<-vect("POLYGON ((-66.940374 45.008097, -66.971762 45, -66.986546 45.008447, -66.940267 45.026992, -66.940374 45.008097))")                                                                                                                                                                                                              
p2<-vect("POLYGON ((-66.828 45.032592, -66.829 45.035, -66.824 45.031, -66.8301 45.0275, -66.834 45.027, -66.828 45.032592))")                                                                                                                                                                                 
p3<-vect("POLYGON ((-67 46, -66.2 46, -66.2 45.65, -66.37 45.45, -66.37 45.105, -66.459 45.05, -66.46 45.17,
-66.78 45.04, -66.76 45.12, -66.89 45.04, -66.916 45.11, -66.89 45.1, -66.82 45.127, -66.96 45.19,
-67 45.15, -67 46))")

wkt.p<-rbind(p1, p2, p3)
crs(wkt.p)<-"+proj=longlat +datum=WGS84"

hol_buff<-buffer(aggregate(wkt.p), -2000)
plot(wkt.p, ylim=c(45, 46))
plot(hol_buff, add=T, col="green")

#dissaggreated polygons do not create odd inset
hol_buff<-buffer(wkt.p, -2000)
plot(hol_buff, add=T, col="red")

#three-polygon spatVector which does NOT generates odd inset buffer (upper limit of p3 changed from 46 to 45.8)
p1<-vect("POLYGON ((-66.940374 45.008097, -66.971762 45, -66.986546 45.008447, -66.940267 45.026992, -66.940374 45.008097))")                                                                                                                                                                                                              
p2<-vect("POLYGON ((-66.828 45.032592, -66.829 45.035, -66.824 45.031, -66.8301 45.0275, -66.834 45.027, -66.828 45.032592))")                                                                                                                                                                                 
p3<-vect("POLYGON ((-67 45.8, -66.2 45.8, -66.2 45.65, -66.37 45.45, -66.37 45.105, -66.459 45.05, -66.46 45.17,
-66.78 45.04, -66.76 45.12, -66.89 45.04, -66.916 45.11, -66.89 45.1, -66.82 45.127, -66.96 45.19,
-67 45.15, -67 45.8))")

wkt.p<-rbind(p1, p2, p3)
crs(wkt.p)<-"+proj=longlat +datum=WGS84"

hol_buff<-buffer(aggregate(wkt.p), -2000)
plot(wkt.p, ylim=c(45, 46))
plot(hol_buff, add=T, col="blue")

enter image description here

#issue no 2 : holes within polygons leads to negative buffer that are outside the polygon boundary


par(mfrow=c(1,2))
#polygon with inner ring
wkt.y<-"POLYGON ((-66.902 45.013, -66.901 45.01, -66.906 45.0114,  -66.9045 45.0155, -66.902 45.0155, -66.902 45.013),
(-66.903371 45.012734, -66.903372 45.012914, -66.903626 45.012913, -66.903624 45.012733, -66.903371 45.012734))"
y<-vect(wkt.y,crs="+proj=longlat +datum=WGS84")
plot(y, main="terra output")
plot(buffer(y, -40), add=T, col="red")

#polygon without inner ring
wkt.w<-"POLYGON ((-66.902 45.013, -66.901 45.01, -66.906 45.0114,  -66.9045 45.0155, -66.902 45.0155, -66.902 45.013))"
w<-vect(wkt.w,crs="+proj=longlat +datum=WGS84")
plot(y, main="desired output")
plot(buffer(w, -40), add=T, col="red")#negative buffer using the simple polygon that shows the desired buffered polygon
plot(y, add=T)#hypothetical input
#removing the inner ring can be done by using fillHoles
w2<-fillHoles(y)


#polygon with inner ring at edge
par(mfrow=c(1,2))
wkt.y<-"POLYGON ((-66.902 45.012,  -66.906 45.0114,  -66.906 45.0155, -66.902 45.0155, -66.902 45.012),
(-66.903371 45.012734, -66.903372 45.012914, -66.903626 45.012913, -66.903624 45.012733, -66.903371 45.012734))"
y<-vect(wkt.y,crs="+proj=longlat +datum=WGS84")
#plot(buffer(y, -200), col="red", )
plot(y, xlim=c(-66.907, -66.901), ylim=c(45.01, 45.017), main="terra output")
plot(buffer(y, -200), add=T, col="red", border="forestgreen", lwd=2)
plot(y, add=T)#overlay original polygon

#polygon without inner ring
wkt.w<-"POLYGON ((-66.902 45.012,  -66.906 45.0114,  -66.906 45.0155, -66.902 45.0155, -66.902 45.012))"
w<-vect(wkt.w,crs="+proj=longlat +datum=WGS84")
plot(y, main="desired output")
plot(buffer(w, -200), add=T, col="red")#negative buffer using the simple polygon that shows the desired buffered polygon
plot(y, add=T)#hypothetical input

enter image description here enter image description here


Solution

  • A much better place to report issues like this is the terra github site.

    If you do not want the holes to grow with a negative buffer then you could remove them first, as you do. But it would be possible to add alternative options, such as shrinking the holes.

    Otherwise, I think that the issues you have reported have been fixed in terra version 1.7-63. That is currently the development version. You can install it with

    install.packages('terra', repos='https://rspatial.r-universe.dev')