geosphere:: areaPolygon 的面积计算错误？

Question

我已经多次成功使用geosphere::areaPolygon()，但现在我遇到了一个问题。

我有一个多边形 sp2，其中包含另一个多边形 sp1。所以 sp2 的面积应该比 sp1 大。当我用areaPolygon() 计算每个区域时，我得到相反的结果。

areasp1 = 10133977
areasp2 = 9858811

我使用 gSymdifference 找到对称的不同多边形 sp3，它有

areasp3 = 275165.4

sp1 : 红色虚线
sp2 : 黑线
sp3 : 绿色虚线

我使用的代码：

library(sp)
library(geosphere)
library(rgeos)
library(raster)

s1 <- data.frame(lon= c(130.4327, 129.85127, 121.00775, 84.9176, 60.40123, 
58.97929, 58.55841, 94.95237),lat= c(30.25074, 29.83075, 23.7992, 28.25964, 
36.89905, 37.72305, 52.58793, 43.47459))

s2 <- data.frame(lon= c(130.4327, 129.85127, 121.00775, 54.35377, 58.55841, 
94.95237),lat= c(30.25074, 29.83075, 23.7992, 31.61798, 52.58793, 43.47459))

sp1 <- SpatialPolygons(list(Polygons(list(Polygon(s1)),1)))
crs(sp1) <- CRS("+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0")

sp2 <- SpatialPolygons(list(Polygons(list(Polygon(s2)),1)))
crs(sp2) <-CRS ("+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0")

areasp1 <- areaPolygon(sp1)/10^6    # to get area in square km
areasp2 <- areaPolygon(sp2)/10^6

sp3 <- gSymdifference(sp1,sp2,checkvalidity = TRUE)
areasp3 <- areaPolygon(sp3)/10^6 

所有多边形都经过gIsValid() 的自相交或其他问题测试，因此与here 提到的问题无关。

知道为什么 sp1 的面积比 sp2 大吗？

注意：如果将areasp2 + areasp3 相加，则几乎等于areasp1。不知道是不是偶然。

Answer 1

这很好地说明了当您将角坐标视为平面坐标时，事物并不总是它们看起来的样子。您的绘图假定节点之间的直线连接。如果坐标是平面的，那将是合理的。但他们不是，在这种情况下，这一切都不同了。 下面我用原始数据重绘了绘图，并在节点之间的最短距离之后添加了顶点。您可以看到对象 2 比看起来要小得多。

library(raster)
library(geosphere)

s1 <- cbind(lon= c(130.4327, 129.85127, 121.00775, 84.9176, 60.40123, 58.97929, 58.55841, 94.95237),lat= c(30.25074, 29.83075, 23.7992, 28.25964, 36.89905, 37.72305, 52.58793, 43.47459))
s2 <- cbind(lon= c(130.4327, 129.85127, 121.00775, 54.35377, 58.55841, 94.95237),lat= c(30.25074, 29.83075, 23.7992, 31.61798, 52.58793, 43.47459))
sp1 <- spPolygons(s1, crs="+proj=longlat +datum=WGS84")
sp2 <- spPolygons(s2, crs="+proj=longlat +datum=WGS84")

# add nodes on shortest path
mp1 <- makePoly(sp1, interval=100000)
mp2 <- makePoly(sp2, interval=100000)

# background countries for map
library(maptools)
data(wrld_simpl)
w <- crop(wrld_simpl, extent(mp2) + 40)

plot(w, col='light gray')
points(s2, pch=20, cex=3)
points(s1, pch=20, cex=2, col='gray')

lines(sp2, col='red', lwd=3)
lines(sp1, col='blue', lwd=1)
lines(mp2, col='red', lwd=4, lty=2)
lines(mp1, col='blue', lwd=2, lty=2)

legend("bottomright", c('sp1', 'mp1', 'sp2', 'mp2'), col=rep(c('blue', 'red'), each=2), lwd=rep(c(2,3), each=2), lty=c(1, 2))

您可以通过将坐标投影到平面参考系统来更好地理解这一点。

library(rgdal)
crs <- CRS("+proj=ortho +lat_0=40 +lon_0=90 +a=6370997 +b=6370997 +units=m")
ww <- spTransform(w, crs)
sp1x <- spTransform(sp1, crs) 
sp2x <- spTransform(sp2, crs) 
mp1x <- spTransform(mp1, crs) 
mp2x <- spTransform(mp2, crs) 

par(mai=c(0,0,0,0))
plot(ww, col='light gray')
lines(sp2x, col='red', lwd=3)
lines(sp1x, col='blue', lwd=1)
lines(mp2x, col='red', lwd=4, lty=2)
lines(mp1x, col='blue', lwd=2, lty=2)

伊斯法罕到台北的最短路线不经过尼泊尔。