2地点の緯度・経度から距離を求める

　ラズパイピコでGPSを使えるようになったので、2地点の緯度・経度から距離を求めます。犬の散歩の距離をGPSで測るガジェットを作りたい。

地球を真球と仮定して計算

　緯度・経度から距離を求める方法はいくつかあるようです。
　まず最初に考えつくのは地球を真球と仮定して、地球表面の2点と地球の中心点が作る角度を使う方法です。緯度をφ、経度をθ、地球の半径をrとすると、地球表面の点は(X, Y, Z)=(r cosφ cosθ, r cosφ sinθ, r sinφ)と極座標で表せます。2地点A、Bの内積はXa Xb + Ya Yb + Za Zb = r^2 cosξとなりますので、ξ = acos((Xa Xb + Ya Yb + Za Zb)/r^2)と計算でき、2地点の距離はD = r ξ となります。
　なお地球を真球と仮定する方法は、google mapでも採用しているそうです。

　私の家と隣家の緯度・経度から、エクセルを使って距離を計算したところ、16.7mと出ました。
　ところがMicroPythonでプログラムを組んでみると、この位の近距離だとξが小さすぎて計算が上手くできず、距離は0mとなってしまいました。
　どうやらMicroPythonのfloat型は、有効数字6桁以上は誤差が生じるらしい。ちなみに本家Pythonでも17桁以上で誤差が生じるらしい。これは10進数の計算をCPU内で2進数で計算しているからで、原理的にある程度以上の桁数では近似値となるようです。17桁なら十分だけど、6桁だと数値計算にはちょっと厳しい・・・。
　本家PythonにはDecimalというモジュールがあり、数値計算の誤差を減らせられるようですが、MicroPythonにはDecimalモジュールをimportできませんでした（探せばどこかにあるかもしれないけれど）。

ヒュベニ（Hubeny）の式で計算

　次にHubenyの式を使ってみました。これはカシミール3Dで採用されており、地球の極部以外の近距離では精度が高いようです。精度についてはこちらのサイトが詳しい。ヒュベニの式では2地点の緯度差、経度差、平均緯度、地球の赤道半径、極半径、離心率、子午線曲率半径、卯酉線曲率半径から、2点間の距離を計算します。プログラムは以下の通りです。

import math

#ヒュベニの式で2点間の距離を求める
def hubenyDist(lat1,lon1,lat2,lon2):
    #赤道半径（WGS84）
    Rx=6378137
    #極半径（WGS84）
    Ry=6356752.314
    #緯度の差（ラジアン）
    Dy=math.fabs(lat1-lat2)/180*math.pi
    #経度の差（ラジアン）
    Dx=math.fabs(lon1-lon2)/180*math.pi
    #平均緯度（ラジアン）
    P=(lat1+lat2)/2/180*math.pi
    #離心率
    E=math.sqrt((Rx*Rx-Ry*Ry)/Rx/Rx)
    W=math.sqrt(1-E*E*math.sin(P)*math.sin(P))
    #子午線曲率半径
    M=Rx*(1-E*E)/W/W/W
    #卯酉線曲率半径
    N=Rx/W
    #2点間の距離
    return math.sqrt(Dy*Dy*M*M+Dx*Dx*N*N*math.cos(P)*math.cos(P))

　計算してみると、隣家との距離は18.4mと出ました。同じ計算をエクセルでやってみると、16.8mでした。うーん、計算方法の選定以前に、MicroPythonの計算誤差をもう少しなんとかしないといけない。。。

　追記：MicroPythonで使えるDecimalモジュールを見つけ、解決しました。