Polygons in R - HRUG March 2016
Michael Sumner
2016-03-03
What is a polygon?
Wikipedia: a polygon is a plane figure that is bounded by a finite chain of straight line segments closing in a loop to form a closed chain or circuit
- Straight line segments?
OR
- Coordinates connected by a boundary “line”?
Does it matter?
The simplest polygon
Aren’t these all the same thing?
Connected vertices
p <- cbind(x = c(0, 0, 0.75, 1, 0.5, 0.8, 0.69, 0),
y = c(0, 1, 1, 0.8, 0.7, 0.6, 0, 0))
plot(p, xlim = xl, ylim = yl)
lines(p)## x y
## pt1 0.00 0.0
## pt2 0.00 1.0
## pt3 0.75 1.0
## pt4 1.00 0.8
## pt5 0.50 0.7
## pt6 0.80 0.6
## pt7 0.69 0.0
## pt8 0.00 0.0
Note that we visited 0,0 twice, when we first started as “pt1” and at the final destination as “pt8”. What if we left out the last 0,0?
npoints <- nrow(p)
plot(p, xlim = xl, ylim = yl, main = "no final link")
lines(p[seq(npoints - 1), ])
Well it’s obvious, we don’t have a polygon so we just end up with a line.
What if we use polygon() rather than lines()?
plot(p, xlim = xl, ylim = yl, main = "default polygon()")
polygon(p)
What happened? We have a closed loop of line segments.
We are seeing right through this polygon. We have to tell polygon() to fill the space it contains with a colour (which might be white).
plot(p, xlim = xl, ylim = yl, main = "grey-fill")
polygon(p, col = "grey")
What if we have two of these kinds of objects? What goes wrong here?
p1 <- cbind(x = c(0, 0, 0.75, 1, 0.5, 0.8, 0.69, 0),
y = c(0, 1, 1, 0.8, 0.7, 0.6, 0, 0))
p2 <- cbind(x = c(0.2, 0.2, 0.3, 0.5, 0.5, 0.2),
y = c(0.2, 0.4, 0.6, 0.4, 0.2, 0.2))
pp <- rbind(p1, p2)
plot(pp, xlim = xl, ylim = yl, main = "yikes!")
polygon(rbind(p1, p2), col = "grey")
We have defined a path and we’ve told the polygon-drawing-tool to follow that path and fill space with a particular rule. (What?).
One way to fix this is to plot the two parts separately, but we have to cheat since the white “hole” is not a hole - we can no longer see the points on the inner ring.
plot(pp, xlim = xl, ylim = yl, main = "cheat's hole")
polygon(p1, col = "grey")
polygon(p2, col = "white")
The proper way to do this is to separate the parts of the path with the missing “NA” where we don’t want the path to stay connected. But, here polygon doesn’t understand that it was meant to draw it as a hole.
pph <- rbind(p1, NA, p2)
plot(pph, xlim = xl, ylim = yl, main = "polygon with hole, or overlapping polygons?")
polygon(pph, col = "grey")
To really get a hole, we need to use polypath, with the right rule.
pph <- rbind(p1, NA, p2)
plot(pph, xlim = xl, ylim = yl, main = "evenodd - clockwise/clockwise")
polypath(pph, col = "grey", rule = "evenodd")
pph <- rbind(p1, NA, p2[nrow(p2):1, ])
plot(pph, xlim = xl, ylim = yl, main = "winding - clockwise/anti-clockwise")
polypath(pph, col = "grey", rule = "winding")
What have we learnt?
- polygon() and polypath() are the same, until we need holes
- the “evenodd” rule counts how many boundaries are crossed from the inside-out
- the “winding” rule counts clock-wise and many anti-clock-wise turns (or vice versa)
Yes it’s more complicated - not just holes, but generally “self-intersecting” polygon paths.
Multi-polygons
This time the second polygon is an “island”.
p1 <- cbind(x = c(0, 0, 0.75, 1, 0.5, 0.8, 0.69, 0),
y = c(0, 1, 1, 0.8, 0.7, 0.6, 0, 0))
p3 <- cbind(x = c(0.2, 0.2, 0.3, 0.5, 0.5, 0.2) + 0.8,
y = c(0.2, 0.4, 0.6, 0.4, 0.2, 0.2))
pp <- rbind(p1, NA, p3)
plot(pp, main = "two polygons")
polypath(pp, col = "grey")
What if the smaller poly should be another colour?
plot(pp, main = "two polygons, or one multi-polygon?")
polypath(pp, col = c("grey", "firebrick"))
What if we had holes and islands, how does the function apply two different colours? (It can’t).
To get different colours we have to call polypath separately.
plot(pp, main = "polypath() multiple")
polypath(p1, col = "grey")
polypath(p3, col = "firebrick")
The sp package makes this easy, right?
library(sp)
## no, this really is not easy
spoly <- SpatialPolygons(list(Polygons(list(Polygon(p1)), "p1"), Polygons( list(Polygon(p3)), "p3")))
plot(spoly, col = c("grey", "firebrick"), main = "sp SpatialPolygons")
(It’s another level of detail if we want attributes on each object.)
Raster to the rescue, at least a little.
library(raster)
rpoly <- spPolygons(p1, p3)
plot(rpoly, col = c("grey", "firebrick"), main = "raster::spPolygons")
See ?spPolygons for details, you can include a data frame of attributes too.
What if we want that hole again? There’s a “nesting” of polygon rings, but our “hole” has to be in reverse order.
rhpoly <- spPolygons(list(p1, p2[nrow(p2):1, ]), p3)
plot(rhpoly, col = c("grey", "firebrick"), main = "SpatialPolygons with holes")
Shared vertices
This time we have a separate polygon “island”, but it’s actually sharing an edge.
p4 <- cbind(x = c(0.69, 0.8, 1.1, 1.23, 0.69),
y = c(0, 0.6, 0.63, 0.3, 0))
pp <- rbind(p1, NA, p2[nrow(p2):1, ])
plot(rbind(pp, p4), cex = 1.3, main = "two polygons, shared edge")
polypath(pp, col = "grey")
polypath(p4, col = "firebrick")
It now makes sense to remove the duplicate X/Y coordinates, in a complicated map there can be many such “shared vertices”.
Our code might look like this:
## unique vertices
coords <- cbind(x = c(0, 0, 0.75, 1, 0.5, 0.8, 0.69, 0.2, 0.2, 0.3, 0.5, 0.5, 1.1, 1.23),
y = c(0, 1, 1, 0.8, 0.7, 0.6, 0, 0.2, 0.4, 0.6, 0.4, 0.2, 0.63, 0.3))
## indexes for each polygon part
i1 <- c(1, 2, 3, 4, 5, 6, 7, 1)
i2 <- rev(c(8, 9, 10, 11, 12, 8)) ## holes are reversed for winding rule
i3 <- c(7, 6, 13, 14, 7)
plot(rbind(p1, p2, p4), main = "polygons, indexed from unique coords")
polypath(coords[c(i1, NA, i2), ], col = "grey")
polypath(coords[c(i3), ], col = "firebrick")
Obviously we want easy tools to do this for us.
- The
sppackage does not work this way. - It’s not easy, but we can make it much easier.
These indexes are not helpful once we subset a map - if we only want the second polygon, we’d have to keep all the vertices for the index to make sense.
How far can we take this?
library(rgl)
library(rglwidget)
rglpoly1 <- extrude3d(cbind(p1, 0), thickness = 0.2)
rglpoly2 <- extrude3d(cbind(p3, 0), thickness = 0.15)
shade3d(rglpoly1, col = "grey", alpha = 0.2)
shade3d(rglpoly2, col = "firebrick", alpha = 0.2)
subid <- currentSubscene3d()
rglwidget(elementId="simpleExtrude", width = 700, height = 700)
Polygon storage
For GIS it is like 2D plotting in R - “lists of matrix coordinates”
For 3D vis and modelling - “vertices indexed by line segments”
Lots of details and gotchas in-between!
We can flip between these models for greater flexibility.
Thanks!
Links to code etc.
Ecosystem Modelling, Australian Antarctic Division
Join Data Science Hobart (DaSH) for a weekly get-together at IMAS: http://datasciencehobart.github.io