On Wed, 6 Aug 2014, Daniel Rodolphe Schlaepfer wrote:
Unfortunately, your suggestion to consider plotting order doesn’t work
in general and not in my case; for instance, if I need my underlaying
map to show filled polygons:
plot(polySP, density=10, angle=45, pbg="white", axes=TRUE)
map(add=TRUE, fill=TRUE, col=rainbow(n=5))
-the hatched polygon is overplotted by the map with this plotting order;
if I reverse the plotting order, then the underlaying map is
overplotted:
map(fill=TRUE, col=rainbow(n=5))
plot(polySP, density=10, angle=45, pbg="white", add=TRUE)
Thus, I still think this issue is not resolved.
Given the graphics devices at out disposal, one has to operate within the
posibilities available. By the way, polypath does not support hatching, so
the plot method in sp is as good as it gets. Had you considered using the
alpha channel to add colour fills over a hatched background?
I can understand that fixing maptools/SpatialPolygons2PolySet may have a
low priority. However, I already suggested a patch in my original email.
Here, again, as the complete function (only one line needs to be
changed, see second instance of POS <- ):
Yes, but from there no indication of how PBSmapping::plotMap() gets around
the same problem.
I've committed your suggestion to R-forge, maptools revision 283. Please
report whether this is what you wanted.
Roger
SpatialPolygons2PolySet <- function (SpP)
{
require(PBSmapping)
pls <- slot(SpP, "polygons")
n <- length(pls)
PID <- NULL
SID <- NULL
POS <- NULL
X <- NULL
Y <- NULL
for (i in 1:n) {
srs <- slot(pls[[i]], "Polygons")
m <- length(srs)
for (j in 1:m) {
crds <- slot(srs[[j]], "coords")
k <- nrow(crds)
PID <- c(PID, rep(i, k))
SID <- c(SID, rep(j, k))
#POS <- c(POS, 1:k)
POS <- if(slot(srs[[j]], "hole")) c(POS, k:1) else c(POS, 1:k)
#Suggested fix
X <- c(X, crds[, 1])
Y <- c(Y, crds[, 2])
}
}
PID <- as.integer(PID)
SID <- as.integer(SID)
POS <- as.integer(POS)
storage.mode(X) <- "double"
storage.mode(Y) <- "double"
pj <- .pbsproj(SpP)
zn <- NULL
if (pj == "UTM") {
zn <- attr(pj, "zone")
attr(pj, "zone") <- NULL
}
res <- as.PolySet(data.frame(PID = PID, SID = SID, POS = POS,
X = X, Y = Y), projection = pj, zone = zn)
res
}
Thanks,
Daniel
On Aug 6, 2014, at 3:02 PM, Roger Bivand <roger.biv...@nhh.no> wrote:
On Wed, 6 Aug 2014, Daniel Rodolphe Schlaepfer wrote:
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")
Obviously not, these are base graphics, and you have to consider the plotting
order:
plot(polySP, density=10, angle=45, pbg="white", axes=TRUE)
map(add=TRUE)
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.
Since you no longer need PolySet objects, why? Maybe the coercion method can be
fixed, but it is not a priority - patch welcomed.
Roger
Thanks,
Daniel
On Aug 6, 2014, at 12:52 PM, Roger Bivand <roger.biv...@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@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.biv...@nhh.no
--
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.biv...@nhh.no
--
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.biv...@nhh.no
_______________________________________________
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
