Python 经度大于 90 时的经纬度中点计算错误

4 投票

3 回答

3429 浏览

提问于 2025-04-16 17:02

我遇到了一个问题，想写一个简单的函数来计算一条线的中点，前提是我有这条线两端的经纬度。简单来说，当经度大于-90度或小于90度时，这个函数能正常工作。但在地球的另一半，它给出的结果就有点随机了。

这段代码是把一个在 http://www.movable-type.co.uk/scripts/latlong.html 上的JavaScript转换成Python的，看起来也符合在这里和这里的修正版本。虽然我承认我不太会用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)

有什么建议吗？

计算机科学编程错误数学函数经纬度转换经纬度地理坐标中点计算地球模型

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

回答于 2025-04-16 由 Python大师

分享举报

首先，抱歉我又来回答一个问题。我找到了一个解决方案，关于如何找到国际日期变更线两侧两个点的中点。我本来想在已有的回答下评论，但我没有足够的声望来这样做。

这个解决方案是通过查看这个工具的Javascript文件找到的，网址是 http://www.movable-type.co.uk/scripts/latlong.html。我使用的是 shapely.geometry.Point。如果你不想安装这个包，使用元组（tuples）也能达到同样的效果。

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))

希望这个信息对你有帮助，也希望大家不要觉得这不合适，因为这是对之前回答中提到的问题的解答。

回答于 2025-04-16 由 Python大师

分享举报

把你的参数设置代码替换成：

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) 的过程中出现了问题。

回答于 2025-04-16 由 Python大师

分享举报

Python 经度大于 90 时的经纬度中点计算错误

3 个回答

撰写回答