如何计算两个坐标点的方位角

1.基于cesium计算方位角

/*
 版权： http://www.cesiumlab.com 免费的cesium基础数据处理平台
*/

function caculHeading(p0, p1) {

    //计算p0位置的enu位置矩阵的旋转部分
    var defrotmat = Cesium.Matrix4.getRotation(Cesium.Transforms.eastNorthUpToFixedFrame(p0), new Cesium.Matrix3());

    //获取矩阵的三个坐标轴
    var xaxis = Cesium.Matrix3.getColumn(defrotmat, 0, new Cesium.Cartesian3());
    var yaxis = Cesium.Matrix3.getColumn(defrotmat, 1, new Cesium.Cartesian3());
    var zaxis = Cesium.Matrix3.getColumn(defrotmat, 2, new Cesium.Cartesian3());

    //两个位置点的射线方向
    var dir = Cesium.Cartesian3.subtract(p1, p0,new Cesium.Cartesian3());

    //计算在 enu 旋转矩阵上的 x y 平面上的投影向量
    dir = Cesium.Cartesian3.cross(dir, zaxis, dir);
    dir = Cesium.Cartesian3.cross(zaxis, dir, dir);

    //归一化
    dir = Cesium.Cartesian3.normalize(dir, dir);

    //计算和x轴夹角  0 ~  pi
    var heading = Cesium.Cartesian3.angleBetween(xaxis, dir);

    //和y轴夹角 判定在y轴的正方向还是负方向
    var ay = Cesium.Cartesian3.angleBetween(yaxis, dir);

    // 保证处于0~2PI
    if (ay > Math.PI * 0.5) {
        heading = 2 * Math.PI - heading;
    }

    console.log("heading:" + Cesium.Math.toDegrees(heading).toFixed(5));
    return heading;
}

//测试

var p0 = Cesium.Cartesian3.fromDegrees(112.1, 30, 0);
var p1 = Cesium.Cartesian3.fromDegrees(112.1, 29, 0);

caculHeading(p0, p1);

2. 根据测量学计算公式计算

     /// <summary>
        /// 在已知经纬度的点上，沿bearing方向角，延长L距离。获得点
        /// </summary>
        /// <param name="lng"></param>
        /// <param name="lat"></param>
        /// <param name="L"></param>
        /// <param name="bearing"></param>
        /// <returns></returns>
        public static PointLatLng GetSecondPointLatlng(double lng, double lat, double L, double bearing)
        {
            double num = L / Axis;
            double num2 = Math.Acos(Math.Cos(1.5707963267948966 - DegreesToRadians(lat)) * Math.Cos(num) + Math.Sin(1.5707963267948966 - DegreesToRadians(lat)) * Math.Sin(num) * Math.Cos(DegreesToRadians(bearing)));
            double rad = Math.Asin(Math.Sin(num) * Math.Sin(DegreesToRadians(bearing)) / Math.Sin(num2));
            double lng2 = lng + RadiansToDegrees(rad);
            double lat2 = 90.0 - RadiansToDegrees(num2);
            return new PointLatLng(lng2, lat2);
        }
        /// <summary>
        /// 计算方位角,P1P2
        /// </summary>
        /// <param name="p1"></param>
        /// <param name="p2"></param>
        /// <returns></returns>
        public static double GetBearing(PointLatLng p1, PointLatLng p2)
        {
            double num = DegreesToRadians(p1.Lat);
            double num2 = DegreesToRadians(p2.Lat);
            double num3 = DegreesToRadians(p2.Lng - p1.Lng);
            double y = Math.Sin(num3) * Math.Cos(num2);
            double x = Math.Cos(num) * Math.Sin(num2) - Math.Sin(num) * Math.Cos(num2) * Math.Cos(num3);
            return (RadiansToDegrees(Math.Atan2(y, x)) + 360.0) % 360.0;
        }

 public static double RadiansToDegrees(double rad)
        {
            return 57.295779513082323 * rad;
        }