如何在R中计算多部分多边形的面积

Question

我无法在 R 中计算具有多部分多边形（一个特征包含两个单独的多边形）的导入 shapefile 的面积。我注意到 ArcMap 给我的 shapefile 区域的值与raster::area 不同。为了弄清楚哪个程序给了我正确的面积，我将 shapefile 分解为单个部分并重新计算了两个单独多边形的面积：

library(raster)

> single_part <- shapefile("../Desktop/test/test_sp.shp")
> area(single_part)
[1] 575924.0 433409.8
> sum(area(single_part))
[1] 1009334
> 
> multi_part <-  shapefile("../Desktop/test/test_mp.shp")
> area(multi_part)
[1] 1018390

我现在意识到我知道这个问题，我应该总是将多边形要素类分解成单个部分，但是有人知道raster::area 如何计算多部分多边形的面积吗？我也尝试使用rgeos::gArea，但得到了相同的结果。有没有办法计算 R 中多部分多边形的面积？

我很想知道，因为它们很常见，我正在尝试从 ArcMap 中的所有分析切换到 R。

如果有帮助，这里是 shapefile 的图像： multipart poly shapefile

编辑于 2018 年 9 月 21 日添加 ------------------------- ---------------

这是link to the shapefile test_mp.shp

据我所知，问题似乎源于 R（相对于 ArcMap）如何解释这些漏洞。查看ArcMap display 和R display 之间的区别。出于某种原因，R 作为 shapefile 的一部分填充了这些孔，这一定是我对该区域进行不同计算的原因。 shapefile 是否有问题，或者我如何导入它？

Answer 1

显然，您名为“multi_part”的对象只有一个（多？）多边形，因为area 返回一个值。我在这里说明如何调查你所追求的：

library(raster)
d <- getData('GADM', country='Isle of Man', level=0)

area(d)
[1] 579672897

分成4个多边形（岛）

dd <- disaggregate(d)
a <- area(dd)
a
[1]     19424.12   2705442.41     25629.79 576922400.90
sum(a)
[1] 579672897

同一个区域，没有理由不一样。除非可能与多边形孔混淆。没有您的数据很难发表评论。

您可以将这些对象写入磁盘（见下文）并查看 ArcGIS 为您提供的区域（但请注意，此示例使用 lon/lat 坐标，我不确定 ArcGIS 是否可以计算这些区域）。

shapefile(d, "man.shp")

这是一个有洞和没有洞的案例：

p1 <- rbind(c(-180,-20), c(-140,55), c(10, 0), c(-140,-60), c(-180,-20))
p2 <- rbind(c(-150,-20), c(-100,-10), c(-110,20), c(-150,-20))

# two (overlapping) polygons (no hole)   
pol1 <- spPolygons(p1, p2, crs="+proj=utm +zone=1 +datum=WGS84")
# single polygon with hole
pol2 <- spPolygons(list(p1, p2), crs="+proj=utm +zone=1 +datum=WGS84")

a <- area(pol1) / 10e+9
b <- area(pol2) / 10e+9

a
#[1] 10925   800
sum(a)
#[1] 11725
a[1]-a[2]
#[1] 10125

b 匹配 a[1] - a[2]，正如预期的那样

b
#[1] 10125

我使用 ArcGIS 获得完全相同的结果，对属性表中的字段使用“计算几何”。