maptools/SpatialPolygons2PolySet produces ascending instead of descending POS numbers for polygons with holes

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

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