Voronoi图多边形包含在地理边界内

Question

我试图在一组固定的地理区域内创建Voronoi多边形(又名Dirichlet镶嵌或Thiessen多边形).但是,我在R中找到一个方法会遇到地图边界内的多边形.我的主要目标是获得准确的面积计算(不仅仅是生成视觉图).例如,以下内容直观地传达了我想要实现的目标:

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)

产生以下内容:

从概念上讲,我想counties与之相交vt,应该提供一组由县界限定的多边形,并为每个多边形进行准确的面积计算.现在,vt$summary为每个多边形提供面积计算,但除了一个内部多边形之外,它们显然被夸大了,并且deldir()似乎只接受其rw参数的矩形包围.我是R的geospacial能力的新手,所以我可以接受超出我上面概述的其他方法.

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格式化五个点,制作Dirichlet tesslation并制作区域:

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多边形和countiesas SpatialPolygons对象,我们就可以借助于实现这一目标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对象SpatialPolygons与map2SpatialPolygons从maptools包.我将它包裹rgeos::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)