当经度大于90度时，Python计算纬度/经度中点的结果不正确

Question

我有一个问题，需要编写一个简短的函数来计算两个端点的经纬度时线段的中点。简单地说，当经度大于-90度或小于90度时，它可以正确地工作。对于另一半地球，它提供了一个有些随机的结果。

该代码是从http://www.movable-type.co.uk/scripts/latlong.html提供的javascript转换而来，并且似乎符合这里和这里的更正版本。与两个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)

任何建议？

- ichneumonad

---

放心吧：http://www.ig.utexas.edu/outreach/googleearth/latlong.html 也是错误的。 - S.Lott

@S.Lott：utexas代码有什么问题吗？ - John Machin

@John Machin：第二个例子给出的答案似乎不在两个点之间。 - S.Lott

@S.Lott：抱歉，我没有理解你所说的是你错误理解的例子。 - John Machin

@John Machin：正确。而且您明显地发现并纠正了我的错误。同时还问我的错误是什么。有点令人困惑。 - S.Lott

Answer 1

将您的参数设置代码替换为：

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)。请注意保留HTML标记。

- John Machin

---

更正为assert p1 [0]等，现在感觉自己像个白痴 - 是的，纬度和经度颠倒了。虽然只针对半个地球提供了正确答案！谢谢约翰。 - ichneumonad

@ichneumonad：请参考我的编辑版本，相比所有的p0[1]等内容，更具有更好的风格。 - John Machin

@ichneumonad：我有点难以相信交换纬度和经度会得到“半个地球”的相同答案。例如：midpoint((1, 2), (3, 4)) 产生 (2.0003044085023722, 2.999390393801055)，但 midpoint((2, 1), (4, 3)) 却产生 (3.0004561487854735, 1.9990851259125342)。 - John Machin

Answer 2

首先，很抱歉我要再留下一個答案。我有一個解決方案，可以解決在國際換日線兩側找到兩點中點的問題。我很想在現有的答案中添加評論，但我沒有足夠的聲望這樣做。

通過查看http://www.movable-type.co.uk/scripts/latlong.html工具的JavaScript文件，找到了解決方案。我正在使用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))

希望这个回答能够有所帮助，不会被视为不恰当，因为它是对之前回答中提出的问题的回答。

- carmenism

Answer 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

- Guiem Bosch