当经度 > 90 时，Python lat/long 中点计算给出错误的结果

Question

当给定两端点的纬度和经度时，我在计算直线中点的短函数时遇到了问题。简单的说就是当经度大于-90度或者小于90度的时候都可以正常工作。对于地球的另一半，它提供了一个有点随机的结果。

该代码是 http://www.movable-type.co.uk/scripts/latlong.html 提供的 javascript 的 python 转换, 并且似乎符合更正后的版本 here和 here .与两个 stackoverflow 版本进行比较时，我承认我不会用 C# 或 Java 编写代码，但我无法发现我的错误所在。

代码如下:

#!/usr/bin/python

import math

def midpoint(p1, p2):
   lat1, lat2 = math.radians(p1[0]), math.radians(p2[0])
   lon1, lon2 = math.radians(p1[1]), math.radians(p2[1])
   dlon = lon2 - lon1
   dx = math.cos(lat2) * math.cos(dlon)
   dy = math.cos(lat2) * math.sin(dlon)
   lat3 = math.atan2(math.sin(lat1) + math.sin(lat2), math.sqrt((math.cos(lat1) + dx) * (math.cos(lat1) + dx) + dy * dy))
   lon3 = lon1 + math.atan2(dy, math.cos(lat1) + dx)
   return(math.degrees(lat3), math.degrees(lon3))

p1 = (6.4, 45)
p2 = (7.3, 43.5)
print "Correct:", midpoint(p1, p2)

p1 = (95.5,41.4)
p2 = (96.3,41.8)
print "Wrong:", midpoint(p1, p2)

有什么建议吗？

Answer 1

将您的 arg 设置代码替换为:

lat1, lon1 = p1
lat2, lon2 = p2
assert -90 <= lat1 <= 90
assert -90 <= lat2 <= 90
assert -180 <= lon1 <= 180
assert -180 <= lon2 <= 180
lat1, lon1, lat2, lon2 = map(math.radians, (lat1, lon1, lat2, lon2))

然后再次运行您的代码。

更新关于涉及纬度/经度的计算的一些希望有用的一般建议:

1. 以度数或弧度输入纬度/经度？
2. 检查输入纬度/经度的有效范围
3. 检查输出纬度/经度的有效范围。经度在国际日期变更线上不连续。

中点例程的最后一部分可以有效地更改以避免长距离使用的潜在问题:

lon3 = lon1 + math.atan2(dy, math.cos(lat1) + dx)
# replacement code follows:
lon3d = math.degrees(lon3)
if lon3d < -180:
    print "oops1", lon3d
    lon3d += 360
elif lon3d > 180:
    print "oops2", lon3d
    lon3d -= 360
return(math.degrees(lat3), lon3d)

例如，找到新西兰奥克兰 (-36.9, 174.8) 和塔希提岛帕皮提 (-17.5, -149.5) 之间的中点会在通往有效答案的途中生成 oops2 194.270430902 (-28.355951246746923, -165.72956909809082)