Importance of Scale and Slivers in RGEOS

rgeos-scale-slivers.r

library(sp)
library(raster)
library(rgeos)
set_RGEOS_polyThreshold(0)
set_RGEOS_warnSlivers(TRUE)
set_RGEOS_dropSlivers(FALSE)

# function to rigidly shift polygon
# only works for simple polygon objects with a single polygon
shift_poly <- function(p, x, y) {
  p@polygons[[1]]@Polygons[[1]]@coords[,'x'] <- 
    p@polygons[[1]]@Polygons[[1]]@coords[,'x'] + x
  p@polygons[[1]]@Polygons[[1]]@coords[,'y'] <- 
    p@polygons[[1]]@Polygons[[1]]@coords[,'y'] + y
  row.names(p) <- paste0('shifted', row.names(p))
  return(p)
}

p1 <- readWKT("POLYGON((0 0,0 1,1 1,1 0,0 0))")
row.names(p1) <- '1'
p2 <- readWKT("POLYGON((1 0,1 1,2 1,2 0,1 0))")
row.names(p2) <- '2'

plot(rbind(p1, p2))
plot(shift_poly(p2, -0.5, 0.1), border='orange', add=T, lty=2, lwd=1)

# default scale of 10^8
# behaviour of topology operations depends on scale!
setScale(1e8)
# shift horizontally by small amount so no longer touching
plot(gUnion(p1, shift_poly(p2, 1e-1, 0)))
plot(gUnion(p1, shift_poly(p2, 1e-8, 0)))
plot(gUnion(p1, shift_poly(p2, 1e-9, 0)))

gIntersects(p1, p2)
gIntersects(p1, shift_poly(p2, 1e-1, 0))
gIntersects(p1, shift_poly(p2, 1e-8, 0))
gIntersects(p1, shift_poly(p2, 1e-9, 0))
# polygons with coordinates different by less that the scale set by setScale()
# are considered to intersect and are merge together by union operations

# p1 and p2 share a boundary exactly so the intersection of the 2 is a line
class(gIntersection(p1, p2))
# shift right polygon horizontally slightly to it is just overlapping or just
# separated from the left polygon and consider the results
gIntersection(p1, shift_poly(p2, -1e-9, 0)) # overlap > scale => polygon
gIntersection(p1, shift_poly(p2, -1e-8, 0)) # overlap < scale => line
gIntersection(p1, p2) # exactly touching => line
gIntersection(p1, shift_poly(p2, 1e-9, 0)) # separation < scale => line
gIntersection(p1, shift_poly(p2, 1e-8, 0)) # separation > scale => NULL

# lower scale to 10^4
setScale(1e4)
plot(gUnion(p1, shift_poly(p2, 1e-1, 0)))
plot(gUnion(p1, shift_poly(p2, 1e-4, 0)))
plot(gUnion(p1, shift_poly(p2, 1e-5, 0)))

gIntersects(p1, p2)
gIntersects(p1, shift_poly(p2, 1e-1, 0))
gIntersects(p1, shift_poly(p2, 1e-4, 0))
gIntersects(p1, shift_poly(p2, 1e-5, 0))

gIntersection(p1, shift_poly(p2, -1e-4, 0)) # overlap > scale => polygon
gIntersection(p1, shift_poly(p2, -1e-5, 0)) # overlap < scale => line
gIntersection(p1, p2) # exactly touching => line
gIntersection(p1, shift_poly(p2, 1e-5, 0)) # separation < scale => line
gIntersection(p1, shift_poly(p2, 1e-4, 0)) # separation > scale => NULL


# consider the effect of setting polyThreshold and turning on sliver warning
# this will identify slivers resulting from topology operations
setScale(1e4)
set_RGEOS_polyThreshold(1e-2)
set_RGEOS_warnSlivers(TRUE)
# shift horizontally by increasing amount so just overlapping
gi <- gIntersection(p1, shift_poly(p2, -1e-5, 0))
class(gi) # overlap < scale => line, i.e. assumes just touching
# warnings raised in next 2 cases because:
#   a. overlap > scale => polygon on intersection
#   b. resulting polygon area < polyThreshold => sliver
gIntersection(p1, shift_poly(p2, -1e-4, 0)) 
gIntersection(p1, shift_poly(p2, -1e-3, 0))
# warnings raised in next 2 cases because:
#   a. overlap > scale => polygon on intersection
#   b. resulting polygon area > polyThreshold => not a sliver
gIntersection(p1, shift_poly(p2, -1e-2, 0))
gIntersection(p1, shift_poly(p2, -1e-1, 0))

# with a lower threshold
set_RGEOS_polyThreshold(1e-3)
# this still raises a warning
gIntersection(p1, shift_poly(p2, -1e-4, 0)) 
# but this doesn't since resulting polygon is at the threshold now
gIntersection(p1, shift_poly(p2, -1e-3, 0))

# not that it isn't the linear overlap that triggers the warning, it is that the
# area of the resulting polygons are below the threshold
# no warning
gi1 <- gIntersection(p1, shift_poly(p2, -1e-3, 0))
# now a warning is raised because a slight shift in the y direction has caused
# the polygons the resulting polygon to be just less than the 1e-3 threshold
gi2 <- gIntersection(p1, shift_poly(p2, -1e-3, 1e-3))
gArea(gi1) / get_RGEOS_polyThreshold()
gArea(gi2) / get_RGEOS_polyThreshold()

# rgeos can also be set to automatically remove these slivers
class(gIntersection(p1, shift_poly(p2, -1e-4, 0))) # SpaitalPolygons
set_RGEOS_dropSlivers(TRUE)
class(gIntersection(p1, shift_poly(p2, -1e-4, 0))) # NULL
set_RGEOS_dropSlivers(FALSE)


# slivers can also be generated as a result of union operations
setScale(1e4)
set_RGEOS_polyThreshold(1e-2)
set_RGEOS_warnSlivers(TRUE)
set_RGEOS_dropSlivers(FALSE)

p3 <- readWKT("POLYGON((0 1,0 2,2 2,2 1,0 1))")
row.names(p3) <- '3'
p4 <- readWKT("POLYGON((0 -1,0 0,2 0,2 -1,0 -1))")
row.names(p4) <- '4'
plot(rbind(p1, p3, p4))
plot(p2, add=T, col='red')

# offset the middle right (i.e. red) square slightly to the right
# not picked up since middle edge misalignment is < scale
pp <- rbind(p1, p3, p4, shift_poly(p2, 1e-5, 0))
plot(gUnaryUnion(pp))
# misalignment picked up since = scale => a very narrow hole in center
# warning raised because hole area is < polyThreshold
pp <- rbind(p1, p3, p4, shift_poly(p2, 1e-4, 0))
plot(gUnaryUnion(pp))
# misalignment picked up since = scale => a very narrow hole in center
# warning not raised because hole area is > polyThreshold
pp <- rbind(p1, p3, p4, shift_poly(p2, 1e-2, 0))
plot(gUnaryUnion(pp))
# the fact that this is a hole and not a vertical line becomes apparent when
# the shift is bigger
pp <- rbind(p1, p3, p4, shift_poly(p2, 1e-1, 0))
plot(gUnaryUnion(pp))
# finally, set_RGEOS_dropSlivers can be used to remove these slivers 
# and fix the topology
set_RGEOS_dropSlivers(TRUE)
pp <- rbind(p1, p3, p4, shift_poly(p2, 1e-4, 0))
plot(gUnaryUnion(pp))
set_RGEOS_dropSlivers(FALSE)