如何计算两个WGS84坐标之间的方位角（相对于北方的角度）？

Question

我有两个WGS84坐标，纬度和经度以度为单位。这些点非常接近，例如仅相隔一米。是否有一种简单的方法来计算连接这些点的方位角，也就是与北方的夹角？幼稚的方法是假设一个笛卡尔坐标系（因为这些点非常接近）并且只使用

sin(a) = abs(L2-L1) / sqrt(sqr(L2-L1) + sqr(B2-B1))

a = 方位角 L1，L2 = 经度 B1，B2 = 纬度

误差会随着坐标远离赤道而变大，因为在那里，两个经度度之间的距离越来越小，而两个纬度度之间的距离保持不变。

我发现了一些非常复杂的公式，我真的不想实施它们，因为它们似乎对于非常接近的点来说是过度的，并且我不需要非常高的精度（两个小数点足以，一个可能已经足够了，因为还存在其他降低精度的因素，例如GPS返回的因素）。

也许我可以根据纬度确定一个近似的经度修正因子，然后使用类似于以下内容的东西：

sin(a) = abs(L2*f-L1*f) / sqrt(sqr(L2*f-L1*f) + sqr(B2-B1))

其中f是修正因子

有什么提示吗？

（我不想使用任何库来完成此操作，特别是那些需要运行时许可证的库。任何MPLed Delphi源代码都很好。）

- dummzeuch

---

说句闲话，你要找的术语是“标题”。 - hobbs

Answer 1

你在文本中提到的公式是用于计算两点之间的大圆距离。以下是我计算两点间夹角的方法：

uses Math, ...;
...

const
  cNO_ANGLE=-999;

...

function getAngleBetweenPoints(X1,Y1,X2,Y2:double):double;
var
  dx,dy:double;
begin
  dx := X2 - X1;
  dy := Y2 - Y1;

  if (dx > 0) then  result := (Pi*0.5) - ArcTan(dy/dx)   else
  if (dx < 0) then  result := (Pi*1.5) - ArcTan(dy/dx)   else
  if (dy > 0) then  result := 0                          else
  if (dy < 0) then  result := Pi                         else
                    result := cNO_ANGLE; // the 2 points are equal

  result := RadToDeg(result);
end;

• 记得处理两个点重合的情况（检查结果是否等于cNO_ANGLE，或修改函数以抛出异常）；

• 此函数假定您在平坦的表面上。对于你提到的小距离来说这没问题，但如果你要计算全球各地城市之间的方向，你可能需要考虑一些将地球形状考虑在内的方法；

• 最好提供已经映射到平面的坐标给这个函数。您可以直接将WGS84纬度输入Y（经度输入X），以获得粗略的近似值。

- Wouter van Nifterick

---

太棒了！我将它翻译成了C#来替换我的先前答案。 - Jader Dias

@Martin Beckett：是的，atan2也可以用Atan2(dy;dx)来实现这个功能，但是当然，负数dx会返回一个负值，需要加上360。而且0也需要单独处理。 - MPelletier

Answer 2

这是C#的解决方案。经过了0度、45度、90度、135度、180度、225度、270度和315度的测试。

编辑：我用Wouter的解决方案翻译成了C#，替换了以前丑陋的解决方案。

public double GetAzimuth(LatLng destination)
{
    var longitudinalDifference = destination.Lng - this.Lng;
    var latitudinalDifference = destination.Lat - this.Lat;
    var azimuth = (Math.PI * .5d) - Math.Atan(latitudinalDifference / longitudinalDifference);
    if (longitudinalDifference > 0) return azimuth;
    else if (longitudinalDifference < 0) return azimuth + Math.PI;
    else if (latitudinalDifference < 0) return Math.PI;
    return 0d;
}

public double GetDegreesAzimuth(LatLng destination)
{
    return RadiansToDegreesConversionFactor * GetAzimuth(destination);
}

- Jader Dias

---

最好使用Math.Atan2而不是Math.Atan，以避免除法，特别是避免除以零。另一个优点是它给出的答案范围为-π..+π，从而避免了您在最后进行的更正。 - Jonathan Leffler

Answer 3

我找到了这个链接：http://williams.best.vwh.net/avform.htm，它在回答中给出了Lat/Lon + Distance + Heading --> Lat/Lon。这看起来非常有前途，尤其是在末尾给出的平面地球近似解法。

- dummzeuch

Answer 4

这只适用于小差异。否则，您不能仅使用“纬度差/经度差”。

- Yosi

Answer 5

有人测试过这个吗？它没有返回正确的答案。

此函数假定您在平坦的表面上。对于您提到的小距离，这一切都很好，但如果您要计算世界各地城市之间的航向，您可能需要考虑一些考虑了地球形状的东西；

您的平坦地球与此无关。正如您所说的那样，错误是因为您正在从一个点计算初始方位角。除非您直接朝着极点前进，否则您与极点的关系将随着距离的增加而改变。无论如何，上述程序都不会返回正确的结果。

- Meta Mussel

Answer 6

我建议根据经度实现一个修正因子。我曾经实现过类似的例程，以返回特定位置x英里内的所有地理编码记录，并遇到了类似的问题。不幸的是，我现在没有代码了，也想不起来我是如何得到修正数字的，但你正在正确的轨道上。

- skamradt