緯度経度による距離計算に、pythonのnumpyを使ってみた

gpxのログから、 とある緯度経度に最も近い点を求めるのに、 pythonのnumpyを使ってみた。

今迄は、手持ちのPOIしか扱うことがなかったので、 点数が少なく、listの計算で不満は無かった。 しかし、地理院からゲットできる公共施設の住所と座標を扱おうとすると、 5000点(全国で、データによっては5万)を越えるので、 とても遅くなる。 これを解消したかった。

やったことは、 ログの緯度経度を、radianなnp.arrayにして、 対象座標との距離を計算して、結果をarrayに収める。 あとは、argmin()で、最小値のindexを求める。

def calc_dist_array(lad_ar, lod_ar, ck_rlad, ck_rlod):
    """Calc distance."""
    #

    d_x_ar = lad_ar - ck_rlad
    d_y_ar = lod_ar - ck_rlod
    dx_ave = (ck_rlad + lad_ar) / 2.0

    w_2_ar = 1.0 - 0.0066943800230118 * np.sin(dx_ave)**2
    m_m_ar = 6335439.32708317 / np.sqrt(w_2_ar) / w_2_ar
    n_n_ar = 6378137.0 / np.sqrt(w_2_ar)
    t_2_ar = n_n_ar * np.cos(dx_ave) * d_y_ar

    return np.sqrt(m_m_ar**2 * d_x_ar**2 + t_2_ar**2)


np_dist = calc_dist_array(
    np.radians(np.array(gpx_log_latigude)),
    np.radians(np.array(gpx_log_longitude)),
    math.radian(target_rat), math.radian(target_lon))

nearest_index = np_dist.argmin()

これを使って、gpsのログが通過する住所リストを作成してみた。

大阪府高槻市田能ｽﾊﾉ下11
京都府亀岡市安町野々神8
京都府南丹市八木町八木東久保29-1
京都府南丹市園部町小桜町47
兵庫県篠山市福住344-1
兵庫県篠山市日置385-1
兵庫県篠山市北新町41
兵庫県篠山市網掛429
兵庫県篠山市宮田240
兵庫県丹波市柏原町柏原1
兵庫県丹波市氷上町成松甲賀1
兵庫県丹波市柏原町柏原1

リストで処理するよりかなり早くなったが、 5万点になるとさすがに使えそうにない。
もっと良い方法があるかもしれない。

例えば、最小を判定できればよいので、 1degree=94km で処理するとか、 フィルタをかけて、5万点から桁落ちさせるとか、 色々工夫の余地はありそう。
あと、番地以下が不要とか、 公共施設の近所を通過しないログだと役立たずとか。

う、不満だらけだ。

まあ目安程度なら使えるかも。
まあ、よしとしておこう。