如何通过简化 R 中的两个 SpatialPolygonsDataFrame 对象来创建新的多边形？

假设我有两组形状文件，它们覆盖同一区域，并且经常但并不总是共享边界，例如美国各县和 PUMA。我想定义一个新的多边形规模，它使用 PUMA 和县作为原子构建块，即两者都不能分割，但我仍然希望有尽可能多的单位。这是一个玩具示例：

library(sp)
# make fake data
# 1) counties:
Cty <- SpatialPolygons(list(
    Polygons(list(Polygon(cbind(x=c(0,2,2,1,0,0),y=c(0,0,2,2,1,0)), hole=FALSE)),"county1"),
    Polygons(list(Polygon(cbind(x=c(2,4,4,3,3,2,2),y=c(0,0,2,2,1,1,0)),hole=FALSE)),"county2"),
    Polygons(list(Polygon(cbind(x=c(4,5,5,4,4),y=c(0,0,3,2,0)),hole=FALSE)),"county3"),
    Polygons(list(Polygon(cbind(x=c(0,1,2,2,0,0),y=c(1,2,2,3,3,1)),hole=FALSE)),"county4"),
    Polygons(list(Polygon(cbind(x=c(2,3,3,4,4,3,3,2,2),y=c(1,1,2,2,3,3,4,4,1)),hole=FALSE)),"county5"),
    Polygons(list(Polygon(cbind(x=c(0,2,2,1,0,0),y=c(3,3,4,5,5,3)),hole=FALSE)),"county6"),
    Polygons(list(Polygon(cbind(x=c(1,2,3,4,1),y=c(5,4,4,5,5)),hole=FALSE)),"county7"),
    Polygons(list(Polygon(cbind(x=c(3,4,4,5,5,4,3,3),y=c(3,3,2,3,5,5,4,3)),hole=FALSE)),"county8")
))

counties <- SpatialPolygonsDataFrame(Cty, data = data.frame(ID=paste0("county",1:8),
            row.names=paste0("county",1:8),
            stringsAsFactors=FALSE)
)
# 2) PUMAs:
Pum <- SpatialPolygons(list(
            Polygons(list(Polygon(cbind(x=c(0,4,4,3,3,2,2,1,0,0),y=c(0,0,2,2,1,1,2,2,1,0)), hole=FALSE)),"puma1"),
            Polygons(list(Polygon(cbind(x=c(4,5,5,4,3,3,4,4),y=c(0,0,5,5,4,3,3,0)),hole=FALSE)),"puma2"),
            Polygons(list(Polygon(cbind(x=c(0,1,2,2,3,3,2,0,0),y=c(1,2,2,1,1,2,3,3,1)),hole=FALSE)),"puma3"),
            Polygons(list(Polygon(cbind(x=c(2,3,4,4,3,3,2,2),y=c(3,2,2,3,3,4,4,3)),hole=FALSE)),"puma4"),
            Polygons(list(Polygon(cbind(x=c(0,1,1,3,4,0,0),y=c(3,3,4,4,5,5,3)),hole=FALSE)),"puma5"),
            Polygons(list(Polygon(cbind(x=c(1,2,2,1,1),y=c(3,3,4,4,3)),hole=FALSE)),"puma6")
    ))
Pumas <- SpatialPolygonsDataFrame(Pum, data = data.frame(ID=paste0("puma",1:6),
            row.names=paste0("puma",1:6),
            stringsAsFactors=FALSE)
)

# desired result:
Cclust <- SpatialPolygons(list(
            Polygons(list(Polygon(cbind(x=c(0,4,4,3,3,2,2,1,0,0),y=c(0,0,2,2,1,1,2,2,1,0)), hole=FALSE)),"ctyclust1"),
            Polygons(list(Polygon(cbind(x=c(4,5,5,4,3,3,4,4),y=c(0,0,5,5,4,3,3,0)),hole=FALSE)),"ctyclust2"),
            Polygons(list(Polygon(cbind(x=c(0,1,2,2,3,3,4,4,3,3,2,2,0,0),y=c(1,2,2,1,1,2,2,3,3,4,4,3,3,1)),hole=FALSE)),"ctyclust3"),
            Polygons(list(Polygon(cbind(x=c(0,2,2,3,4,0,0),y=c(3,3,4,4,5,5,3)),hole=FALSE)),"ctyclust4")
    ))
CtyClusters <- SpatialPolygonsDataFrame(Cclust, data = data.frame(ID = paste0("ctyclust", 1:4),
            row.names = paste0("ctyclust", 1:4),
            stringsAsFactors=FALSE)
)

# take a look
par(mfrow = c(1, 2))
plot(counties, border = gray(.3), lwd = 4)
plot(Pumas, add = TRUE, border = "#EEBB00", lty = 2, lwd = 2)
legend(-.5, -.3, lty = c(1, 2), lwd = c(4, 2), col = c(gray(.3), "#EEBB00"),
    legend = c("county line", "puma line"), xpd = TRUE, bty = "n")
text(coordinates(counties), counties@data$ID,col = gray(.3))
text(coordinates(Pumas), Pumas@data$ID, col = "#EEBB00",cex=1.5)
title("building blocks")
#desired result:
plot(CtyClusters)
title("desired result")
text(-.5, -.5, "maximum units possible,\nwithout breaking either PUMAs or counties",
    xpd = TRUE, pos = 4)

I've naively tried many of the g* functions in the rgeos package to achieve this and have not made headway. Does anyone know of a nice function or awesome trick for this task? Thank you!

[我也愿意接受关于更好标题的建议]

这是一个相对简洁的解决方案：

• Uses rgeos::gRelate()识别重叠但不完全包含/覆盖每个县的美洲狮。要了解它的作用，请运行example(gRelate)并看到这个维基百科页面。 （致蒂姆·里夫）

• Uses RBGL::connectedComp()以确定应该合并的美洲狮群体。 （有关安装的指导RBGL包，请参阅我的回答这个问题.)

• Uses rgeos::gUnionCascaded()合并指定的美洲狮。

  library(rgeos)
  library(RBGL)
  
  ## Identify groups of Pumas that should be merged
  x <- gRelate(counties, Pumas, byid=TRUE)
  x <- matrix(grepl("2.2......", x), ncol=ncol(x), dimnames=dimnames(x))
  k <- x %*% t(x)
  l <- connectedComp(as(k, "graphNEL"))
  
  ## Extend gUnionCascaded so that each SpatialPolygon gets its own ID.
  gMerge <- function(ii) {
      x <- gUnionCascaded(Pumas[ii,])
      spChFIDs(x, paste(ii, collapse="_"))
  }
  
  ## Merge Pumas as needed
  res <- do.call(rbind, sapply(l, gMerge))
  
  plot(res)