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.biv...@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 <mdsum...@gmail.com>
> > wrote:
> > >
> > >
> > >
> > > On Thu, 5 Nov 2015 at 01:10 Roger Bivand <roger.biv...@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.biv...@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@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@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.biv...@nhh.no
> > > > > >
> > > > > > _______________________________________________
> > > > > > R-sig-Geo mailing list
> > > > > > R-sig-Geo@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.biv...@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.biv...@nhh.no
>
