I am trying to create Voronoi polygons (aka Dirichlet tessellations or Thiessen polygons) within a fixed geographic region for a set of points. However, I am having trouble finding a method in R that will bound the polygons within the map borders. My main goal is to get accurate area calculations (not simply to produce a visual plot). For example, the following visually communicates what I'm trying to achieve:
library(maps)
library(deldir)
data(countyMapEnv)
counties <- map('county', c('maryland,carroll','maryland,frederick', 'maryland,montgomery', 'maryland,howard'), interior=FALSE)
x <- c(-77.208703, -77.456582, -77.090600, -77.035668, -77.197144)
y <- c(39.188603, 39.347019, 39.672818, 39.501898, 39.389203)
points(x,y)
vt <- deldir(x, y, rw=counties$range)
plot(vt, wlines="tess", lty="solid", add=TRUE)
which produces the following:
Conceptually I want to intersect counties
with vt
which should provide a set of polygons bounded by the county borders and accurate area calculations for each. Right now, vt$summary
provides area calculations for each polygon, but they are obviously overstated for all but the one interior polygon, and deldir()
appears to only accept rectangular enclosings for its rw
argument. I am new to R's geospacial capabilities, so am open to other approaches beyond what I outlined above.
You should be able to use the
spatstat
functiondirichlet
for this.The first task is to get the counties converted into a spatstat observation window of class
owin
(code partially based on the answer by @jbaums):Then you just put the five points in the
ppp
format, make the Dirichlet tesslation and calctulate the areas:Once we have both the Voronoi polygons and the
counties
asSpatialPolygons
objects, we can achieve this with the help ofgIntersection
.First, let's load some necessary libraries and prepare your data.
Now we can convert our
counties
map
object toSpatialPolygons
withmap2SpatialPolygons
from themaptools
package. I've wrapped it inrgeos::gUnaryUnion
to combine the four polygons into a single polygon (otherwise we'd have internal boundaries plotted down the track). I've also added the relevant projection.For converting the
deldir
object to aSpatialPolygons
object, there's a nice function that I referred to here (hat-tip to Carson Farmer) and which @Spacedman subsequently modified (to clip to a given extent) and posted here.Now we can go ahead and use this to create our Voronoi
SpatialPolygons
.All that's left to do now is intersect the two geometries - the bread and butter of
rgeos
: