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

Grigory Alexandrovich alexandrovich at mathematik.uni-marburg.de
Fri Apr 3 14:48:59 CEST 2015



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)

Grigory



More information about the R-sig-Geo mailing list