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

当给定每一端点的纬度和经度时，我有一个短函数来计算线的中点的问题。简单来说就是在经度大于-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)

有什么建议么？

最佳答案

将您的 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))

并再次运行您的代码。

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

以度或弧度为单位输入纬度/经度？

检查输入纬度/经度的有效范围

检查 OUTPUT lat/lon 的有效范围。经度在国际日期变更线上不连续。

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

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)

关于当经度 > 90 时，Python 纬度/经度中点计算给出错误结果，我们在Stack Overflow上找到一个类似的问题：https://stackoverflow.com/questions/5895832/