计算地球上两个坐标之间的距离
Calculate distance between two coordinates on a globe
我得到两个坐标对,形式为 90°0′0″N 0°0′0″E
作为字符串,我想计算半径为 R=6371km 的球体上这些点之间的距离。
我在网上找到了两个公式 here,"haversine" 和 "spherical law of cosines",但它们似乎不起作用。对于应该 return 2*pi*R / 4
的 90° 角,haversine 运行正确但余弦失败并且 return 0。具有更多随机坐标的不同点 returns 假值两种算法:haversine 太高,cosine 太低。
是我的实现有误还是我选择了错误的算法?
我应该如何进行这些计算(坐标对与地球表面的距离)?
(是的,我知道我还没有检查 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
您在 a
的计算中似乎混淆了 x
和 y
。你应该取纬度的余弦 (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))
我得到两个坐标对,形式为 90°0′0″N 0°0′0″E
作为字符串,我想计算半径为 R=6371km 的球体上这些点之间的距离。
我在网上找到了两个公式 here,"haversine" 和 "spherical law of cosines",但它们似乎不起作用。对于应该 return 2*pi*R / 4
的 90° 角,haversine 运行正确但余弦失败并且 return 0。具有更多随机坐标的不同点 returns 假值两种算法:haversine 太高,cosine 太低。
是我的实现有误还是我选择了错误的算法?
我应该如何进行这些计算(坐标对与地球表面的距离)?
(是的,我知道我还没有检查 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
您在 a
的计算中似乎混淆了 x
和 y
。你应该取纬度的余弦 (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))