r - 2 つの SpatialPolygonsDataFrame オブジェクトのベクトルオーバーレイを実行する方法は?

Question

私は 2 つの GIS レイヤーを持っており、それらSoilsを呼び出して、 s ( s)Parcelsとして保存しており、ここで説明する意味でそれらを「オーバーレイ」したいと考えています。SpatialPolygonsDataFrameSPDF

オーバーレイ操作の結果は、次のような新しい SPDF になります。

コンポーネントには、2 つのレイヤーの交差によって形成されたSpatialPolygonsポリゴンが含まれます。(オーバーヘッドプロジェクターに 2 つのマイラーを重ねることによって形成されるすべてのアトミックポリゴンを考えてみてください)。
このdata.frameコンポーネントは、各アトミックポリゴンが含まれるSoilsおよびポリゴンの属性を記録します。Parcels

私の質問:これを行う既存の R 関数はありますか? SpatialPolygons(コンポーネントを正しく取得し、2 つのレイヤーの交差から形成されるアトミックポリゴンを計算する関数を知りたいとさえ思っています。) rgeosには、少なくとも (1) を実行する関数が必要だと思いますが、ないようです...

これは、私が求めているものをより明確にするのに役立つ図であり、その後に、図に示されているレイヤーSoilsとParcelsレイヤーを作成するコードが続きます。

ここに画像の説明を入力

library(rgeos)

## Just a utility function to construct the example layers.
flattenSpatialPolygons <- function(SP) {
    nm <- deparse(substitute(SP))
    AA <- unlist(lapply(SP@polygons, function(X) X@Polygons))
    SpatialPolygons(lapply(seq_along(AA),
                           function(X) Polygons(AA[X], ID=paste0(nm, X))))
}

## Example Soils layer
Soils <-
local({
    A <- readWKT("MULTIPOLYGON(((3 .5,7 1,7 2,3 1.5,3 0.5), (3 1.5,7 2,7 3,3 2.5,3 1.5)))")
    AA <- flattenSpatialPolygons(A)
    SpatialPolygonsDataFrame(AA,
           data.frame(soilType = paste0("Soil_", LETTERS[seq_along(AA)]),
                      row.names = getSpPPolygonsIDSlots(AA)))
})

## Example Parcels layer
Parcels <-
local({
    B <- readWKT("MULTIPOLYGON(((0 0,2 0,2 3,0 3,0 0),(2 0,4 0,4 3,2 3,2 0)),((4 0,6 0,6 3,4 3,4 0)))")
    BB <- flattenSpatialPolygons(B)
    SpatialPolygonsDataFrame(BB,
           data.frame(soilType = paste0("Parcel_", seq_along(BB)),
                      row.names = getSpPPolygonsIDSlots(BB)))
})

score 7 · Accepted Answer

2014 年 1 月以降、ラスターパッケージには、union()これを簡単にする関数が含まれています。

library(raster)
Intersects <- Soils + Parcels  ## Shorthand for union(Soils, Parcels)

## Check that it work
data.frame(Intersects)
## soilType.1 soilType.2
## 1     Soil_A       <NA>
## 2     Soil_B       <NA>
## 3       <NA>   Parcel_1
## 4       <NA>   Parcel_2
## 5       <NA>   Parcel_3
## 6     Soil_A   Parcel_2
## 7     Soil_A   Parcel_3
## 8     Soil_B   Parcel_2
## 9     Soil_B   Parcel_3
plot(Intersects, col = blues9)

ここに画像の説明を入力

score 4 · Accepted Answer

raster::union()2014 年 1 月以降、ここに投稿されたソリューションは、上記の公式に受け入れられた回答で示されている関数に完全に取って代わられました。

これにより、2 つの異なるレイヤーを重ね合わせて形成された「アトミック」ポリゴンが取得されSpatialPolygons、質問のパート 1 が解決されます。

library(rgeos)

gFragment <- function(X, Y) {
    aa <- gIntersection(X, Y, byid = TRUE)
    bb <- gDifference(X, gUnionCascaded(Y), byid = T)
    cc <- gDifference(Y, gUnionCascaded(X), byid = T)
    ## Note: testing for NULL is needed in case any of aa, bb, or cc is empty,
    ## as when X & Y don't intersect, or when one is fully contained in the other
    SpatialPolygons(c(if(is.null(aa)) NULL else aa@polygons,
                      if(is.null(bb)) NULL else bb@polygons,
                      if(is.null(cc)) NULL else cc@polygons)) 
}

## Try it out
Fragments <- gFragment(Parcels, Soils)
plot(Fragments, col=blues9)

そして、これは、上記で出力された各「アトミック」ポリゴンによってオーバーレイされた2つの入力レイヤー内のポリゴンのID（存在する場合）を抽出し、gFragment()私の質問のパート2を解決します。

