说我有两组形状文件覆盖同一区域，并且经常但不总是共享边界，例如美国县和 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) 

请您参考如下方法:

我认为您可以通过一组智能的遏制测试来做到这一点。这得到了你的两个部分，简单的配对案例，其中 puma1包含 county1和 county2 , 和 puma2包含 county8和 county3 .

library(rgeos) 
 
## pumas by counties 
pbyc <- gContains(Pumas, counties, byid = TRUE) 
 
## row / col pairs of where contains test is TRUE for Pumas 
pbyc.pairs <-  cbind(row(pbyc)[pbyc], col(pbyc)[pbyc]) 
 
par(mfrow = c(nrow(pbyc.pairs), 1)) 
 
for (i in 1:nrow(pbyc.pairs)) { 
plot(counties, col = "white") 
 
plot(gUnion(counties[pbyc.pairs[i,1], ], Pumas[pbyc.pairs[i,2], ]), col = "red", add = TRUE) 
 
}