很多人常常会有这样的疑问，知道两个点的经纬度，如：22:59:59.8888 113:59:59.8888 和 22:58:59.9999 113:59:59.9999， 怎么计算两者之间的距离。

这其实是测绘中的大地问题解算—大地问题反算，知道BL，求距离和方位角。

方法二比较简单，是网络上面找到的，出处已经不详，10米范围内和投影后的平面距离就差1cm。可用。

方法1是一个朋友写的，看起来比较专业。具体公式请参考武汉大学出版的《大地测量学》。

今天刚刚用到，所以和大家分享一下研究成果。

方法1：

/// <summary>
        /// 大地主题反算—已知2点经纬度，求距离和方位角. 84椭球。help by Wjj.
        /// </summary>
        /// <param name="latStart">起点纬度，单位：弧度</param>
        /// <param name="lonStart">起点经度，单位：弧度</param>
        /// <param name="latEnd">终点纬度，单位：弧度</param>
        /// <param name="lonEnd">终点经度，单位：弧度</param>
        /// <param name="dist">输出：大地线的距离，单位：m</param>
        /// <param name="angle">输出：大地线方位角，单位：弧度</param>
        private void GeodesicLine(double latStart, double lonStart, double latEnd, double lonEnd, ref double dist, ref double angle)
        {
            // 输入：170181690毫秒(47.272691666666667度)，475273606, 170181698, 475273509，结果是2.05米
            //WGS-84 ellipsphere
            double e2 = 0.0066943799013;
            double e = Math.Pow(e2 / (1 - e2), 0.5);
            double a = 6378137.0;

            double Bm = 0.5 * (latStart + latEnd);
            double deltaL = lonEnd - lonStart;
            double deltaB = latEnd - latStart;
            double t = Math.Tan(Bm);
            double p = e * Math.Cos(Bm);
            double W = Math.Pow((1 - e2 * Math.Pow(Math.Sin(Bm), 2)), 0.5);
            double Nm = a / W;
            double Vm = Math.Pow(W * (1 + e * e), 0.5);
            //
            double r01 = Nm * Math.Cos(Bm);
            double r21 = r01 * (1 - p * p - 9 * p * p * t * t) / 24;
            double r03 = Nm * Math.Pow(Math.Cos(Bm), 3) * t * t / 24;
            double S10 = Nm / (Vm * Vm);
            double S12 = Nm * Math.Pow(Math.Cos(Bm), 2) * (-2 - 3 * t * t + 3 * t * t * p * p) / 24;
            double S30 = Nm * (p * p - t * t * p * p) / 8;
            //
            double T1 = r01 * deltaL + r21 * deltaB * deltaB * deltaL + r03 * Math.Pow(deltaL, 3);
            double T2 = S10 * deltaB + S12 * deltaB * deltaL * deltaL + S30 * Math.Pow(deltaB, 3);
            //
            double t01 = t * Math.Cos(Bm);
            double t21 = Math.Cos(Bm) * t * (3 + 2 * p * p - 2 * Math.Pow(p, 4)) / 24;
            double t03 = Math.Pow(Math.Cos(Bm), 3) * t * (1 + p * p) / 12;
            double deltaA = t01 * deltaL + t21 * deltaB * deltaB * deltaL + t03 * Math.Pow(deltaL, 3);
            //
            double Am = Math.Atan2(T1, T2);
            dist = T2 / Math.Cos(Am); // 距离
            angle = Am - 0.5 * deltaA; // 角度，单位: 弧度
        }

使用：

double latStart = (1.0 * 170181690 * 0.001 / 3600) * Math.PI / 180.0;
                double lonStart = (1.0 * 475273606 * 0.001 / 3600) * Math.PI / 180.0;

                double latEnd = (1.0 * 170181698 * 0.001 / 3600) * Math.PI / 180.0;
                double lonEnd = (1.0 * 475273509 * 0.001 / 3600) * Math.PI / 180.0;
                double dist = 0;
                double angle = 0;
                GeodesicLine(latStart, lonStart, latEnd, lonEnd, ref dist, ref angle);
                this.Text = dist.ToString("F2") + " angle= " + (angle * 180.0 / Math.PI).ToString("F4");

方法2：(简单方法)

        /// <summary>
        /// 大地主题反算，知道两点的纬度和经度，求大地线的长度(和投影后的平面距离至差1cm左右). Help by LinXB
        /// </summary>
        /// <param name="p1">起点的纬度和经度，单位：度</param>
        /// <param name="p2">终点的纬度和经度，单位：度</param>
        /// <returns>输出：距离，单位：m</returns>
        private double CalcDistanceBy2LatLon(Coord p1, Coord p2)
        {
            // 输入：170181690毫秒，475273606, 170181698, 475273509，结果是2.05米
            double x, y, result;
            double R = 6378137.0; //84椭球

            x = (p2.lon - p1.lon) * Math.PI * R * Math.Cos(((p1.lat + p1.lat) / 2) * Math.PI / 180) / 180;
            y = Math.PI * (p2.lat - p1.lat) * R / 180;

            result = Math.Sqrt(x * x + y * y);
            return result;
        }

        /// <summary>
        /// 经纬度，单位：度
        /// </summary>
        public struct Coord
        {
            /// <summary>
            /// 纬度，单位：度
            /// </summary>
            public double lat;
            /// <summary>
            /// 经度，单位：度
            /// </summary>
            public double lon;
        }

使用：

Coord pt1 = new Coord();
                Coord pt2 = new Coord();
                pt1.lat = 1.0 * 170181690 / 1000 / 3600;
                pt1.lon = 1.0 * 475273606 / 1000 / 3600;
                pt2.lat = 1.0 * 170181698 / 1000 / 3600;
                pt2.lon = 1.0 * 475273509 / 1000 / 3600;
                double dist2 = CalcDistanceBy2LatLon(pt1, pt2);
                this.Text = dist2.ToString("F5");

两种方法在10米距离内，差别就1cm，和把两个点的经纬度用高斯投影变成平面坐标后求解的差1cm以内。

不过如果距离达到600米左右，方法1的结果大概是599米，方法2的结果大概是601米。

原因暂时没有分析。