原來使用超圖軟件進行道路長度的統計，后來決定在項目中使用開源的工具，兩點之間距離量算的方法，網上很容易找到，基本上是把地球作為一個標準球體，使用固定半徑來進行計算，這樣的每一個計算結果與實際長度差距不是很大，一般是可以接受的，但是計算對象較多時統計的結果與超圖的計算結果對比就會發現有較大的差距，這是因為超圖的長度計算中地球是作為橢球體處理的。

　　下面提供一種橢球體表面兩點間的距離量算，“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 Distancepublic 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 Distancepublic 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

總結

以上是生活随笔為你收集整理的空间两点间距离计算的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。