2014年8月9日土曜日

速度計算

　Vikingはまずい、とわかった段階で考えた。いっそ自分で計算してやるか。Vikingにパッチをあててやろうか、と——。結局、パッチはやめた。個人的なパッチなのでVikingがバージョンアップするたびに、対応することになる。面倒だ。
　で、Emacsでやることにした。
　GPXには位置(経度、緯度)と時刻(グリニッジ時)がある。そこから二点間の速度を計算するのを自前でやるのはそれほど、むずかしくないはずだ。問題は位置情報から距離をもとめるやり方だ。
　一瞬、ピタゴラスの定理で求められるんじゃね？ と思ったのだが、もちろんちがう。まちがっている。だいたい単位がちがうじゃないか。
　で、調べてみると、いくつか方法はあるようで。

　どうやらVikingは「２地点間の距離と方位角 - 高精度計算サイト」と同じやり方らしく、似た値になる。ほかに有名な「カシミール 3D」で使っているという「ヒュベニの公式」というものもあるらしい。さらにRのライブラリにあるだろうと、見当をつけてみたところ、2Dでの距離関数も見つけた。
　さいわい「二地点の緯度・経度からその距離を計算する（日本は山だらけ〜）」というサイトに「ヒュベニの公式」のR実装があった。それを参考にelispで組み直した(引数の順番は変えた)。

(defun deg2rad(x)
  (/ (* x pi) 180.0))

(defun dist(lon1 lat1 lon2 lat2)
  "緯度1 経度1 緯度2 経度2"
  (let* ((a 6378137.0)
         (b 6356752.314140)
         (dy (deg2rad (- lon1 lon2)))
         (dx (deg2rad (- lat1 lat2)))
         (my (deg2rad (/ (+ lon1 lon2) 2)))
         (e2 (/ (- (expt a 2) (expt b 2)) (expt a 2)))
         (Mnum (* a (- 1 e2)))
         (W (sqrt (- 1 (* e2 (expt (sin my) 2)))))
         (M (/ Mnum (expt W 3)))
         (N (/ a W)))
    (sqrt (+ (expt (* dy  M) 2) (expt (* dx N (cos my)) 2)))))

(dist 36.10056 140.09111 35.65500 139.74472)
58502.45893124113

　まあ、だいじょうぶみたいだ。

ラベル: , , ,