[R-sig-Geo] gUnaryUnion Not Dissolving Correctly

Roger Bivand Roger.Bivand at nhh.no
Wed Nov 4 17:57:22 CET 2015 

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.

Roger

>
> 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
>>>
>>
>

-- 
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