怎么用Python计算地球上任意两个用经纬度表示的点的弧面距离?

举报
天元浪子 发表于 2021/07/26 22:54:14 2021/07/26
【摘要】 这是来自知乎上的问题。问我的时候,恰好我在写一个和向量计算相关的文章,于是灵光乍现,顺手写了这样一个答案。该算法未经严格验证,请谨慎参考。具体思路如下。 将两个点的经纬度换算成空间坐标;计算地心与两个点所成的两个向量的点积;点积除以两个向量的模(也就是地球半径)之积,结果就是向量夹角的余弦;反余弦值对应着两点所在大圆(即经过两点的地球表面最大的圆)的弧度;弧度乘以地球...

这是来自知乎上的问题。问我的时候,恰好我在写一个和向量计算相关的文章,于是灵光乍现,顺手写了这样一个答案。该算法未经严格验证,请谨慎参考。具体思路如下。

  1. 将两个点的经纬度换算成空间坐标;
  2. 计算地心与两个点所成的两个向量的点积;
  3. 点积除以两个向量的模(也就是地球半径)之积,结果就是向量夹角的余弦;
  4. 反余弦值对应着两点所在大圆(即经过两点的地球表面最大的圆)的弧度;
  5. 弧度乘以地球半径,即得弧长。

代码如下。

>>> def get_arc(p0, p1, r=1):
	z0 = r*np.sin(np.radians(p0[1]))
	x0 = r*np.cos(np.radians(p0[1]))*np.cos(np.radians(p0[0]))
	y0 = r*np.cos(np.radians(p0[1]))*np.sin(np.radians(p0[0]))
	z1 = r*np.sin(np.radians(p1[1]))
	x1 = r*np.cos(np.radians(p1[1]))*np.cos(np.radians(p1[0]))
	y1 = r*np.cos(np.radians(p1[1]))*np.sin(np.radians(p1[0]))
	theta = np.arccos(np.dot((x0,y0,z0),(x1,y1,z1))/(r*r))
	return theta * r

>>> r = 6377.830 # 使用赤道半径,单位:km
>>> 北京 = (116.5,40.0)
>>> 济南 = (117.0,36.4)
>>> 悉尼 = (151.2,-33.9)
>>> 纽约 = (-74.0,40.5)
>>> 巴西利亚 = (-48.0,-15.9)
>>> get_arc(北京, 济南, r=r)
403.10853123743505
>>> get_arc(北京, 悉尼, r=r)
8966.07982007138
>>> get_arc(北京, 纽约, r=r)
11012.732277406261
>>> get_arc(北京, 巴西利亚, r=r)
16962.014473669606

  
 
  • 1
  • 2
  • 3
  • 4
  • 5
  • 6
  • 7
  • 8
  • 9
  • 10
  • 11
  • 12
  • 13
  • 14
  • 15
  • 16
  • 17
  • 18
  • 19
  • 20
  • 21
  • 22
  • 23
  • 24

文章来源: xufive.blog.csdn.net,作者:天元浪子,版权归原作者所有,如需转载,请联系作者。

原文链接:xufive.blog.csdn.net/article/details/114881245

【版权声明】本文为华为云社区用户转载文章,如果您发现本社区中有涉嫌抄袭的内容,欢迎发送邮件进行举报,并提供相关证据,一经查实,本社区将立刻删除涉嫌侵权内容,举报邮箱: cloudbbs@huaweicloud.com
  • 点赞
  • 收藏
  • 关注作者

评论(0

0/1000
抱歉,系统识别当前为高风险访问,暂不支持该操作

全部回复

上滑加载中

设置昵称

在此一键设置昵称,即可参与社区互动!

*长度不超过10个汉字或20个英文字符,设置后3个月内不可修改。

*长度不超过10个汉字或20个英文字符,设置后3个月内不可修改。