在 R 中找到最佳匹配的重叠多边形

Question

我有两个 shapefile，我使用 readOGR() 作为 SpatialPolygonsDataFrame 对象读入 R。两者都是具有不同内部边界的新西兰地图。一个有大约 70 个多边形，代表领土当局的边界；另一个有大约 1900 个代表区域单位。

我的目标 - 一个更大项目的一个令人讨厌的基本部分 - 是使用这些地图生成一个参考表，可以查找一个区域单位并返回它主要位于哪个领土权威。我可以使用 over() 来查找哪些多边形重叠，但在许多情况下，区域单位似乎，至少在一小部分，在多个领土当局内 - 尽管查看个别案例表明，通常 90% 以上的区域单位在一个单一的领土当局内。

是否有一个现成的手段来完成 over() 的作用，但它不仅可以识别（或什至不能）所有重叠的多边形，而且在每种情况下，几个重叠的多边形中哪个重叠最多？

Answer 1

我认为地理信息系统 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 包。）

Answer 2

这是完成这项工作的代码，借鉴了@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(au@data, 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 <- ta@data[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, col=ta@data$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()