两点之间的方位角

Question

我一直在使用geopy包，它表现得非常出色，但有些结果是不一致的或者存在较大的偏差，我怀疑问题出在我的方位角计算上：

def gb(x,y,center_x,center_y):
dx=x-center_x
dy=y-center_y
if ((dy>=0)and((dx>0)or(dx<0))):
    return math.degrees(math.atan2(dy,dx))
elif (dy<=0)and((dx>0)or (dx<0)):
    return (math.degrees(math.atan2(dy,dx))+360)
else:
    return (math.degrees(math.atan2(dy,dx))+360)%360

我需要计算方位角，使中心点(center_x和center_y)作为支点。之后我将使用geopy来反向推导GPS坐标:

latlon = VincentyDistance(miles=dist).destination(Point(lat1, lon1), bearing)

有没有人能指出我可能做错了什么？

- 242Eld

---

如果您的问题已经得到回答解决，您应该点击答案旁边的大“勾”接受答案。否则，请提供反馈，以便我们可以尝试进一步帮助您。 - John Machin

Answer 1

有人能指出我哪里做错了吗？

1. 没有展示“不一致或存在相对较大位移”的实际结果或你的预期结果，因此回答者必须依靠猜测。

2. 没有说明输入（x、y等）的单位以及如何获得用于“destination”计算中的“dist”。在计算下面的“bearing2”时，我假设正向x是以英里为单位的东部，在英里中正向y是北部。如果您编辑问题以修复（1）和（2），将会有很大帮助。

3. 编码风格不太有利于让人们阅读...请查看这个。

4. 在学校的三角学中，角度是逆时针从X轴（东）开始度量的。在导航中，方位角是从Y轴（北）开始顺时针测量的。请参见下面的代码。关于航向的例子，请点击此链接，滚动到“从起始点给定距离和方位角的目标点”部分，注意该示例讨论的方位角大约为96或97度，然后单击“查看地图”，您会注意到航向稍微偏向东南（东为90度）。

代码：

from math import degrees, atan2
    def gb(x, y, center_x, center_y):
        angle = degrees(atan2(y - center_y, x - center_x))
        bearing1 = (angle + 360) % 360
        bearing2 = (90 - angle) % 360
        print "gb: x=%2d y=%2d angle=%6.1f bearing1=%5.1f bearing2=%5.1f" % (x, y, angle, bearing1, bearing2)

    for pt in ((0, 1),(1,1),(1,0),(1,-1),(0,-1),(-1,-1),(-1, 0),(-1,1)):
        gb(pt[0], pt[1], 0, 0)

输出：

gb: x= 0 y= 1 angle=  90.0 bearing1= 90.0 bearing2=  0.0
gb: x= 1 y= 1 angle=  45.0 bearing1= 45.0 bearing2= 45.0
gb: x= 1 y= 0 angle=   0.0 bearing1=  0.0 bearing2= 90.0
gb: x= 1 y=-1 angle= -45.0 bearing1=315.0 bearing2=135.0
gb: x= 0 y=-1 angle= -90.0 bearing1=270.0 bearing2=180.0
gb: x=-1 y=-1 angle=-135.0 bearing1=225.0 bearing2=225.0
gb: x=-1 y= 0 angle= 180.0 bearing1=180.0 bearing2=270.0
gb: x=-1 y= 1 angle= 135.0 bearing1=135.0 bearing2=315.0

- John Machin

---

你好 John，关于我的编码风格，你说得对，谢谢你提醒我。至于代码，我相信你提供的实现方式更好，我获取方位角的代码只是从维基百科上找来的。我会再次运行我的脚本并测试结果，谢谢！ - 242Eld

@242Eld：请提供您提到的维基百科文章的URL。另外，您还没有回答我的答案中的第（1）和第（2）点。 - John Machin

维基百科的文章是希伯来语。由于我正在进行反向地理编码，因此无法使用代码描述示例。但问题领域是，在从必应地图解析HTML页面后，我需要反向地理编码过程，以便检索所需的多边形。结果很好，但不令人满意，因为多边形与原始位置相比存在200-1500米的位移。我使用像素中的x y坐标，将其转换为英里，使用给定的比例尺。反向地理编码是使用geopy包的VincentyDistance完成的。 - 242Eld

这并不是回答原问题，但是 john-machin 要求提供包含数学知识的维基站点，此站点 提到了和 242Eld 使用相同的数学知识，如果有人需要可以去查看。 - Otzen

Answer 2

我不太确定你在代码中想要做什么，但我看到一些可能需要清理的奇怪现象。

1. 在条件语句之前，你先测试了 dy>=0，然后测试了 dy<=0。如果 dy==0 且 dx==0 ，你的代码应该做什么？
2. 你的测试 ((dy>=0)and((dx>0)or(dx<0))) 相当于 (dy>=0 and dx!=0)，这是你想要的吗？
3. 在所有条件语句中，你基本上都在做同样的事情。在这种情况下，return math.degrees(math.atan2(dy,dx))+360)%360 是否可以在每种情况下都起作用？如果是这样，你就不需要使用 if 语句了。

- amccormack

---

+1 特别是针对3，但它可能不能回答问题。您知道例如 math.atan2 是否在 dx==0 附近给出相对不准确的结果（它是否包含一些逻辑来避免除以接近零的问题）？如果需要 if 语句，则可能是为了围绕 45 度线反射，以获得更好的数值精度 - 但我认为库应该已经处理了这个问题。 - user180247

@Steve314，说得好。我查了Python2.7文档，它没有关于输入0的警告。此外，我尝试了math.atan2(0,30)和math.atan2(30,0)，都没有返回错误，如果您给它无效的输入，我会期望内置库返回错误。 - amccormack