空间两点间距离计算

　　原来使用超图软件进行道路长度的统计，后来决定在项目中使用开源的工具，两点之间距离量算的方法，网上很容易找到，基本上是把地球作为一个标准球体，使用固定半径来进行计算，这样的每一个计算结果与实际长度差距不是很大，一般是可以接受的，但是计算对象较多时统计的结果与超图的计算结果对比就会发现有较大的差距，这是因为超图的长度计算中地球是作为椭球体处理的。

　　下面提供一种椭球体表面两点间的距离量算，“Vincenty's formulae”，代码非原创，有趣的朋友可以去原作者页面研究一下，源码下载页面 http://www.codeproject.com/Articles/19939/GPS-Receivers-Geodesy-and-Geocaching-Vincenty-s-Fo

WGS84坐标系下，地球赤道半径、极半径、平均半径长度及椭球体扁率

        private static  readonly double _radiusA = 6378137;
        private static  readonly double _radiusB = 6356752.314;
        private static readonly double _radius = 6371000.0;
        private static  readonly double _flatten = 1 / 298.257223563;

        private static double Deg2Rad(double d)
        {
            return d * Math.PI / 180;
        }

使用固定地球半径进行距离计算的函数

　　　　 // Greater Circle Distance
        public static double Distance(double long1, double lat1, double long2, double lat2)
        {
            double dLat = (lat2 - lat1) * (Math.PI / 180);
            double dLon = (long2 - long1) * (Math.PI / 180);
            double a = Math.Sin(dLat / 2) * Math.Sin(dLat / 2) +
                    Math.Cos(lat1 * (Math.PI / 180)) * Math.Cos(lat2 * (Math.PI / 180)) *
                    Math.Sin(dLon / 2) * Math.Sin(dLon / 2);
            double c = 2 * Math.Atan2(Math.Sqrt(a), Math.Sqrt(1 - a));
            return _radiusA * c;
        }

Vincenty's formulae

　　　　 // Vincenty Formula Distance
        public static double VincentyDistance(double lat1, double lon1, double lat2, double lon2)
        {

            double phi1 = Deg2Rad(lat1);
            double lambda1 = Deg2Rad(lon1);
            double phi2 = Deg2Rad(lat2);
            double lambda2 = Deg2Rad(lon2);

            double a2 = _radiusA * _radiusA;
            double b2 = _radiusB * _radiusB;
            double a2b2b2 = (a2 - b2) / b2;

            double omega = lambda2 - lambda1;

            double tanphi1 = Math.Tan(phi1);
            double tanU1 = (1.0 - _flatten) * tanphi1;
            double U1 = Math.Atan(tanU1);
            double sinU1 = Math.Sin(U1);
            double cosU1 = Math.Cos(U1);

            double tanphi2 = Math.Tan(phi2);
            double tanU2 = (1.0 - _flatten) * tanphi2;
            double U2 = Math.Atan(tanU2);
            double sinU2 = Math.Sin(U2);
            double cosU2 = Math.Cos(U2);

            double sinU1sinU2 = sinU1 * sinU2;
            double cosU1sinU2 = cosU1 * sinU2;
            double sinU1cosU2 = sinU1 * cosU2;
            double cosU1cosU2 = cosU1 * cosU2;

            double lambda = omega;

            double A = 0.0;
            double B = 0.0;
            double sigma = 0.0;
            double deltasigma = 0.0;
            double lambda0;
            bool converged = false;

            for (int i = 0; i < 20; i++)
            {
                lambda0 = lambda;

                double sinlambda = Math.Sin(lambda);
                double coslambda = Math.Cos(lambda);
                double sin2sigma = (cosU2 * sinlambda * cosU2 * sinlambda) + Math.Pow(cosU1sinU2 - sinU1cosU2 * coslambda, 2.0);
                double sinsigma = Math.Sqrt(sin2sigma);
                double cossigma = sinU1sinU2 + (cosU1cosU2 * coslambda);
                sigma = Math.Atan2(sinsigma, cossigma);
                double sinalpha = (sin2sigma == 0) ? 0.0 : cosU1cosU2 * sinlambda / sinsigma;
                double alpha = Math.Asin(sinalpha);
                double cosalpha = Math.Cos(alpha);
                double cos2alpha = cosalpha * cosalpha;
                double cos2sigmam = cos2alpha == 0.0 ? 0.0 : cossigma - 2 * sinU1sinU2 / cos2alpha;
                double u2 = cos2alpha * a2b2b2;

                double cos2sigmam2 = cos2sigmam * cos2sigmam;

                A = 1.0 + u2 / 16384 * (4096 + u2 * (-768 + u2 * (320 - 175 * u2)));
                B = u2 / 1024 * (256 + u2 * (-128 + u2 * (74 - 47 * u2)));
                deltasigma = B * sinsigma * (cos2sigmam + B / 4 * (cossigma * (-1 + 2 * cos2sigmam2) - B / 6 * cos2sigmam * (-3 + 4 * sin2sigma) * (-3 + 4 * cos2sigmam2)));
                double C = _flatten / 16 * cos2alpha * (4 + _flatten * (4 - 3 * cos2alpha));
                lambda = omega + (1 - C) * _flatten * sinalpha * (sigma + C * sinsigma * (cos2sigmam + C * cossigma * (-1 + 2 * cos2sigmam2)));

                double change = Math.Abs((lambda - lambda0) / lambda);

                if ((i > 1) && (change < 0.0000000000001))
                {
                    converged = true;
                    break;
                }
            }

            return _radiusB * A * (sigma - deltasigma);
        }

 

 

转载于:https://www.cnblogs.com/aqi-silence/p/3151143.html