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. 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 <- ):
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
_______________________________________________
R-sig-Geo mailing list
R-sig-Geo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-geo
