[R-sig-Geo] gUnaryUnion Not Dissolving Correctly

Matt Strimas-Mackey strimas at zoology.ubc.ca
Wed Nov 4 19:10:37 CET 2015 

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.

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:

library(sp)
library(raster)
library(rgeos)

# just a subset of full shapefile
set.seed(1)
size <- sqrt(2 * 1e8 / sqrt(3))
study_area <- shapefile('papua.shp')
hex_points <- spsample(study_area[2,], type="hexagonal", cellsize=size)
hex_grid <- HexPoints2SpatialPolygons(hex_points)
hex_union <- gUnaryUnion(hex_grid)
plot(hex_grid, col='lightgrey')
plot(hex_union, border='orange', lwd=3, add=T)

# results chage according to seed
set.seed(100)
size <- sqrt(2 * 1e8 / sqrt(3))
study_area <- shapefile('papua.shp')
hex_points <- spsample(study_area[2,], type="hexagonal", cellsize=size)
hex_grid <- HexPoints2SpatialPolygons(hex_points)
hex_union <- gUnaryUnion(hex_grid)
plot(hex_grid, col='lightgrey')
plot(hex_union, border='orange', lwd=3, add=T)

On Wed, Nov 4, 2015 at 8:57 AM, Roger Bivand <Roger.Bivand at nhh.no> wrote:
> 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