【问题标题】:Calculate distance between two coordinates on a globe计算地球上两个坐标之间的距离
【发布时间】:2016-07-16 19:29:15
【问题描述】:

我得到两个90°0′0″N 0°0′0″E 形式的坐标对作为字符串,并想计算半径为 R=6371km 的球体上这些点之间的距离。

我在互联网here 上找到了两个公式,“haversine”和“sphere of cosines”,但它们似乎不起作用。对于应返回 2*pi*R / 4 的 90° 角,harsine 操作正确,但余弦失败并返回 0。具有更多随机坐标的不同点在两种算法中都返回错误值:harsine 太高而余弦太低。

是我的实现错误还是我选择了不正确的算法?

我应该如何进行这些计算(坐标对到地球表面上的距离)?

(是的,我知道我还没有检查 N/S 和 E/W,但测试的坐标都在东北半球。)

这是我的 Python 3 代码:

import math, re
R = 6371
PAT = r'(\d+)°(\d+)′(\d+)″([NSEW])'

def distance(first, second):
    def coords_to_rads(s):  
        return [math.radians(int(d) +int(m)/60 +int(s)/3600) \
                for d, m, s, nswe in re.findall(PAT, s)]

    y1, x1 = coords_to_rads(first)
    y2, x2 = coords_to_rads(second)  
    dx = x1 - x2  
    dy = y1 - y2  

    print("coord string:", first, "|", second)
    print("coord radians:", y1, x1, "|", y2, x2)
    print("x/y-distances:", dy, dx)

    a = math.sin(dx/2)**2 + math.cos(x1) * math.cos(x2) * math.sin(dy/2)**2  
    c = 2 * math.atan2(math.sqrt(a), math.sqrt(1-a))  
    haversine =  R * c  

    law_of_cosines = math.acos( math.sin(x1) * math.sin(x2) + \
                                math.cos(x1) * math.cos(x2) ) * R

    print("HS:", round(haversine, 2), "LOC:", round(law_of_cosines, 2))

    return haversine
    #return law_of_cosines

if __name__ == '__main__':
    def test(result, correct):
        print("result: ", result)
        print("correct:", correct)

    test(distance("90°0′0″N 0°0′0″E", "0°0′0″N, 0°0′0″E"), 10007.5)
    test(distance("51°28′48″N 0°0′0″E", "46°12′0″N, 6°9′0″E"), 739.2)
    test(distance("90°0′0″N 0°0′0″E", "90°0′0″S, 0°0′0″W"), 20015.1)
    test(distance("33°51′31″S, 151°12′51″E", "40°46′22″N 73°59′3″W"), 15990.2)

这是一些输出:

coord string: 90°0′0″N 0°0′0″E | 0°0′0″N, 0°0′0″E
coord radians: 1.5707963267948966 0.0 | 0.0 0.0
x/y-distances: 1.5707963267948966 0.0
HS: 10007.54 LOC: 0.0
result: 10007.543398010286
correct: 10007.5

coord string: 51°28′48″N 0°0′0″E | 46°12′0″N, 6°9′0″E
coord radians: 0.8984954989266809 0.0 | 0.8063421144213803 0.10733774899765128
x/y-distances: 0.09215338450530064 -0.10733774899765128
HS: 900.57 LOC: 683.85
result: 900.5669567853056
correct: 739.2

【问题讨论】:

  • 好问题,但我想知道它是否应该转到Code Review?我知道很多问题,其中你已经有一个工作或几乎工作的代码示例被认为是离题这里和那里的主题,但我不确定网站之间的确切“截止”。
  • @DavidZ 据我了解,CR 只想要寻找优化的工作代码 sn-ps。所以我把它贴在这里。但这不会是我第一个从 SO 迁移到 CR 或返回...
  • 对我来说很有意义。无论如何,我检查了您与Wikipedia 链接的站点,似乎在haversine 公式的定义中存在不一致(尽管我不确定它是否真的不一致)。你调查过吗?
  • @DavidZ On-Topic 的 3 条黄金法则Code Review 问题:代码是否按预期工作?它是实际代码(不是存根、示例或假设)吗? OP 是否乐于收到有关所述代码任何方面的反馈?如果所有 3 点都是,那么 CR 就是去处。
  • 在这种情况下,代码显然由于“未按预期工作”而失败。

标签: python coordinates distance coordinate-systems


【解决方案1】:

您在计算a 时似乎混淆了xy。您应该取纬度的余弦 (y),而不是经度 (x)。

我通过将您的 distance 更改为 angular_distance(即不要乘以 R)并添加一些额外的测试来发现这一点:

test(angular_distance("90°0′0″N 0°0′0″E", "89°0′0″N, 0°0′0″E"), math.radians(1))
test(angular_distance("90°0′0″N 0°0′0″E", "80°0′0″N, 0°0′0″E"), math.radians(10))
test(angular_distance("90°0′0″N 0°0′0″E", "50°0′0″N, 0°0′0″E"), math.radians(40))
test(angular_distance("90°0′0″N 0°0′0″E", "50°0′0″N, 20°0′0″E"), math.radians(40))

【讨论】:

  • 谢谢,你是对的。我只是更改了y1, x1 = coords_to_rads(first)y2, x2 = coords_to_rads(second) 行并交换了xs 和ys。