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.

Mord

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

Snowber

(+1) Ha, very nice indeed - a lot less stuffing around. Feel free to add the info about converting the map to a SpatialPolygon from my answer. – Tropology

I guess so, but at the same time the answer is missing the correct polygons, but I hope it sufficient for @Bryan. – Snowber

OK. I will try to include the info about converting the map to a SpatialPolygon from @Tropology answer. – Snowber

@EgeRubak and @jbaums, this is great, thanks so much. I had tried the spatstat approach but I was missing the conversion of the map to the owin class, so these are both super helpful. However, any idea why the first tile of the tessellation has two sets of bdry (i.e., y$tiles[[1]]$bdry lists one tiny polygon in upper left corner in addition to the main one)? I noticed this when I was attempting to loop through each tile to do area overlap calculations with a different set of polygons. – Mord

Also, I edited the above to include @jbaums' gUnaryUnion usage, as I was getting weird county lines when plotting another set of data. – Mord

I moved the gUnaryUnion command to the correct line (hopefully). I can't reproduce the strange y$tiles[[1]]$bdry behaviour you are talking about. Does it happen with the code above, and what is your spatstat version? – Snowber

I'm using spatstat 1.37-0. However, the y$tiles[[1]]$bdry issue seems to be solved by using gUnaryUnion. If you run it without gUnaryUnion you'll see the extra polygon in that first tile. – Mord

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)

Tropology

How do I fill colors to each tile of the merged map? The deldir package takes method plot.tile.list but it fills colors over the whole square. – Fluorene

Awesome, exactly what I needed! – Xylography