maptools/SpatialPolygons2PolySet produces ascending instead of descending POS numbers for polygons with holes
Thank you for the suggestion. Unfortunately, this does not work for my application: I plot the polygon on top of other another plot where I need the holes of the polygon to be transparent. For instance, I don?t want the African countries in the following examples to be covered by white: library(maps) map() plot(polySP, density=10, angle=45, pbg=?white") However, this was not the point that I try to ask about. My point is that maptools/SpatialPolygons2PolySet does not translate into correct PolySet objects as illustrated in my previous example. Thanks, Daniel
On Aug 6, 2014, at 12:52 PM, Roger Bivand <Roger.Bivand at nhh.no> wrote:
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