目安時間：約
2分
地図を扱っていると、2点間の距離を計算することがよくある。
平面上での2点間の距離は、学生時代に数学で習った式で計算することができる。
ただ、それでは正確な距離とは言えない。
地球の丸みを考慮した距離を計算したいときにはどうすればよいか。
方法はたくさんあるが、最も簡単でそれなりの精度が計算できる、
ヒュベニ（Hubeny）の公式を用いることが多い。
詳細な計算式や概念は他のサイトを探せばたくさん見つかるため
ここでは割愛させていただく。
例えば以下のサイトはとても分かりやすい。
参考サイトではJavaの実装が書かれていたが、
私はC++を使っているため、C++でHubenyクラスを作成した。
|
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 |
#ifndef HUBENY_H #define HUBENY_H class Hubeny { public: const double L_RAD = 6378137.0; // WGS84 const double S_LAD = 6356752.314245; // WGS84 inline double toRadian(const double &value) const { return value * M_PI / 180.0; } inline void hubeny_distance(const double &x1, const double &y1, const double &x2, const double &y2, double &d, double &an) const { double e = pow(sqrt((pow(L_RAD, 2) - pow(S_LAD, 2)) / pow(L_RAD, 2)), 2); double dx = toRadian(x2) - toRadian(x1); double dy = toRadian(y2) - toRadian(y1); double u = (toRadian(y1) + toRadian(y2)) / 2.0; double W = sqrt(1.0 - e * pow(sin(u), 2)); double M = (L_RAD * (1.0 - e)) / pow(W, 3); double N = L_RAD / W; d = sqrt(pow(dy * M, 2) + pow(dx * N * cos(u), 2)); an = atan2(dy * M, dx * N * cos(u)) * 180 / M_PI + 360; if (an >= 360) an = an - 360; return; } }; #endif // HUBENY_H |
使い方は簡単である。
|
1 2 3 4 5 |
double d; double an; Hubeny h; h.hubeny_distance(139.74472, 35.65500, 140.09111, 36.10056, d, an); std::cout << "d:" << d << ", an:" << an << std::endl; |
計算すると以下の値となり、参考サイトと同じ値となった。
|
1 |
d:58502.5, an:57.678 |
今度時間があればヒュベニ（Hubeny）以外の計算も試してみたいところ。
