gUnaryUnion Not Dissolving Correctly

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:

library(rgdal)
study_area2 <- spTransform(study_area, CRS("+proj=aea +lat_1=-7.42
 +lat_2=-0.62 +lat_0=-9.119 +lon_0=128.91 +x_0=0 y_0=0 +datum=WGS84
 +units=km +no_defs +ellps=WGS84 +towgs84=0,0,0"))
size <- sqrt(2 * 1e2 / sqrt(3)) # reduce to suit
set.seed(1)
hex_points <- spsample(study_area2[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)

An advantage of these scaling issues being made visible is that we see that 
computers really do not operate in continuous space, and that computational 
geometry actually matters, I suppose.

I'll check the underlying code generating the hexagons to see why the nodes 
(where hexagon boundaries meet) appear to slide apart at the edge of machine 
precision.

Roger

 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

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

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