4

パッケージのgTouches関数は、rgeos「ジオメトリに共通の境界点が少なくとも 1 つあるが、内部点がない」かどうかをテストします。内部点に関連する基準なしで、 「ジオメトリに共通の境界点が少なくとも 1 つある」かどうかをテストする方法を探しています。

基本的なセットアップは次のとおりです。ほとんどが互いに埋め込まれている 2 つのシェープファイルがあります。大きな領域の境界にある小さな領域を持つファイル内のポリゴンを見つけたいです。ここに私がやろうとしていることを説明するグラフがあります:

plot(map2, col=NA, border='black', lwd=0.4)
plot(map1, col=NA, border='#666666', lwd=0.2, add=TRUE)

ここに画像の説明を入力

この図は、ニューヨーク州スタテン アイランドの国勢調査ブロックを示しています。大きな領域の 1 つにある緑色の強調表示は、識別したいブロックを示しています。より大きな領域の境界を共有または交差するもののみ (太い線)。大きなエリアの真ん中にあるブロックではありません。パッケージ内のwithgTouches(map2,map1, byid=TRUE)および他の機能を使用してこれを実行しようとしましたrgeosが、成功しませんでした。おそらく、「ジオメトリには共通の境界点が少なくとも 1 つありますが、内部の点はありません」という基準があるため、gTouches返されるだけです。FALSE基本的には、内部に関係なく「ジオメトリに共通の境界点が少なくとも 1 つある」かどうかをテストする関数を探しています。

フォローアップの質問は、相互境界線の長さを取得できるかどうかです。

データ:こことここから2 つのマップをダウンロードできます。どちらもrdsファイルなので、次のように開くことができます：

library('rgdal')
library('rgeos')
library('sp')
map1 = readRDS('map1.rds')
map2 = readRDS('map2.rds')
4

1 に答える 1

5

gIntersects()(学区の任意の部分と交差するすべての小さな多角形をgContainsProperly()検索する) と (学区の境界内に完全に含まれ、境界と交差しないすべての小さな多角形を検索する)の組み合わせを使用できます。次に、結果として得られる 2 つの論理行列を単純に組み合わせて、目的のポリゴンを識別します。

## Identify polygons that intersect but aren't fully contained within the
## school district whose polygon is given by SD = map2[13,]
SD <- map2[13,]
ii <- gIntersects(SD, map1, byid=TRUE) &
      !gContainsProperly(SD, map1, byid=TRUE)
ii <- apply(ii, 1, any)  ## Handy construct if both layers contain >1 polygon

## Plot that area, to show that results are  correct
plot(SD, col=NA, border='black')                  ## Establish plotted area
plot(map1, col=NA, border='#666666', lwd=0.2, add=TRUE)
plot(map1[ii,], col="lightgreen", add=TRUE)
plot(SD, col=NA, border='black', lwd=2, add=TRUE) ## Put SD boundary on top

ここに画像の説明を入力

編集 ：

しかし、それは正しくありません。上の地図に見られるように、学区の南西と南東の内部に沿って識別されるべきであった多くの小さなポリゴンが識別されていません。このような結果は、 rgeos操作でかなり頻繁に発生し、レイヤーのペア (または GEOS エンジンによるそれらの中間表現) のわずかな位置ずれから発生します。

gBuffer()解決策は、トポロジ クエリを実行する前に、レイヤの 1 つを少量バッファ アウトするために使用することです。ここでは、座標はメートル単位であり、試行錯誤の結果、問題を解決するには 20 メートルのバッファーでほぼ十分であることがわかりました。

## Expand every polygon in map1 by a 20-meter wide buffer 
map1_buff <- gBuffer(map1, byid=TRUE, width=20)

## and then use the buffered version of map1 in the topological queries
ii <- gIntersects(SD, map1_buff, byid=TRUE) &
      !gContainsProperly(SD, map1_buff, byid=TRUE)
ii <- apply(ii, 1, any)  ## Handy construct if both layers contain >1 polygon

## Plot that area, to show that results are  correct
plot(SD, col=NA, border='black')                  ## Establish plotted area
plot(map1, col=NA, border='#666666', lwd=0.2, add=TRUE)
plot(map1[ii,], col="lightgreen", add=TRUE)
plot(SD, col=NA, border='black', lwd=2, add=TRUE)

ここに画像の説明を入力

これでも海岸沿いのいくつかのポリゴンが欠けていますが、完全な解決策として、詳細レベルがよりよく一致するマップのペアを取得する必要があるかもしれません。バッファ サイズを大幅に大きくすると、この分析で誤検知が発生し始めます。たとえば、学区の北西の角にある真に内側のポリゴンがいくつか検出されます。

于 2013-11-14T18:42:55.970 に答える