[R-sig-Geo] Hole attribute in Polygon lost when made into Polygons object?

[Michael Sumner]

On Wed, 27 May 2015 at 08:26 MacQueen, Don <macqueen1 at llnl.gov> wrote:

> I must be missing something basic
> (either in understanding, or possibly being blind to a typo!)
>
> Here is the example data:
>
> m1 <-structure(c(2.8, 8, 8.2, 3.9, 1.6,
>     9.2, 6.8, 3, 1.1, 4.2),
>    .Dim = c(5L, 2L))
>
>
> Pnh <- Polygon(m1, hole=FALSE)
> Ph  <-  Polygon(m1, hole=TRUE)
>
> Pnhs <- Polygons(list(Pnh), ID=1)
> Phs <- Polygons(list(Ph), ID=1)
>
>
>
> Then:
> > Pnh at hole
> [1] FALSE
> > Ph at hole
> [1] TRUE
>
> > Pnhs at Polygons[[1]]@hole
> [1] FALSE
> > Phs at Polygons[[1]]@hole
> [1] FALSE
>
> It appears that the hole attribute has been lost??
>
>
>


I believe that the hole attribute only makes sense when it's relative to at
least one other Polygon object. If you combine your Pnh and Ph in the same
Polygons object, the hole status is preserved (despite this being a broken
object, see below):

 sapply(Polygons(list(Pnh, Ph), ID = 1)@Polygons, slot, "hole")
[1] FALSE  TRUE

If you create a Polygons object with just one Polygon that is a 'hole',
there's nothing for it to be a hole of (there's no 'background'  in sp).  I
presume this happens internally when a new Polygons is initialized with
just one ring, it makes no sense for the hole attribute  to be TRUE (the sp
source is here src/sp_xports.c but I'm not au fait with that).

Note that both of your Polygon objects are identical, so you get a nonsense
object:


rgeos::gIsValid(SpatialPolygons(list(Polygons(list(Pnh, Ph), ID = 1))))
[1] FALSE
Warning message:
In RGEOSUnaryPredFunc(spgeom, byid, "rgeos_isvalid") :
  Self-intersection at or near point 1.6000000000000001 4.2000000000000002

Also, you can see the different ways the topology is treated by creating a
hole, and then offset that to incorrectly overlap
it's background. With usePolypath the topology of the situation is more
clearly represented
x1 <- Polygon(m1, hole=FALSE)
x2  <-  Polygon((m1 + 3) / 1.5, hole=TRUE)
p1 <- SpatialPolygons(list(Polygons(list(x1, x2), ID = 1)))

op <- par(mfrow = c(2, 2))
plot(p1, col = "grey", usePolypath = TRUE, main = "ok, polypath TRUE")
plot(p1, col = "grey", usePolypath = FALSE, main = "ok, polypath FALSE")

## now with a hole, but incorrectly overlapping

y1 <- Polygon(m1, hole=FALSE)
y2  <-  Polygon((m1 + 4) / 1.5, hole=TRUE)
p2 <- SpatialPolygons(list(Polygons(list(y1, y2), ID = 1)))

plot(p2, col = "grey", usePolypath = TRUE, main = "broken, polypath TRUE")
plot(p2, col = "grey", usePolypath = FALSE, main = "broken, polypath FALSE")
par(op)

Does that help?

Cheers, Mike.



> sp_1.1-0
> R 3.1.2
>
> Thanks
> -Don
>
> --
> Don MacQueen
>
> Lawrence Livermore National Laboratory
> 7000 East Ave., L-627
> Livermore, CA 94550
> 925-423-1062
>
> _______________________________________________
> R-sig-Geo mailing list
> R-sig-Geo at r-project.org
> https://stat.ethz.ch/mailman/listinfo/r-sig-geo
>

	[[alternative HTML version deleted]]