geosphere::areaPolygon 的面积计算错误？
Wrong area calculation with geosphere:: areaPolygon?
我已经成功使用 geosphere::areaPolygon() 多次，但现在我遇到了一个问题。
我有一个包含另一个多边形 sp1 的多边形 sp2。所以 sp2 的面积应该比 sp1 大。当我用 areaPolygon() 计算每个区域时，我得到相反的结果。
areasp1 = 10133977<br>
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。不知道是不是偶然
这很好地说明了当您将 angular 坐标视为平面坐标时，事情并不总是它们看起来的样子。您的情节假定节点之间存在直线连接。如果坐标是平面的，那将是合理的。但事实并非如此，这在这种情况下产生了重大影响。
下面我用原始数据重绘了绘图，并在节点之间的最短距离后添加了顶点。您可以看到对象 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)
伊斯法罕到台北的最短路线不经过尼泊尔。
我已经成功使用 geosphere::areaPolygon() 多次，但现在我遇到了一个问题。
我有一个包含另一个多边形 sp1 的多边形 sp2。所以 sp2 的面积应该比 sp1 大。当我用 areaPolygon() 计算每个区域时，我得到相反的结果。
areasp1 = 10133977<br>
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。不知道是不是偶然
这很好地说明了当您将 angular 坐标视为平面坐标时，事情并不总是它们看起来的样子。您的情节假定节点之间存在直线连接。如果坐标是平面的，那将是合理的。但事实并非如此，这在这种情况下产生了重大影响。 下面我用原始数据重绘了绘图，并在节点之间的最短距离后添加了顶点。您可以看到对象 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)
伊斯法罕到台北的最短路线不经过尼泊尔。