[R-sig-Geo] Holes in polygons disappear after performing checkPolygonsHoles from maptools

Roger Bivand Roger.Bivand at nhh.no
Fri Apr 3 20:42:56 CEST 2015


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