[R-sig-Geo] Holes in polygons disappear after performing checkPolygonsHoles from maptools
Grigory Alexandrovich
alexandrovich at mathematik.uni-marburg.de
Fri Apr 3 22:19:51 CEST 2015
Am 03.04.2015 um 20:42 schrieb Roger Bivand:
> On Fri, 3 Apr 2015, Grigory Alexandrovich wrote:
>
>>
>>
>> Am 03.04.2015 um 13:16 schrieb Roger Bivand:
>>> On Fri, 3 Apr 2015, Grigory Alexandrovich wrote:
>>>
>>> > Am 02.04.2015 um 22:24 schrieb Roger Bivand:
>>> > > On Thu, 2 Apr 2015, Roger Bivand wrote:
>>> > > > > > On Thu, 2 Apr 2015, Edzer Pebesma wrote:
>>> > > > > > > > > > > > > > > > On 04/01/2015 08:03 PM, Grigory
>>> Alexandrovich wrote:
>>> > > > > > Hello all,
>>> > > > > > > I encountered the following unexpected result using
>>> > > > > checkPolygonsHoles:
>>> > > > > > > # attach the worldmap as SpatialPolygonsDataFrame from
>>> the
>>> > > > > package > maptools
>>> > > > > > library(sp)
>>> > > > > > library(maptools)
>>> > > > > > data(wrld_simpl)
>>> > > > > > > # get a polygon with a hole
>>> > > > > > shape_with_hole <- wrld_simpl[5,]
>>> > > > > > > # plot it (hole is left white, surrounded by blue color)
>>> > > > > > plot(shape_with_hole, col = "blue")
>>> > > > > > > # perform checkPolygonsHoles
>>> > > > > > shape_with_hole at polygons <-
>>> lapply(shape_with_hole at polygons,
>>> > > > > > checkPolygonsHoles)
>>> > > > > > > # plot again, now holes aren't recognized as such
>>> > > > > > plot(shape_with_hole, col = "blue")
>>> > > > > > > # and even the original SpatialPolygonsDataFrame
>>> object is
>>> > > > > changed !?
>>> > > > > > plot(wrld_simpl[5,], col = "blue")
>>> > > > > > > One irritating side effect here is that the original
>>> object
>>> > > > > wrld_simpl
>>> > > > > > is also changed. This result looks to me like a bug, or
>>> have I
>>> > > > > missed
>>> > > > > > something?
>>> > > > > > > > > > It looks as if it is related to rgeos, as
>>> > > > > > > > > avoidGEOS = FALSE
>>> > > > > library(sp)
>>> > > > > library(maptools)
>>> > > > > data(wrld_simpl)
>>> > > > > x = wrld_simpl
>>> > > > > > > > > shape_with_hole <- wrld_simpl[5,]
>>> > > > > > > > > if (avoidGEOS)
>>> > > > > gpclibPermit()
>>> > > > > shape_with_hole at polygons <- lapply(shape_with_hole at polygons,
>>> > > > > checkPolygonsHoles,
>>> > > > > avoidGEOS = avoidGEOS)
>>> > > > > > > > > data(wrld_simpl)
>>> > > > > all.equal(x, wrld_simpl)
>>> > > > > > > No, but:
>>> > > > > > > library(sp)
>>> > > > library(maptools)
>>> > > > data(wrld_simpl)
>>> > > > x = wrld_simpl
>>> > > > shape_with_hole <- wrld_simpl[5,]
>>> > > > pls <- slot(shape_with_hole, "polygons")
>>> > > > pls1 <- lapply(pls, checkPolygonsHoles, avoidGEOS = FALSE)
>>> > > > isTRUE(all.equal(pls, pls1))
>>> > > > # TRUE
>>> > > > gpclibPermit()
>>> > > > pls1 <- lapply(pls, checkPolygonsHoles, avoidGEOS = TRUE)
>>> > > > isTRUE(all.equal(pls, pls1))
>>> > > > # FALSE
>>> > > > > > > where the Polygon objects are re-ordered and plotted
>>> out of order, > > > so
>>> > > > the deprecated gpclib code may be broken here. However, this
>>> isn't
>>> > > > the problem,
>>> > > > > The gpclib code is OK, it re-orders the Polygon objects,
>>> but the plot
>>> > > order is correct. The problem is as described below.
>>> > > > > Roger
>>> > > > > > but rather (perhaps) graphics::polypath. If it is used
>>> (by default > > > it
>>> > > > is), and with the default path fill mode:
>>> > > > > > > plot(x[5,], col="blue", usePolypath=TRUE, rule="winding")
>>> > > > > > > does not whiten the holes. However:
>>> > > > > > > plot(x[5,], col="blue", usePolypath=TRUE, rule="evenodd")
>>> > > > > > > does. Using:
>>> > > > > > > plot(x[5,], col="blue", pbg="white", usePolypath=FALSE)
>>> > > > > > > also works by overpainting the island with the holes in
>>> white.
>>> > > > > > > Please summarise the your SO thread.
>>> > > > > > > By the way, in R it is unusual for changes in one
>>> object to > > > propagate
>>> > > > in an unintended way to another object, and this is not
>>> happening
>>> > > > here. The visual impression is being driven by polypath doing
>>> > > > different things depending on its rule setting. Look at
>>> ?polypath > > > for
>>> > > > details.
>>> > > > > > > Well, this doesn't clarify, but it explains what is
>>> happeing, I hope > > > ...
>>> > > > > > > Roger
>>> > > > > > > > > > > > returns a difference, but TRUE if avoidGEOS =
>>> TRUE.
>>> > > > > > > > > > | Thank you for your answers guys. But I am still
>>> not fully convinced
>>> > that it is an issue of the plot function,
>>> > because the object shape_with_hole edited with checkPolygonsHoles
>>> > before, continues to behave in a strange way:
>>> > > # we check which polygons are flagged as holes.
>>> > # The flags are still set properly,
>>> > # although the plot() function didn't recognize them:
>>> > sapply(shape_with_hole at polygons[[1]]@Polygons, slot, "hole")
>>> > > [1] FALSE TRUE TRUE TRUE
>>> > > # load library rgdal for reprojection
>>> > library(rgdal)
>>> > > # reproject with spTransform, just for testing
>>> > shape_with_hole <- spTransform(shape_with_hole, CRS("+proj=longlat
>>> > +ellps=WGS84 +datum=WGS84"))
>>> > > # after reprojection all flags are set to FALSE
>>> > sapply(shape_with_hole at polygons[[1]]@Polygons, slot, "hole")
>>> > > [1] FALSE FALSE FALSE FALSE
>>> > > > So it is still something wrong with shape_with_hole, or?
>>>
>>> Please provide the output of sessionInfo(), rgeos::version_GEOS(), and
>>> rgdal::getGDALVersionInfo() with rgdal::getPROJ4VersionInfo(). I
>>> cannot
>>> reproduce the degradation of hole status following a no-op use of
>>> spTransform() with:
>>>
>>> > sessionInfo()
>>> R version 3.1.3 (2015-03-09)
>>> Platform: x86_64-unknown-linux-gnu (64-bit)
>>> Running under: Fedora 21 (Twenty One)
>>> ...
>>> attached base packages:
>>> [1] stats graphics grDevices utils datasets methods base
>>>
>>> other attached packages:
>>> [1] rgdal_0.9-3 maptools_0.8-34 sp_1.0-17
>>>
>>> loaded via a namespace (and not attached):
>>> [1] foreign_0.8-63 grid_3.1.3 lattice_0.20-31 rgeos_0.3-8
>>> [5] tools_3.1.3
>>> > rgeos::version_GEOS()
>>> [1] "3.4.2-CAPI-1.8.2 r3921"
>>> > rgdal::getGDALVersionInfo()
>>> [1] "GDAL 1.11.2, released 2015/02/10"
>>> > rgdal::getPROJ4VersionInfo()
>>> [1] "Rel. 4.9.1, 04 March 2015, [PJ_VERSION: 491]"
>>>
>>> Roger
>>>
>>>
>>> > > Best
>>> > > Grigory|
>>> > >
>> My session infos:
>>
>> sessionInfo()
>> R version 3.1.1 (2014-07-10)
>> Platform: x86_64-w64-mingw32/x64 (64-bit)
>>
>> ...
>> attached base packages:
>> [1] stats graphics grDevices utils datasets methods base
>>
>> other attached packages:
>> [1] rgdal_0.9-1 maptools_0.8-34 sp_1.0-17
>>
>> loaded via a namespace (and not attached):
>> [1] foreign_0.8-61 grid_3.1.1 lattice_0.20-29 tools_3.1.1
>>
>>
>> rgeos::version_GEOS()
>> [1] "3.4.2-CAPI-1.8.2 r3921"
>>
>> rgdal::getGDALVersionInfo()
>> [1] "GDAL 1.11.0, released 2014/04/16"
>>
>> rgdal::getPROJ4VersionInfo()
>> [1] "Rel. 4.8.0, 6 March 2012, [PJ_VERSION: 480]"
>>
>>
>> But I also obtained the same result on another windows machine with R
>> 3.1.2 (dont have session info for it)
>
> How often do you run update.packages() ? This can be done any time in
> the 3.1 release cycle, as packages are compatible with R releases at
> that level unless otherwise stated. rgdal_0.9-1 was built linking
> against sp_1.0-15, as we don't lock-step packages depending on sp with
> its releases. So some of the code for creating sp objects being used
> in rgeos and rgdal at the C level is possibly out of date. Could you
> please try update.packages() with regard to the last problem which I
> cannot reproduce?
>
> Roger
>
>>
>> Grigory
>>
>>
>>
>
Roger, actually I didn't run update.packages() on the corresponding
spatial packages yet, since I installed them maybe two-three months ago.
Anyway, I just installed the newest R Version (3.1.3) and the required
packages and the problem doesn't appear any more.
So one should take care using older versions of the rgeos/rgdal packages.
Thank you and Edzer for the concurrence and your contributions to R's
spatial analysis toolkit.
Grigory
More information about the R-sig-Geo
mailing list