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

Posted 2023-02-23

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

【英文标题】：Python lat/long midpoint calculation gives wrong result when longitude > 90

【发布时间】：2011-05-05 10:00:31

【问题描述】：

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

该代码是http://www.movable-type.co.uk/scripts/latlong.html 提供的javascript 的python 转换，并且似乎符合here 和here 的更正版本。在比较两个 *** 版本时，我承认我不会用 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)

有什么建议吗？

【问题讨论】：

安慰一下：ig.utexas.edu/outreach/googleearth/latlong.html 也是错的。

@S.Lott：utexas 代码有什么问题？

@S.Lott：短距离有什么相同的错误？？

@John Machin：第二个示例提供的答案似乎不在两点之间。

@S.Lott：如果您的意思是 OP 的第二个示例 p1=(95.5,41.4) 等：纬度的 95.5 度是无效的——它将比北极更北。阅读我的答案- OP 得到了 lat 和 lon 颠倒并承认了。如果你的意思是别的，请解释一下。

【参考方案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))

然后再次运行您的代码。

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

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

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

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)

【讨论】：

更正为断言 p1[0] 等，现在我觉得自己像个白痴 - 是的，纬度和经度是错误的。尽管如此，仍然为半个地球提供了正确的答案！谢谢约翰。

@ichneumonad：参考我的编辑版本，以获得比所有 p0[1] 东西更好的建议

@ichneumonad：我很难相信交换 lat 和 lon 会为“半个星球”提供相同的答案。例如：midpoint((1, 2), (3, 4)) 产生 (2.0003044085023722, 2.999390393801055) 但 midpoint((2, 1), (4, 3)) 产生 (3.0004561487854735, 1.9990851259125342)

【参考方案2】：

首先，很抱歉我要留下另一个答案。对于该答案中提到的问题，我有一个解决方案，涉及如何找到日期线两侧两点的中点。我很想简单地对现有答案添加评论，但我没有这样做的声誉。

通过查看支持该工具的 Javascript 文件找到了解决方案，地址为 http://www.movable-type.co.uk/scripts/latlong.html。我正在使用 shapely.geometry.Point。如果您不想安装此软件包，则使用元组代替也可以。

def midpoint(pointA, pointB):
    lonA = math.radians(pointA.x)
    lonB = math.radians(pointB.x)
    latA = math.radians(pointA.y)
    latB = math.radians(pointB.y)

    dLon = lonB - lonA

    Bx = math.cos(latB) * math.cos(dLon)
    By = math.cos(latB) * math.sin(dLon)

    latC = math.atan2(math.sin(latA) + math.sin(latB),
                  math.sqrt((math.cos(latA) + Bx) * (math.cos(latA) + Bx) + By * By))
    lonC = lonA + math.atan2(By, math.cos(latA) + Bx)
    lonC = (lonC + 3 * math.pi) % (2 * math.pi) - math.pi

    return Point(math.degrees(lonC), math.degrees(latC))

我希望这是有帮助的，并且不会被认为是不恰当的，因为它是对上一个答案中提出的问题的答案。

【参考方案3】：

不是直接回答 >90 的问题，但它对我的情况有所帮助，所以我把它留在这里以防万一它可以帮助其他人。

我只需要将中点放在一张地图中，以度数为单位。我没有使用弧度转换，并且通过使用简单的欧几里得距离在地图中正确描绘了它。示例函数：

def midpoint_euclidean(x1,y1,x2,y2):
    dist_x = abs(x1-x2) / 2.
    dist_y = abs(y1-y2) / 2.
    res_x = x1 - dist_x if x1 > x2 else x2 - dist_x
    res_y = y1 - dist_y if y1 > y2 else y2 - dist_y
    return res_x, res_y