2

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

我有一个包含另一个多边形 sp1 的多边形 sp2。所以 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()提到的问题无关。

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

注意：如果你求和areasp2+ areasp3，它几乎等于areasp1。我不知道这是不是偶然。

4

1 回答 1

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)

在此处输入图像描述

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

于 2018-02-07T04:16:19.423 回答