【C++】緯度経度を使って2点間の距離を計算してみた

地球上の2点間の距離を計算する

2点間の距離の計算と聞くと、最初に「三平方の定理」を思い浮かべる人が多いと思います。
しかし、緯度と経度の2点間の距離を計算するときは三平方の定理を使っても正しい結果を得ることができません。
これは地球が平面ではなく、楕円形であるためです。

そこで、国土交通省の距離計算にも使われている「ヒュベニの公式」を使います。

ヒュベニの公式

ヒュベニの公式とは、以下のような計算式です。

2点間の距離D = \sqrt{(Dy･M)^2+(dx･N・cosP)^2}\\\
\\\
\\\
\\\
D_y：2点の緯度（ラジアン）の差\\\
\\\
D_x：2点の経度（ラジアン）の差\\\
\\\
P：2点の緯度の平均\\\
\\\
M=\frac{Rx(1-E^2)}{W^3}：子午線曲率半径\\\
N=\frac{Rx}{W}：卯酉（ぼうゆう）線曲線半径\\\
W=\sqrt{1-E^2･sin^2P}\\\
E=\sqrt{\frac{Rx^2-Ry^2}{Rx^2}}：離心率\\\
\\\
Rx：長半径（赤道半径）\\\
\\\
Ry：短半径（極半径）\\\

長半径、短半径というのは、地球を楕円体として計算する際必要となる値で、Wikipediaによると以下のように定められています。

測地系	発表年	長半径（m）	短半径（m）
Bessel	1841年	6,377,397.155	6,356,079
GRS80	1980年	6,378,137.000	6,356,752.314 140
WGS84	1984年	6,378,137.000	6,356,752.314 245

測地系というのは、緯度と経度によって地球上の位置を表すために決められたシステムのことで、各国で色々な測地系が存在するみたいです。
ただ、世界的には国際測地学協会などが定めたGRS80と、アメリカ国防省が策定したWGS84が使われています。

今回は、GoogleMapやGPSで使われているWGS84の値を使用しました。

C++で2点間の距離を計算する

#include <iostream>
#include <string>
#include <cmath>


// 緯度、経度をラジアンに変換する
inline double toRadian(double value){
    return value * M_PI / 180.0;
}

// 2点間の距離を算出
double calculateDistance(std::pair<double, double> &x, std::pair<double, double> &y){
    double Rx = toRadian(x.first - y.first);     // 2点の緯度の差
    double Ry = toRadian(x.second - y.second);   // 2点の経度の差
    double latAve = (toRadian(x.first + y.first)) / 2.0;   // 2点の緯度の平均

    // 離心率を求める
    double longRad = 6378137.000;       // 長半径（赤道半径）
    double shortRad = 6356752.314245;   // 短半径（極半径）
    double E = std::sqrt((std::pow(longRad, 2) - std::pow(shortRad, 2)) / std::pow(longRad, 2));   // 離心率

    double W = std::sqrt(1.0 - std::pow(E, 2) * std::pow(sin(latAve), 2));

    double N = longRad / W;   // 卯酉線曲率半径
    double M = (longRad * (1 - std::pow(E, 2))) / std::pow(W, 3);   // 子午線曲率半径

    // 2点間の距離を求める（単位 m）
    double distance = std::sqrt(std::pow(Rx * M, 2) + std::pow(Ry * N * cos(latAve), 2));
    return distance;
}

int main (){
    std::pair<double, double> tokyo(35.6812362, 139.7671248);
    std::pair<double, double> shinjuku(35.6896067, 139.7005713);
    double pointDistance = calculateDistance(tokyo, shinjuku);

    std::cout << pointDistance << std::endl;
}

出力結果：6095.62

解説

まずインクルードについて説明します。
C言語やC++では、初期状態では標準関数も含め関数が全く定義されていないため、定義を後付けしなければなりません。
そのため#includeで任意のファイルを書かれた場所に展開します。
今回は以下の3つのファイルを定義しました。

#include <iostream>
#include <string>
#include <cmath>

<iostream>は「I/O入出力ストリーム」です。
定義することによって、cinの標準入力やcoutの標準出力を行うことができます。
<string>はString型を扱えるようにするファイルで、数値型を文字列型に変更するstd::to_stringなどもここに定義されています。
<cmath>はルートの計算や定義された円周率を扱えるようにするファイルです。

さて、続いて本題の説明に入ります。
今回のソースコードでは以下のようなtoRadianメソッドを作成しました。

inline double toRadian(double value){
    return value * M_PI / 180.0;
}

ラジアンとは、**「円（扇形）の弧の長さ（L）÷ 円の半径（r）」**によって求められる値のことです。
単位は「ラジアン」で、以下の式によって求められます。

θ = θ × \frac{π}{180}

2点間の距離を算出するcalculateDistanceメソッドでは、ヒュベニの公式に則って計算しています。
計算結果が正しいかどうかはこちらのツールで確認しました。

double calculateDistance(std::pair<double, double> &x, std::pair<double, double> &y){

calculateDistanceメソッドでポインタを引数として渡しているのは、元の変数の実体のみを渡すことで余計なメモリを食わないからです。
変数を引数として渡してしまうと、渡された変数のコピーが渡されることになり、元の変数の2倍メモリを食うことになります。
メソッドに何を渡すかによって、メモリの消費量が変わるんですね。

感想

今回C++を初めて触りましたが、String型を使うのにもインクルードが必要と知って驚きました。
またポインタの使い方も分からず苦戦しましたが、ゼロから学ぶC++というサイトを読んで何となく理解しました。
せっかく少し勉強したので、今後も何か作ってみたいなと思います。

参考

ヒュベニの公式について、以下のサイトを参考にしました。
TrailNote ： 2地点間の距離の計算
ヒュベニの間略式を用いた幾何補正の応用
【swift】緯度経度から距離を求める方法 | 株式会社イーガオ
ヒュベニの公式（Hybeny's Distance Formula）による緯度・経度の2点間の距離の算出 - butter-tigerのブログ
ラジアンとは何か？角度をラジアンに変換する方法が理解できる練習問題付き｜高校生向け受験応援メディア「受験のミカタ」