[R-sig-Geo] maptools/SpatialPolygons2PolySet produces ascending instead of descending POS numbers for polygons with holes

[Roger Bivand]

On Wed, 6 Aug 2014, Daniel Rodolphe Schlaepfer wrote:

> Hello all,
>
> I try to convert a sp SpatialPolygons object with a hole to a PBS 
> PolySet object - my ultimate goal is to plot the polygon with hatching 
> considering the hole.

The problem is not that you need a PolySet representation, but that you 
need to set the polygon background explicitly to a value other than 
"transparent" using the pbg= argument:

plot(polySP, density=10, angle=45, pbg="white")

or adjust par("bg") to suit. When hatching is used, polypath is not used, 
so automatic handling of holes in the plot method is not available.

Roger

>
> My problem is that maptools/SpatialPolygons2PolySet only produces 
> PolySet objects with increasing POS even for polygons with holes; 
> however, the documentation of PolySet indicates that “We adopt the 
> convention that POS goes from 1 to n along an outer boundary, but from n 
> to 1 along an inner boundary, regardless of rotational direction.”
>
>
> #Create simple doughnut-shaped polygon
> library(sp)
> library(maptools)
>
> coords1 <- matrix(c(108, -54, -108, -54, -108, 54, 108, 54, 108, -54), ncol=2, byrow=TRUE)
> coords2 <- matrix(c(36, -18, -36, -18, -36, 18, 36, 18, 36, -18), ncol=2, byrow=TRUE)
>
> polySP <- SpatialPolygons(list(Polygons(list(Polygon(coords1, hole=FALSE), Polygon(coords2, hole=TRUE)), ID=1)), proj4string=CRS("+proj=longlat +datum=WGS84"))
>
> #Convert sp-SpatialPolygons to PBS-PolySet with maptools function
> polyPBS <- SpatialPolygons2PolySet(polySP)
>
> #-> POS for SID == 2 are increasing and thus do not reflect PBS standards for a polygon with a hole
>   PID SID POS    X   Y
> 1    1   1   1  108 -54
> 2    1   1   2 -108 -54
> 3    1   1   3 -108  54
> 4    1   1   4  108  54
> 5    1   1   5  108 -54
> 6    1   2   1   36 -18
> 7    1   2   2   36  18
> 8    1   2   3  -36  18
> 9    1   2   4  -36 -18
> 10   1   2   5   36 -18
>
> I believe that it would require only a small change to maptools/SpatialPolygons2PolySet to produce PolySet objects that meet the PBS standards also for polygons with holes, i.e., replace the line
> 	POS <- c(POS, 1:k)
> with
> 	POS <- if(slot(srs[[j]], "hole")) c(POS, k:1) else c(POS, 1:k)
> inside the loops: for (i in 1:n) … for (j in 1:m) …
>
>
>
> Sincerely,
> Daniel Schlaepfer
>
>
> My session infos:
> sessionInfo()
> 	R version 3.1.1 (2014-07-10)
> 	Platform: x86_64-apple-darwin13.3.0 (64-bit)
>
> 	locale:
> 	[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
>
> 	attached base packages:
> 	[1] stats     graphics  grDevices
> 	[4] utils     datasets  methods
> 	[7] base
>
> 	other attached packages:
> 	[1] PBSmapping_2.67.60 maptools_0.8-30
> 	[3] sp_1.0-15
>
> 	loaded via a namespace (and not attached):
> 	[1] foreign_0.8-61  grid_3.1.1
> 	[3] lattice_0.20-29
>
>
> -------------------------------------------------------
> Daniel Schlaepfer, PhD
> University of Wyoming
> Laramie, WY 82071
>
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo at 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.Bivand at nhh.no