[R-sig-Geo] gUnaryUnion Not Dissolving Correctly

[Roger Bivand]

Matt:

Thanks very much for this. I think that it would make an excellent 
vignette for rgeos, maybe you could consider using Markdown on your 
existing script to explain to the reader what they might expect? I'd also 
bring in Colin Rundel, as it was his insight that uncovered the scale "can 
of worms" five years ago, if you would like to go ahead.

Maybe also use maptools::elide methods instead of shift_poly() - maptools 
is Suggests: in rgeos, so is assumed to be able to be available 
(conditionally) when the package is built, installed and checked.

Actually, this is somewhat like a spatial version of FAQ 7.31:

https://cran.r-project.org/doc/FAQ/R-FAQ.html#Why-doesn_0027t-R-think-these-numbers-are-equal_003f

Best wishes,

Roger

On Fri, 6 Nov 2015, Matt Strimas-Mackey wrote:

> Thanks for explaining this in such detail! I have a greater
> appreciation for the importance of thinking about these topological
> issues and for the role machine precision plays.
>
> I constructed a series of simple examples to demonstrate to myself how
> these sorts of problems arise and how setScale,
> set_RGEOS_polyThreshold, set_RGEOS_dropSlivers, and
> set_RGEOS_warnSlivers work. In case someone else stumbles on this
> thread looking for similar clarification, I've pasted these examples
> below, and they also appear in this gist:
> https://gist.github.com/mstrimas/1b4a4b93a9d4a158bce4
>
> 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 at polygons[[1]]@Polygons[[1]]@coords[,'x'] <-
>    p at polygons[[1]]@Polygons[[1]]@coords[,'x'] + x
>  p at polygons[[1]]@Polygons[[1]]@coords[,'y'] <-
>    p at 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)
>

-- 
Roger Bivand
Department of Economics, Norwegian School of Economics,
Helleveien 30, N-5045 Bergen, Norway.
voice: +47 55 95 93 55; fax +47 55 95 91 00
e-mail: Roger.Bivand at nhh.no