Skip to content

Instantly share code, notes, and snippets.

# mstrimas/rgeos-scale-slivers.r Created Nov 6, 2015

Importance of Scale and Slivers in RGEOS
 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[]@Polygons[]@coords[,'x'] <- p@polygons[]@Polygons[]@coords[,'x'] + x p@polygons[]@Polygons[]@coords[,'y'] <- p@polygons[]@Polygons[]@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)
to join this conversation on GitHub. Already have an account? Sign in to comment
You can’t perform that action at this time.