gUnaryUnion Not Dissolving Correctly

gUnaryUnion Not Dissolving Correctly

I'm working with a regular hexagonal grid stored as SPDF. At some
point I subset this SPDF, then want to combine all adjacent hexagons
together so that each contiguous set of hexagons is a single polygon.
I'm doing this last step using gUnaryUnion (or gUnionCascaded, not
clear what the different is actually). The problem is that in some
cases boundaries between clearly adjacent polygons are not dissolved.

So, exactly the same apart from the order. I originally thought this
difference in order might be the problem, but this doesn't seem to be
an issue with in this example, where order is also flipped:
sss <- rasterToPolygons(raster(nrow=2, ncol=2, xmn=0, xmx=2, ymn=0,
ymx=2, vals=1:4))
lapply(sss@polygons, function(x) slot(x, 'Polygons')[[1]]@coords)
plot(sss)
plot(gUnaryUnion(sss))

JTS, GEOS, and consequently rgeos by default shift all coordinates to an
integer grid after multiplying by a scale factor (finding integer matches
is much easier than real matches). If the scaling is too detailed (in some
cases), the operations do not give the expected outcomes.

There is work in progress in GEOS and JTS to provide other scaling options
and models, and to permit iteration over scaling values until a "clean"
result is obtained (for some meanings of clean).

gUnionCascaded() was the only possible function for GEOS < 3.3.0, from
GEOS 3.3.0 gUnaryUnion() is available and the prefered and more efficient
route. This is explained on the help page.

Re: gUnaryUnion Not Dissolving Correctly

Thanks for all this detail Roger, is there a way to "re-build" a spatial
object so that the given scale setting is applied? Are there any general
rounding or "orthogonalize" functions in the Spatial suite?

Re: gUnaryUnion Not Dissolving Correctly

Administrator

On Wed, 4 Nov 2015, Michael Sumner wrote:

> Thanks for all this detail Roger, is there a way to "re-build" a spatial
> object so that the given scale setting is applied? Are there any general
> rounding or "orthogonalize" functions in the Spatial suite?

No, not really. In this case, the very detailed coordinate measurements
may have made things worse, or possibly using Polygon rather than Polygons
objects, or not building the Polygons object with a proper comment
attribute, I don't know. rgeos::gIsValid(spoly) is FALSE, and

says the same (based on rgeos). I suggest working with Emmanuel Blondel
(cleangeo maintainer) to extend cleangeo. GEOS is looking at allowing
users to manipulate precision models, not just scale, but I'm uncertain
about that.

Running spoly into GRASS and back out (GRASS builds topology on import)
shows a different error, the object seems to be problematic.

Re: gUnaryUnion Not Dissolving Correctly

> On Wed, 4 Nov 2015, Michael Sumner wrote:
>
> > Thanks for all this detail Roger, is there a way to "re-build" a spatial
> > object so that the given scale setting is applied? Are there any general
> > rounding or "orthogonalize" functions in the Spatial suite?
>
> No, not really. In this case, the very detailed coordinate measurements
> may have made things worse, or possibly using Polygon rather than Polygons
> objects, or not building the Polygons object with a proper comment
> attribute, I don't know. rgeos::gIsValid(spoly) is FALSE, and
>
> cleangeo::clgeo_GeometryReport(clgeo_CollectionReport(clgeo_Clean(spoly)))
>
> says the same (based on rgeos). I suggest working with Emmanuel Blondel
> (cleangeo maintainer) to extend cleangeo. GEOS is looking at allowing
> users to manipulate precision models, not just scale, but I'm uncertain
> about that.
>
> Running spoly into GRASS and back out (GRASS builds topology on import)
> shows a different error, the object seems to be problematic.
>
>

It can be fixed by changing "1287248.96712942" to "1287248.96712943", so
really the creator should not have had those values on input. There's no
easy way out without topology.

This brought out some interesting issues for work I've been doing using the
"near-Delaunay triangulation" in RTriangle, and that requires a normalized
set of vertices (no duplicate vertices) which on its own presents
interesting problems. I have a related issue where a parallel (latitude)
line needs many vertices to look "smooth" on a polar projection, but when
building a mesh with triangles it's really best to allow relatively coarse
segmented boundaries rather than have many elements at parallels. The
Triangle library does not consider these hexagon coordinates to be
duplicates, so there are two vertical segments between the two bottom polys
at

Re: gUnaryUnion Not Dissolving Correctly

Thanks, lots of useful info here. I've never seen the setScale()
function; I don't think it's mentioned in the gUnaryUnion help. This
saves me a lot of headache!

For what it's worth, the invalid geometry is an artifact of the
reproducible example I created. The original hexagonal grid is
produced with
g <- spsample(study_area, type="hexagonal", cellsize=size)
hex_grid <- HexPoints2SpatialPolygons(g)

And this object passes gIsValid() and clgeo_GeometryReport() without
any problems, yet still has the dissolving issue. Regardless, all is
solved with setScale().

Re: gUnaryUnion Not Dissolving Correctly

Administrator

On Wed, 4 Nov 2015, Matt Strimas-Mackey wrote:

> Thanks, lots of useful info here. I've never seen the setScale()
> function; I don't think it's mentioned in the gUnaryUnion help. This
> saves me a lot of headache!
>
> For what it's worth, the invalid geometry is an artifact of the
> reproducible example I created. The original hexagonal grid is
> produced with
> g <- spsample(study_area, type="hexagonal", cellsize=size)
> hex_grid <- HexPoints2SpatialPolygons(g)

OK, thanks - this is useful. Could you please make available the study
area object in some way - so that we can re-create g and see how
HexPoints2SpatialPolygons() creates the artefact Mike spotted (although
this looks like numeric fuzz - 'changing "1287248.96712942" to
"1287248.96712943"' is a change in the 15th digit, which is on the
precision edge of the "double" storage mode. If we can revisit functions
creating SpatialPolygons objects to ensure that they are GEOS-compatible
for the default scale of 1e+8, we'll be more secure.

Since it's a large area, generating the grid takes a long time, so
I've also included code for a small subset of the original
shapefile--one small offshore island.

Finally, some more odd behaviour. I noticed each time I run this code,
the dissolve mistakes change, i.e. different boundaries are
erroneously kept. However, using set.seed() makes the errors the same
each time for a given seed, and changing the seed yields a different
set of errors. Example in the code in the dropbox link and copied
here:

Re: gUnaryUnion Not Dissolving Correctly

Administrator

On Wed, 4 Nov 2015, Matt Strimas-Mackey wrote:

> The original study area shapefile is a boundary of the Indonesia half
> of New Guinea. The file as well as the code to construct the hexagonal
> grids are here:
> https://www.dropbox.com/sh/ff8v08p3ambqcbs/AAAPBlGP4fthdmZhrto7oIuCa?dl=0>
> Since it's a large area, generating the grid takes a long time, so
> I've also included code for a small subset of the original
> shapefile--one small offshore island.

An equivalent scaling fix is to change the coordinate units from metres to
kilometres:

Re: gUnaryUnion Not Dissolving Correctly

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

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

# 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()

Re: gUnaryUnion Not Dissolving Correctly

Administrator

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.