R：计算重叠多边形面积答案

【问题标题】：R: Calculating overlap polygon areaR：计算重叠多边形面积
【发布时间】：2014-11-22 16:05:33
【问题描述】：

我在尝试计算多边形之间的重叠区域时遇到了一些困难。例如，对于瑞士，我有关于该国各州行政边界和ethnic groups 的数据。

## Libraries
library(sp)
library(spdep)
library(maptools)
library(rgeos)

## Data
print(load(url("http://gadm.org/data/rda/CHE_adm1.RData")))
greg<-readShapeSpatial("raw_data/GREG.shp",
                     proj4string=CRS(" +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0")) 
switzerland<-greg[greg$COW==225,]

## Identical projections
proj<-CRS(" +proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0")
switzerland<-spTransform(switzerland,proj)

绘制数据，我们得到如下所示的内容：

## Plot data
par(mar=c(1,1,1,1))
plot(gadm,col="grey80")
plot(switzerland,add=TRUE,lwd=2,border="red")

我们可以看到，民族的边界并不完全遵循国界，但已经足够好了。我要做的是为每个州计算民族的数量，同时考虑到州内民族的面积。所以对于东部的格劳宾登州，我想知道德国瑞士人、意大利瑞士人、Rhaetoromaniens 等占据的区域。

在stackoverflow上阅读了一些类似的问题，我认为gIntersection是要使用的命令，但这给了我以下错误：

int<-gIntersection(gadm,switzerland,byid=TRUE) # Updated


Error in RGEOSBinTopoFunc(spgeom1, spgeom2, byid, id, drop_lower_td, "rgeos_intersection") : 
TopologyException: no outgoing dirEdge found at 7.3306802801147688 47.439399101791921
In addition: Warning message:
In RGEOSBinTopoFunc(spgeom1, spgeom2, byid, id, drop_lower_td, "rgeos_intersection") :
spgeom1 and spgeom2 have different proj4 strings

不完全确定导致第二次警告任务的原因

identicalCRS(gadm,switzerland)
[1] TRUE

有人对我如何计算各州和种族之间的重叠有什么建议吗？

更新：可能的解决方案

这可能是一个可能的解决方案，尽管不同 proj4 字符串上的警告仍然存在。另请注意，由于民族形状文件未正确遵循国界，因此存在一些测量误差（例如在阿尔瓜）。

## Possible solution
int<-gIntersection(gadm,switzerland,byid=T) # Multiple polygons per province
n<-names(int)
n<-data.frame(t(data.frame(strsplit(n," ",fixed=TRUE)))) 

colnames(n)[1:2]<-c("id0","ethnic.group")
n$area<-sapply(int@polygons, function(x) x@area)
a<-data.frame(id0=row.names(gadm),total.area=gadm$Shape_Area)
df<-merge(n,a,all.x=TRUE)

df$share.area<-df$area/df$total.area*100

【问题讨论】：

我认为你需要使用rgdal 包。你查看spTransform的帮助页面了吗？它说sp 只是定义了泛型，但要实际更改您需要rgdal 的坐标。
我确实检查了spTransform 帮助页面并认为我使用的语法是正确的。加载rgdal 后，我确实得到了同样的错误。我可以在这里忽略一些东西吗？
我不能说是什么导致了这个错误，因为我无权访问 GREG.* 文件。它们在某处可用吗？
是的，如果您点击 ethnic groups 文本中的蓝色链接，您将进入需要下载“GREG shapefile”的数据网页。
当您尝试使用整个gIntersects 整个gadm 和整个switzerland 对象时，似乎会出现错误。我尝试了res<-matrix(FALSE,nrow=nrow(gadm@data),ncol=nrow(switzerland@data)); for (i in 1:nrow(gadm@data)) for (j in 1:nrow(switzerland@data)) res[i,j]<-gIntersects(gadm[i,],switzerland[j,])，一切正常，除了警告部分（我怀疑这不应该是问题）。上面的行表示第 i 个州是否存在第 j 个民族。这就是你要找的东西吗？

标签： r polygon spatial

【解决方案1】：

这是一种与您的方法略有不同（但仅略有不同）的方法。

对switzerland@data 的检查显示，虽然有 11 个FeatureIDs（代表种族），但只有 4 个唯一命名的种族（德国、意大利、法国、瑞士和雷托罗曼尼亚人）。所以下面的结果是基于名称，而不是 ID。

library(rgeos)    # for gIntersection(...), etc.
library(rgdal)    # for readOGR(...)

setwd("<directory to accept your files>")
CH.1903 <- "+proj=somerc +lat_0=46.95240555555556 +lon_0=7.439583333333333 +k_0=1 +x_0=600000 +y_0=200000 +ellps=bessel +towgs84=674.374,15.056,405.346,0,0,0,0 +units=m +no_defs"

print(load(url("http://gadm.org/data/rda/CHE_adm1.RData")))
gadm <- spTransform(gadm, CRS(CH.1903))
download.file("http://www.icr.ethz.ch/data/other/greg/GREG.zip","GREG.zip")
unzip("GREG.zip")
greg <- readOGR(dsn=".",layer="GREG")
greg <- spTransform(greg[greg$COW==225,],CRS(CH.1903))

gadm.ids <- gadm@data$ID_1               # Canton Nr.
greg.ids <- unique(greg@data$G1SHORTNAM) # ethnicity
get.area <- Vectorize(function(adm,reg){
  int <- gIntersection(gadm[gadm$ID_1==adm,],greg[greg$G1SHORTNAM==reg,],byid=TRUE)
  if (length(int)==0) return(0)
  gArea(int)
})
result <- outer(gadm.ids,greg.ids,FUN=get.area)
rownames(result) <- gadm.ids
colnames(result)  <- greg.ids
result <- as.data.frame(result)
totals <- rowSums(result)
result <- result/totals
result$totals <- totals/1e6
result$land.area <- sapply(rownames(result),function(p)gArea(gadm[gadm$ID_1==p,])/1e6)
result
#     German Swiss French Swiss Italian Swiss Rhaetoromanians     totals  land.area
# 531  1.000000000   0.00000000    0.00000000    0.000000e+00 1363.27027 1403.22192
# 532  1.000000000   0.00000000    0.00000000    0.000000e+00  244.32279  244.32279
# 533  1.000000000   0.00000000    0.00000000    0.000000e+00  172.40341  172.40341
# 534  1.000000000   0.00000000    0.00000000    0.000000e+00  522.12943  525.73449
# 535  1.000000000   0.00000000    0.00000000    0.000000e+00   70.03116   84.06481
# 536  0.902128658   0.09787134    0.00000000    0.000000e+00 5927.90740 5927.90847
# 537  0.188415744   0.81158426    0.00000000    0.000000e+00 1654.28729 1654.28729

在这里，我们将两个 shapefile 转换为 CH.1903，这是一个以瑞士为中心的墨卡托投影，单位为米。然后我们确定Canton Nrs。和种族，并使用 outer(...) 在两个列表中循环，使用 gArea(...) 计算以平方公里 (sq.m/1e6) 为单位的交叉区域。最终结果每个州有一行，每个种族的百分比基于土地面积。 $totals 是每个州的相交面积之和，$land.area 是该州的总地理土地面积。因此，这可以让您了解由于种族 shapfile 和 gadm shapefile 之间不完全重叠而导致的错误。

【讨论】：