delphi - 如何计算两个 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 源代码都很棒。)

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;

• 记得处理2个点相等的情况(检查结果是否等于cNO_ANGLE，或者修改函数抛出异常)；

• 此函数假设您位于平坦的表面上。对于您提到的小距离来说，这一切都很好，但如果您要计算世界各地城市之间的航向，您可能需要研究一些采用地球形状进行计数的东西；

• 最好为该函数提供已映射到平面的坐标。不过，您可以将 WGS84 Latitude 直接输入 Y(将 lon 输入 X)以获得粗略的近似值。