getAttributes <- function(Fragments, Layer1, Layer2, eps = 0) {
    ## Function to extract attributes from polygon layers
    ## overlain by fragments
    OVER <- function(AA, BB) {
        X <- gRelate(AA, BB, byid = TRUE, pattern="2********")
        ii <- sapply(seq_len(ncol(X)),
                    function(i) {
                        A <- which(X[,i])
                        if(!length(A)) NA else A
                    })
        rownames(X)[ii]
    }
    ## First need to (very slightly) trim Fragments b/c otherwise they
    ## tend to (very slightly) overlap adjacent ring(s)
    Frags <- gBuffer(Fragments, width = -eps, byid = TRUE)
    ## Construct data.frame of attributes
    df <- data.frame(OVER(Frags, Layer1), 
                     OVER(Frags, Layer2),
                     row.names = names(Fragments))
    names(df) <- c(deparse(substitute(Layer1)), deparse(substitute(Layer2)))
    ## Add data.frame to SpatialPolygons object
    SpatialPolygonsDataFrame(Fragments, data=df)
}

FragmentsDF <- getAttributes(Fragments = Fragments,
                             Layer1 = Parcels,
                             Layer2 = Soils)

## A few ways to examine the results
data.frame(FragmentsDF, row.names=NULL)
#   Parcels Soils
# 1      B2    A1
# 2      B2    A2
# 3      B3    A1
# 4      B3    A2
# 5      B1  <NA>
# 6      B2  <NA>
# 7      B3  <NA>
# 8    <NA>    A1
# 9    <NA>    A2
spplot(FragmentsDF, zcol="Soils", col.regions=blues9[3:4])
spplot(FragmentsDF, zcol="Parcels", col.regions=grey.colors(3))

編集：

入力ポリゴンのいずれかに id の名前が付いている場合、このコードは失敗する可能性があることに注意してください"1"。その場合、1 つの回避策は、おそらくのような方法で ID の名前を変更することですParcels <- spChFIDs(Parcels, paste0("pp", row.names(Parcels)))。

score 2 · Accepted Answer

これが私のクラックです。これは、パーセル (パーセル -> 土壌) データを含むピースのリストを提供するだけです。Soils オブジェクトからアトリビュートを取得し、「差異」を使用して同様の作業を行い、その逆 (Soils->Parcels) を行って、あらゆる種類のオーバーラップリレーションを作成する必要があります。

intersects <- list()

## find all intersections (NULLs do nothing to the result)
for (i in 1:nrow(Soils)) {
  for (j in 1:nrow(Parcels)) {
    intersects[[sprintf("%sx%s", i, j)]] <- gIntersection(Soils[i,],
                                                          Parcels[j,]) 
  }
}

result <- list()
## let's try Parcels, transfer data attributes to new pieces
for (i in 1:nrow(Parcels)) {
  for (j in seq_along(intersects))
   if(gContains(Parcels[i,], intersects[[j]])) {
     result <- c(result, SpatialPolygonsDataFrame(intersects[[j]],     as.data.frame(Parcels[i,]), match.ID = FALSE))

   }
}


## plot
plot(Parcels, xlim = range(c(bbox(Parcels)[1,], bbox(Soils[1,]))),
     ylim = range(c(bbox(Parcels)[2,], bbox(Soils[2,]))))
plot(Soils, add = TRUE)

cols <- colorRampPalette(c("lightblue", "darkblue"))(length(result))
for (i in 1:length(result)) plot(result[[i]], col = cols[i], add = TRUE)
for (i in 1:length(result)) text(coordinates(result[[i]]), label =     as.data.frame(result[[i]])[,"soilType"])

score 0 · Accepted Answer

基本的な考え方は次のとおりです (これは、区画に対してネストされたループとして実行し、次に土壌に対して実行します。関数が記述されている方法でベクトル化できるかどうかはわかりませんg*):

i <- 2
j <- 2
pieces <- list()
pieces[["int"]] <- gIntersection(Parcels[i,],Soils[j,])
pieces[["diff1"]] <- gDifference(Parcels[i,],Soils[j,])
pieces[["diff2"]] <- gDifference(Soils[j,],Parcels[i,])
plot(Parcels)
plot(Soils,add=TRUE)
plot(pieces[["int"]],col="red",add=TRUE)
plot(pieces[["diff1"]],col="blue",add=TRUE)
plot(pieces[["diff2"]],col="green",add=TRUE)

プロット

それはあなたが始めるのに十分なはずです。残りは、断片を追跡し、それらすべてを 1 つの大きな SPDF にマージしながらループするだけです。

よりベクトル化された代替アプローチは次のとおりです。

pieces <- list()
pieces[["int"]] <- gIntersection(Parcels,Soils,byid=TRUE)
pieces[["diff1"]] <- gDifference(Parcels,Soils,byid=TRUE)
pieces[["diff2"]] <- gDifference(Soils,Parcels,byid=TRUE)

これにより、何らかの理由で、実際に存在するよりも多くのピースが得られます。次に、戻ってそれらをバインドし、余分な部分を選別する必要がありgEqualsます。

r - 2 つの SpatialPolygonsDataFrame オブジェクトのベクトル オーバーレイを実行する方法は?

4 に答える 4

Related

Reference

r - 2 つの SpatialPolygonsDataFrame オブジェクトのベクトルオーバーレイを実行する方法は?