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

Answer 2

Once we have both the Voronoi polygons and the counties as SpatialPolygons objects, we can achieve this with the help of gIntersection.

First, let's load some necessary libraries and prepare your data.

library(maptools)
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

p <- data.frame(x=c(-77.208703, -77.456582, -77.090600,  -77.035668, -77.197144),
                y=c(39.188603, 39.347019, 39.672818, 39.501898, 39.389203))

Now we can convert our counties map object to SpatialPolygons with map2SpatialPolygons from the maptools package. I've wrapped it in rgeos::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.

counties.sp <- gUnaryUnion(
  map2SpatialPolygons(counties, IDs=counties$names,
                      proj4string=CRS("+proj=longlat +datum=WGS84")))

For converting the deldir object to a SpatialPolygons 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.

voronoipolygons <- function(x, poly) {
  require(deldir)
  if (.hasSlot(x, 'coords')) {
    crds <- x@coords  
  } else crds <- x
  bb = bbox(poly)
  rw = as.numeric(t(bbox(poly)))
  z <- deldir(crds[,1], crds[,2],rw=rw)
  w <- tile.list(z)
  polys <- vector(mode='list', length=length(w))
  require(sp)
  for (i in seq(along=polys)) {
    pcrds <- cbind(w[[i]]$x, w[[i]]$y)
    pcrds <- rbind(pcrds, pcrds[1,])
    polys[[i]] <- Polygons(list(Polygon(pcrds)), ID=as.character(i))
  }
  SP <- SpatialPolygons(polys)

  SpatialPolygonsDataFrame(
    SP, data.frame(x=crds[,1], y=crds[,2], 
                   row.names=sapply(slot(SP, 'polygons'), 
                                    function(x) slot(x, 'ID'))))  
}

Now we can go ahead and use this to create our Voronoi SpatialPolygons.

v <- voronoipolygons(p, counties.sp)
proj4string(v) <- proj4string(counties.sp)

All that's left to do now is intersect the two geometries - the bread and butter of rgeos:

final <- gIntersection(counties.sp, v, byid=TRUE)

plot(final)
points(p, pch=20)