久しぶりに開発途中のネタです。
POI の最寄検索や、検索結果表示で現在値からの距離を表示させたくて、任意2点の緯度経度から距離を計算する方法を調べてみました。
地球が球体なら何となく頭をひねれば分かるかもしれませんが、楕円体の緯度経度から距離を求めるのは私には見当もつきません。
ということでググってみると、楕円体を考慮していて、世界中の緯度経度で通用しそうな計算式が2つ見つかりました。
- ヒュベニの公式 (参考HP: 日本は山だらけ～ 技術研究本部)
- Lambert-Andoyerの公式 (参考HP: 測地線航海算法(Geodesic Sailing))
何となく航海で使用する下の公式のほうが正確な気はしますが、上のヒュベニの公式の方が演算が簡単そうな気がします。
また、qgmapでは多少の誤差も許容できるので数％程度(?)の精度で十分です。そこで近似式も候補に入れてみます。
- 簡易近似式 (参考HP: 緯度・経度からの距離計算)
これらの計算式がどの程度の正確さなのか、ランダムに選んだ緯度経度で計算して誤差と計算速度を検証してみました。
とはいえ、検証しようにも正解がわからないので、とりあえず国土地理院のサイトで2点間の距離を計算できるようですので、そのサイトとの誤差を調べました。つまり、実際に正確かどうかの検証ではなくて、国土地理院のサイトの公式に近いかどうかの検証です (^^;
で、結果は以下のとおり。
アルゴリズムの違いによる国土地理院サイトとの誤差
縦軸は国土地理院サイトで計算した距離(a)と、公式で計算した値(b)の誤差の割合(|a-b|/a) で、横軸が2点間の距離を示しています。
さすが航海算法と呼ばれるだけありLambert-Andoyerの公式は距離によらず 10^-5 以下の誤差率になっています。
ヒュベニの公式は50km以下ではLambert-Andoyerの公式と同等ですが、距離が離れるほど誤差が大きくなっていくようです。
簡易近似式は1000km以下では誤差0.1%前後で、1000kmを越えるとヒュベニの公式と同様誤差が増えていきます。
また、それぞれの公式を Zaurus で計算させた時の時間も測定しました。ちなみに、スペシャルカーネルを使っているので、浮動小数演算が高速化される FastFPE パッチが適用されています。
測地線航海算法のHPにはLambert-Andoyer と式的には同じですが、計算量が少し少ない小野の公式というものも載っていたので、こちらも計測しました。
| 公式 |
1万回の計算時間 |
| 近似式 | 0.2838 sec |
| ヒュベニの公式 | 0.5043 sec |
| Lambert-Andoyerの公式 | 2.8202 sec |
| 小野の公式 | 2.6518 sec |
やはり近似式は演算量が少ないので高速ですね。1万回の計算が0.3秒なら充分実用的です。
ヒュベニの公式は近似式に比べ 1.8倍、Lambert-Andoyerの公式が10倍、小野の公式が9.5倍の時間がかかりました。
やはり精度が高い下の2つの公式は三角関数の演算などが多いため遅いですね。ただ、同じ結果が得られる小野の公式は少し早いようです。
結論
近似式が思いのほか精度が高かったので近似式を使おうと思います。
もともと精度も1%程度あれば十分かなと思っていたので、qgmap には充分ですね。
開発は案外こういう自己満足的な調査が楽しかったりするんですよねー。
ううん凄い。ソフトの開発ってのはこういうことをやっていくんですね。
「何を作っているんでしょうかっ」という番組は，無機質なものから有用なものに変身して行く過程が面白く，よく見ているのですが，今回の書き込みもまたそんなワクワク感があります。
chidaさん、お久しぶりです。
グラフとか書いて凄そうに見えますが、こんなことをしなくても開発はできてしまいます。ある意味冗長な(無駄な？)調査かもしれないですね。
ただ、私はものづくりのこういう過程が好きなんですよねー。
たとえ他の人には気づかれなくても、「こんなところにこだわった」という自己満足的な達成感というか (^^;
こういう知的好奇心も自己満足も満たせるからプログラムってやめられないんですよね (^^)
どうもいつもお世話になっております。
ぷちのいずさんのサイトがあるから、
未だにザウルスが現役でいられます。
今ではもう使ってる人はかなりいなくなったと思いますが、
私の使用目的ではザウルス以上の機器はありません。
qgmap、どんどん便利になっていきますね！
応援するくらいしかできませんけど、がんばってください！
応援コメントありがとうございます!
エンジニア冥利に尽きますね (^^)
こういうコメントってすごく嬉しいんですよ。
qgmap はまだまだ開発途上なので、まだまだ開発は続けるつもりですので
これからも応援よろしくお願いします！
こんにちは。私は以前、社内の地質調査データの検索ソフトを作ろうと思い、緯度経度から２点間の距離を計算する方法を調べて、あまりの難しさからあきらめてしました。構想としては、住所を入力したら、Google EarthのAPIを利用して緯度経度を取得し、その点から「半径1km以内」にある緯度経度を持つデータをデータベースから取得するというものでした。地球が丸いなんて無視して、そのポイントに近いデータを順番に表示させれば全く問題ないんですけど「半径1km以内」にこだわって挫折したのを思い出しました^^
hiderin さん、お久しぶりです。
返信遅くなってすいません。
最寄りの POI 検索は、まさに hiderin さんがやろうとしていたことですね。
私もはじめは地図上での距離で近い順に表示させればいいかと思っていたんですが、どうしても検索結果に現在地からの距離も出したいと思い、本格的に調べ始めました。
地理系はこの辺が難しいですよね。
返事を頂きまして、ありがとうございます。素人の私の場合、それは夢物語に近いのですが、bucchiさんなら困難を一つずつ乗り越えて行かれるのでしょうね。
話がずれますけど、私の住む京都市では、市内中心部の住所は、通り名が挿入された住所が一般的に使われており、住所から緯度経度を取得する際に、これが障害になります。例えば郵便番号604-0944から変換すると京都市中京区橘町となる住所があるのですが、一般的には京都市中京区”押小路通柳馬場東入”橘町というような住所になっています。
データベースに登録してある住所からGoogle Earth COM APIを利用して緯度経度を取得したとき、それらのデータは、ほとんど正しく取得できず、手作業での登録となったのでした^^
いえいえ、地図・地理関係は私も完全に素人ですよ。もう Google 先生に頼りっぱなしです (^^;
京都ではそういう習慣があるんですね。勉強になります。
qgmap にも GM_Lite みたいに住所検索をいれたいと思っていて、国土交通省が公開しているデータベースを元に作成しようと思っているんですが、そのデータを見てみると、確かに “京都市中京区橘町” となっていますね。通り名は入っていないようです。
こういうローカルなルールは京都以外にもありそうですね。
そういう情報がまとまっていれば対応もできそうですが、現状ではちょっと難しそうですね。