[R-sig-Geo] maptools/SpatialPolygons2PolySet produces ascending instead of descending POS numbers for polygons with holes
Daniel Rodolphe Schlaepfer
dschlaep at uwyo.edu
Thu Aug 7 09:52:39 CEST 2014
Since polypath does not support hatching, I had to resort to using PBSmapping. However, I don’t know how PBSmapping is getting around the problem of correctly plotting of hatching for polygons with holes, but it does it correctly and that is why I use it, e.g.,
if(packageVersion(“maptools”) >= “0.8.31"){
polyPBS <- SpatialPolygons2PolySet(polySP)
map(fill=TRUE, col=rainbow(n=5))
addPolys(polyPBS, xlim=c(-180, 180), ylim=c(-90, 90), density=10, angle=45, col="black")
}
Thank you for including the fix for SpatialPolygons2PolySet to maptools revision 283. My question is solved.
Best
-Daniel
On Aug 6, 2014, at 8:54 PM, Roger Bivand <Roger.Bivand at nhh.no> wrote:
> 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.Bivand at 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.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
>>>>
>>>>
>>>
>>> --
>>> 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
>>
>>
>
> --
> 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
More information about the R-sig-Geo
mailing list