r - Voronoi diagram polygons enclosed in geographic borders

Question

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.

Answer 1

You should be able to use the spatstat function dirichlet 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):

library(maps)
library(maptools)
library(spatstat)
library(rgeos)

counties <- map('county', c('maryland,carroll', 'maryland,frederick', 
                            'maryland,montgomery', 'maryland,howard'), 
                fill=TRUE, plot=FALSE)
# fill=TRUE is necessary for converting this map object to SpatialPolygons
countries <- gUnaryUnion(map2SpatialPolygons(counties, IDs=counties$names,
                                 proj4string=CRS("+proj=longlat +datum=WGS84")))
W <- as(countries, "owin")

Then you just put the five points in the ppp format, make the Dirichlet tesslation and calctulate the areas:

X <- ppp(x=c(-77.208703, -77.456582, -77.090600,  -77.035668, -77.197144),
         y=c(39.188603, 39.347019, 39.672818, 39.501898, 39.389203), window = W)

y <- dirichlet(X) # Dirichlet tesselation
plot(y) # Plot tesselation
plot(X, add = TRUE) # Add points
tile.areas(y) #Areas