Vincentyさん
かなりの高精度で距離を計測することができます。
用例
地点１と地点２の距離を求めたいなら、それぞれの緯度経度をdoubleで渡してあげます。戻り値もdoubleです。単位はメートル。
double distance = vincenty_distance(地点１の緯度, 地点１の経度, 地点２の緯度, 地点２の経度)
|
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 |
/// <summary> /// Vincentyの公式による２点間の距離計測 /// </summary> /// <param name="d_lat0"></param> /// <param name="d_long0"></param> /// <param name="d_lat1"></param> /// <param name="d_long1"></param> /// <returns></returns> public double vincenty_distance(double d_lat0, double d_long0, double d_lat1, double d_long1) { double earth_minor_radius = 6356752.0; double earth_major_radius = 6378137.0; double earth_flattening = 0.0033528106811823; //double acceptable_error = 1E-12;//0.006mm double acceptable_error = 1E-10;//0.6mm d_lat0 = deg2rad(d_lat0); d_long0 = deg2rad(d_long0); d_lat1 = deg2rad(d_lat1); d_long1 = deg2rad(d_long1); // 測地緯度から化成緯度を求める double rlat1 = Math.Atan(earth_minor_radius / earth_major_radius * Math.Tan(d_lat0)); double rlat2 = Math.Atan(earth_minor_radius / earth_major_radius * Math.Tan(d_lat1)); double sin_rlat1 = Math.Sin(rlat1); double cos_rlat1 = Math.Cos(rlat1); double sin_rlat2 = Math.Sin(rlat2); double cos_rlat2 = Math.Cos(rlat2); // 経度差を求める double dlng = d_long0 - d_long1; double lambda = d_long0 - d_long1; double prev_lambda = lambda + 1; double sin_sigma = 0; double cos_sigma = 0; double cos2_alpha = 0; double cos_2sigma = 0; double tmp = 0; double sigma = 0; // ラムダ値を計算する。 (値が収束するまで繰り返す) for (int i = 0; i < 100 && Math.Abs(lambda - prev_lambda) > acceptable_error; i++) { double sin_lambda = Math.Sin(lambda); double cos_lambda = Math.Cos(lambda); sin_sigma = Math.Sqrt(Math.Pow(cos_rlat2 * sin_lambda, 2) + Math.Pow(cos_rlat1 * sin_rlat2 - sin_rlat1 * cos_rlat2 * cos_lambda, 2)); cos_sigma = sin_rlat1 * sin_rlat2 + cos_rlat1 * cos_rlat2 * cos_lambda; /*double*/ sigma = Math.Atan2(sin_sigma, cos_sigma); double sin_alpha = cos_rlat1 * cos_rlat2 * sin_lambda / (sin_sigma == (float)0 ? 1 : sin_sigma); cos2_alpha = 1 - sin_alpha * sin_alpha; cos_2sigma = cos_sigma - 2 * sin_rlat1 * sin_rlat2 / (cos2_alpha == (float)0 ? 1 : cos2_alpha); double c = earth_flattening / 16.0 * cos2_alpha * (4 + earth_flattening * (4 - 3 * cos2_alpha)); tmp = sigma + c * sin_sigma * (cos_2sigma + c * cos_sigma * (-1 + 2 * Math.Pow(cos_2sigma, 2))); prev_lambda = lambda; lambda = dlng + (1 - c) * earth_flattening * sin_alpha * tmp; } // 必要な値を事前に計算 double u2 = cos2_alpha * (Math.Pow(earth_major_radius, 2) - Math.Pow(earth_minor_radius, 2)) / Math.Pow(earth_minor_radius, 2); double a = 1 + u2 / 16384 * (4096 + u2 * (-768 + u2 * (320 - 175 * u2))); double b = u2 / 1024 * (256 + u2 * (-128 + u2 * (74 - 47 * u2))); // デルタ値を求める tmp = 1.0 / 6.0 * b * cos_2sigma * (-3 + 4 * Math.Pow(sin_sigma, 2)) * (-3 + 4 * Math.Pow(cos_2sigma, 2)); double delta = b * sin_sigma * (cos_2sigma + 0.25 * b * (cos_sigma * (-1 + 2 * Math.Pow(cos_2sigma, 2)) - tmp)); double d_distance = earth_minor_radius * a * (sigma - delta); return d_distance; } |
One thought on “C#で距離を計測する：Vincentyの公式”
Comments are closed.