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

Daniel Rodolphe Schlaepfer dschlaep at uwyo.edu
Wed Aug 6 15:44:44 CEST 2014 

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.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

More information about the R-sig-Geo mailing list