7

たとえば、米国の郡や PUMA など、同じ地域をカバーするシェイプ ファイルが 2 セットあるとします。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)

ここに画像の説明を入力 これを達成するために rgeos パッケージの g* 関数の多くを素朴に試してみましたが、前進していません。このタスクの優れた機能または素晴らしいトリックを知っている人はいますか? ありがとうございました！

[より良いタイトルの提案も受け付けています]

4

3 に答える 3

3

封じ込めのためのスマートな一連のテストでこれを行うことができると思います。これにより、2つのパーツが取得されます。containsと、、containsと。のpuma1単純county1なペアのケースです。county2puma2county8county3

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)

}

そこのプロットはばかげて冗長ですが、それは始まりを示していると思います。どのテストに最も小さなパーツが蓄積されているかを見つけて、それらを結合する必要があります。申し訳ありませんが、私は終わらせるために努力をしていませんが、これはうまくいくと思います。

于 2012-09-06T09:41:50.463 に答える
3

これは比較的簡潔なソリューションです。

  • rgeos::gRelate()重複しているが、各郡を完全に包含/カバーしていないピューマを識別するために使用しexample(gRelate)ます。(ティム・リフへのht)

  • RBGL::connectedComp()マージする必要があるピューマのグループを識別するために使用します。( RBGLパッケージのインストールに関する指針については、この SO の質問に対する私の回答を参照してください。)

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

ここに画像の説明を入力

于 2012-09-08T02:36:31.217 に答える
1

多くの試行錯誤の後、@mdsummer によるヒントに沿って、次の洗練されていないソリューションを思いつきましたが、チェックを追加し、冗長なマージされたポリゴンを削除し、チェックしました。うわぁ。誰かが私がやったことを取り、それをよりきれいにすることができれば、私はこれではなくその答えを受け入れます。これは可能な限り避けたいです:

library(rgeos)
pbyc <- gCovers(Pumas, counties, byid = TRUE) | 
        gContains(Pumas, counties, byid = TRUE) | 
        gOverlaps(Pumas, counties, byid = TRUE) |

        t(gCovers(counties, Pumas, byid = TRUE) | 
            gContains(counties, Pumas, byid = TRUE) |  
            gOverlaps(counties, Pumas, byid = TRUE))


## row / col pairs of where test is TRUE for Pumas or counties
pbyc.pairs <-  cbind(row(pbyc)[pbyc], col(pbyc)[pbyc])

Potentials <- apply(pbyc.pairs, 1, function(x,counties,Pumas){
     gUnion(counties[x[1], ], Pumas[x[2], ])
    }, counties = counties, Pumas= Pumas)
for (i in 1:length(Potentials)){
  Potentials[[i]]@polygons[[1]]@ID <- paste0("p",i)
}
Potentials <- do.call("rbind",Potentials)

# remove redundant polygons:
Redundants <- gEquals(Potentials, byid = TRUE)
Redundants <- row(Redundants)[Redundants & lower.tri(Redundants)]
Potentials <- Potentials[-c(Redundants),]

# for each Potential summary polygon, see which counties and Pumas are contained:
keep.i <- vector(length=length(Potentials))

for (i in 1:length(Potentials)){
  ctyblocki <- gUnionCascaded(counties[c(gCovers(Potentials[i, ], counties, byid = TRUE)), ])
  Pumablocki <- gUnionCascaded(Pumas[c(gCovers(Potentials[i, ], Pumas, byid = TRUE)), ]) 
  keep.i[i] <- gEquals(ctyblocki, Potentials[i, ]) & gEquals(Pumablocki, Potentials[i, ])    
}
# what do we have left?
NewUnits <- Potentials[keep.i, ]

plot(NewUnits)

ここに画像の説明を入力

于 2012-09-06T17:00:37.480 に答える