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

您应该能够为此使用该spatstat功能dirichlet。

第一个任务是将县转换为类的 spatstat 观察窗口owin（代码部分基于@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")

然后您只需将五个点放入ppp格式中，进行狄利克雷曲面细分并计算面积：

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

一旦我们拥有了 Voronoi 多边形和countiesasSpatialPolygons对象，我们就可以在gIntersection.

首先，让我们加载一些必要的库并准备您的数据。

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))

现在我们可以将我们的counties map对象转换为包中的SpatialPolygonswith 。我已经将它包裹起来以将四个多边形组合成一个多边形（否则我们将在轨道上绘制内部边界）。我还添加了相关的投影。map2SpatialPolygonsmaptoolsrgeos::gUnaryUnion

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

为了将deldir对象转换为对象，我在此处SpatialPolygons提到了一个不错的函数（对Carson Farmer的帽子提示），@Spacedman 随后对其进行了修改（以剪辑到给定范围）并在此处发布。

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'))))  
}

现在我们可以继续使用它来创建我们的 Voronoi SpatialPolygons。

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

现在剩下要做的就是将两个几何图形相交 - 的面包和黄油rgeos：

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

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