[R-sig-Geo] gUnaryUnion Not Dissolving Correctly
Matt Strimas-Mackey
strimas at zoology.ubc.ca
Wed Nov 4 16:38:41 CET 2015
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().
Thanks!
M
On Wed, Nov 4, 2015 at 6:30 AM, Michael Sumner <mdsumner at gmail.com> wrote:
>
>
> 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
>>
>
More information about the R-sig-Geo
mailing list