在R中找到最匹配的重疊多邊形

我有兩個使用readOGR（）作爲SpatialPolygonsDataFrame對象讀入R的shapefile。兩者都是具有不同內部邊界的新西蘭地圖。一個大約有70個多邊形代表領土權力邊界;另一個大約有1900個代表區域單位。在R中找到最匹配的重疊多邊形

我的目標 - 一個更大項目的令人討厭的基本部分 - 是使用這些地圖來生成一個參考表，它可以查找一個區域單位並返回它最主要的地域權限。我可以使用over（）來發現哪些多邊形重疊，但在許多情況下，區域單位似乎至少在多個領土當局內 - 儘管查看個別案例表明通常90％以上的區域單位屬於單一地區當局。

是否有準備好的手段意味着做什麼over（），但它不僅可以識別所有重疊的多邊形（甚至不是），而且可以識別幾個重疊多邊形中哪一個在每種情況下重疊最多？

來源

2013-01-08 Peter Ellis

我想你想要的東西已經覆蓋的地理信息系統SE：

https://gis.stackexchange.com/questions/40517/using-r-to-calculate-the-area-of-multiple-polygons-on-a-map-that-intersect-with?rq=1

特別是如果你的領土多邊形T1，T2，T3等和你的多邊形您要分類是A，可能想在A和T1的gIntersection上使用gArea，然後A和T2，然後A和T3等，然後選擇哪個具有最大面積。（您將需要rgeos包）。

來源

2013-01-08 04:32:23 Silverfish

感謝代碼 - gArea和gIntersection在一起而形成的缺失環節我。這看起來應該起作用。如果是這樣，我會接受這個答案。 –

很高興聽到它。如果你能得到它的工作將是很好，如果你可以留下一個有效的代碼示例，因爲我不相信你是唯一一個試圖執行這樣一個顯而易見的重要任務，並努力尋找工具來做到這一點！ – Silverfish

謝謝@Silverfish - 我已經添加了一個可以完成這項工作的答案，並且應該適用於其他人。 –

這裏是做了工作，借鑑@蠹蟲的回答

library(sp) 
library(rgeos) 
library(rgdal) 

### 
# Read in Area Unit (AU) boundaries 
au <- readOGR("C:/Users/Peter Ellis/Documents/NZ", layer="AU12") 

# Read in Territorial Authority (TA) boundaries 
ta <- readOGR("C:/Users/Peter Ellis/Documents/NZ", layer="TA12") 

### 
# First cut - works ok when only one TA per area unit 
x1 <- over(au, ta) 
au_to_ta <- data.frame([email protected], TAid = x1) 

### 
# Second cut - find those with multiple intersections 
# and replace TAid with that with the greatest area. 

x2 <- over(au, ta, returnList=TRUE) 

# This next loop takes around 10 minutes to run: 
for (i in 1:nrow(au_to_ta)){ 
    tmp <- length(x2[[i]]) 
    if (tmp>1){ 
     areas <- numeric(tmp) 
     for (j in 1:tmp){ 
      areas[j] <- gArea(gIntersection(au[i,], ta[x2[[i]][j],])) 
      } 
#  Next line adds some tiny random jittering because 
#  there is a case (Kawerau) which is an exact tie 
#  in intersection area with two TAs - Rotorua and Whakatane 

     areas <- areas * rnorm(tmp,1,0.0001) 

     au_to_ta[i, "TAid"] <- x2[[i]][which(areas==max(areas))] 
    } 

} 


# Add names of TAs 
au_to_ta$TA <- [email protected][au_to_ta$TAid, "NAME"] 

#### 
# Draw map to check came out ok 
png("check NZ maps for TAs.png", 900, 600, res=80) 
par(mfrow=c(1,2), fg="grey") 
plot(ta, [email protected]$NAME) 

title(main="Original TA boundaries") 
par(fg=NA) 
plot(au, col=au_to_ta$TAid) 
title(main="TA boundaries from aggregated\nArea Unit boundaries") 
dev.off()

enter image description here

來源

2013-01-09 05:58:13

在R中找到最匹配的重疊多邊形

回答

相關問題