[R-sig-Geo] gUnaryUnion Not Dissolving Correctly

[Matt Strimas-Mackey]

I've rewritten the script as an RMarkdown document and posted it on github:
https://github.com/mstrimas/rgeos-scale
The rendered html is here:
https://htmlpreview.github.io/?https://github.com/mstrimas/rgeos-scale/blob/master/rgeos-scale.html

If you or Colin Rundel want to modify this or use it in any way feel
free. It would be great if it helped clarify some of these issues for
others.

M

On Fri, Nov 6, 2015 at 6:29 AM, Roger Bivand <Roger.Bivand at nhh.no> wrote:
> 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
>