[R-sig-Geo] gUnaryUnion Not Dissolving Correctly

Wed Nov 4 15:30:42 CET 2015

On Thu, 5 Nov 2015 at 01:10 Roger Bivand <Roger.Bivand at nhh.no> wrote:

> 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

 points(coordinates(as(as(spoly, "SpatialLines"), "SpatialPoints"))[c(3,
4), ])

Thanks for the reminder about cleangeo, I'll have closer look.

cheers, Mike.

> Roger
>
> >
> > Cheers, Mike.
> >
> > On Wed, 4 Nov 2015 at 18:16 Roger Bivand <Roger.Bivand at nhh.no> wrote:
> >
> >> On Tue, 3 Nov 2015, Matt Strimas-Mackey wrote:
> >>
> >>> 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.
> >>>
> >>> Example:
> >>>
> >>> ## Create three adjacent hexagons
> >>> library(sp)
> >>> library(rgeos)
> >>> p1 <- Polygon(cbind(
> >>>      c(1276503.26781119, 1281876.11747031, 1287248.96712942,
> >>>        1287248.96712942, 1281876.11747031, 1276503.26781119,
> >> 1276503.26781119),
> >>>      c(204391.40834643, 207493.42454344, 204391.40834643,
> >> 198187.37595242,
> >>>        195085.35975541, 198187.37595242, 204391.40834643)))
> >>> p2 <- Polygon(cbind(
> >>>      c(1287248.96712943, 1292621.81678854, 1297994.66644766,
> >>>        1297994.66644766, 1292621.81678854, 1287248.96712943,
> >> 1287248.96712943),
> >>>      c(204391.40834643, 207493.42454344, 204391.40834643,
> >> 198187.37595242,
> >>>        195085.35975541, 198187.37595242, 204391.40834643)))
> >>> p3 <- Polygon(cbind(
> >>>      c(1281876.11747031, 1287248.96712943, 1292621.81678854,
> >>>        1292621.81678854, 1287248.96712943, 1281876.11747031,
> >> 1281876.11747031),
> >>>      c(213697.45693745, 216799.47313446, 213697.45693745,
> >> 207493.42454344,
> >>>        204391.40834643, 207493.42454344, 213697.45693745)))
> >>> spoly <- SpatialPolygons(list(Polygons(list(p1, p2, p3), 's1')))
> >>> plot(gUnaryUnion(spoly))
> >>
> >> No, this is just numerical fuzz:
> >>
> >> plot(spoly)
> >> plot(gUnaryUnion(spoly), border="green", lty=3, lwd=2, add=TRUE)
> >> oS <- getScale()
> >> # default 1e+8
> >> setScale(1e+4)
> >> plot(gUnaryUnion(spoly), border="orange", lwd=2, add=TRUE)
> >> setScale(oS)
> >>
> >> 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.
> >>
> >> Hope this clarifies,
> >>
> >> Roger
> >>
> >>>
> >>> Note that p2 and p3 are dissolved together, but p1 is separate. The
> >>> shared edge of p1 and p2 is:
> >>> p1:
> >>> [2,] 1281876 207493.4
> >>> [3,] 1287249 204391.4
> >>> p2:
> >>> [5,] 1287249 204391.4
> >>> [6,] 1281876 207493.4
> >>>
> >>> 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 at polygons, function(x) slot(x, 'Polygons')[[1]]@coords)
> >>> plot(sss)
> >>> plot(gUnaryUnion(sss))
> >>>
> >>> Session Info:
> >>> R version 3.2.2 (2015-08-14)
> >>> Platform: x86_64-apple-darwin13.4.0 (64-bit)
> >>> Running under: OS X 10.10.5 (Yosemite)
> >>>
> >>> Message when rgeos is loaded:
> >>>
> >>> rgeos version: 0.3-14, (SVN revision 511)
> >>> GEOS runtime version: 3.5.0-CAPI-1.9.0 r0
> >>> Linking to sp version: 1.2-0
> >>> Polygon checking: TRUE
> >>>
> >>> Any help on how to get these polygons to dissolve is appreciated.
> >>>
> >>> M
> >>>
> >>> _______________________________________________
> >>> R-sig-Geo mailing list
> >>> R-sig-Geo at r-project.org
> >>> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
> >>>
> >>
> >> --
> >> 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
> >>
> >> _______________________________________________
> >> R-sig-Geo mailing list
> >> R-sig-Geo at r-project.org
> >> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
> >>
> >
>
> --
> 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
>
>

	[[alternative HTML version deleted]]