大圆计算的差异答案

【问题标题】：Difference in Great Circle calculation大圆计算的差异
【发布时间】：2016-08-23 09:38:22
【问题描述】：

我在使用 geosphere 包中的 perimeter 函数时遇到了一些问题。

根据第一原理，1 度弧等于 1 纳米所以弧度 90° = 90 * 60 = 5400nm

当我使用 distGeo 和赤道上的两个点时，我得到的答案与该图相对应。

library(geosphere)
distGeo(c(0,0),c(0,90), a=6378137, f=1/298.257223563)/1852
#[1] 5400.629

如果我尝试用周长做同样的计算，我会得到两倍的距离！

#Two points on equator
xy <- rbind(c(0,0), c(90,0))     
xy
#     [,1] [,2]
#[1,]    0    0
#[2,]   90    0
perimeter(xy)/1852
#[1] 10819.39
perimeter(xy, a=6378137, f=1/298.257223563)/1852
#[1] 10819.39

如果这次我们沿着本初子午线在赤道和北极之间取另一个点，结果是......

distGeo(c(0,0),c(90,0), a=6378137, f=1/298.257223563)/1852
#[1] 5409.694

这显示了椭球的效果。

当我对这些新点使用周长时，我仍然遇到和以前一样的问题！

#From the equator to the north Pole along the Prime meridian
xy <- rbind(c(0,0), c(0,90)) 
xy
#     [,1] [,2]
#[1,]    0    0
#[2,]    0   90
perimeter(xy)/1852
#[1] 10801.26
perimeter(xy, a=6378137, f=1/298.257223563)/1852
#[1] 10801.26

那我做错了什么？

史蒂夫

【问题讨论】：

标签： r

【解决方案1】：

你没有做错任何事。函数perimeter 计算封闭多边形周围的周长距离，如果端点不是第一个，它将关闭该多边形，因此它将您的两个点解释为没有“高度”的多边形。

为了说明（结果以米为单位，可以看到小的添加）：

## your results
xy <- rbind(c(0,0), c(90,0))
perimeter(xy)
##[1] 20037508

## "closing the polygon" with c(0,0) gives same answer
xy <- rbind(c(0,0), c(90,0), c(0,0)) 
perimeter(xy)
##[1] 20037508

## adding one meter (0.00001 deg latitude) north as another point, 
## perimeter will close the resulting triangle, result increases by one meter!
xy <- rbind(c(0,0), c(90,0), c(90, 0.00001))
perimeter(xy)
##[1] 20037509

## closing the triangle yourself, same answer
xy <- rbind(c(0,0), c(90,0), c(90,0.00001), c(0,0)) 
perimeter(xy)
##[1] 20037509

【讨论】：

【解决方案2】：

感谢您的指导。

我最终做的是滞后于经纬度字段

xy[,3] = lag(xy[,1],k=-1)
xy[,4] = lag(xy[,2],k=-1)

然后使用 distGEO 计算分段

xy$dDist<-distGeo(as.vector(xy[,1:2]), 
    as.vector(xy[,3:4]),
    a=6378137,
    f=1/298.257223563)/1852

然后取总和

sum(xy$dDist,na.rm = TRUE)

它似乎工作！

【讨论】：