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

Question

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

该代码是http://www.movable-type.co.uk/scripts/latlong.html上提供的 JavaScript 的 Python 转换，并且似乎符合此处和此处的更正版本。在与两个 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)