山行時のGPSログをPythonで分析する(2) 国土地理院APIからの標高取得
GPSの標高データは経緯度データよりも誤差が生じやすいと聞きます。実際どの程度の誤差があるのかを調べてみました。
比較対象とするのは、国土地理院が提供している標高値で、国土地理院のAPI[1]国土地理院, サーバサイドで経緯度から標高を求めるプログラムから取得できる値です。
Pythonからの利用例[2]GIS奮闘記, Pythonで国土地理院の標高API を使ってみよう！, 2019を参考にして、以下のコードでGPSログに対応する標高を取得します。ここで、latとlonは、それぞれGPSログから取得した緯度と経度の配列(ndarray)です(前回記事参照)。
import request
import numpy as np
ele_map = []
for i in range(len(lat)):
url = "http://cyberjapandata2.gsi.go.jp/general/dem/scripts/getelevation.php" \
"?lat=%s&lon=%s&outtype=%s" %(lat[i], lon[i], "JSON")
resp = requests.get(url, timeout=10)
data = resp.json()
ele_map.append(data['elevation'])
ele_map = np.array(ele_map)当該GPSログは約5000点のデータポイントがありますので、全データの標高を取得するには少々時間がかかります。15分くらい放っておくと終わりました。
距離と標高をグラフ化し、両標高を比較します。ここで、dis_cumは距離の累積和、eleはgpxログから取得した標高、ele_mapは先程取得した国土地理院地の標高で、いずれもndarrayです。
plt.plot(dis_cum, ele, label='GPS data', color='k')
plt.plot(dis_cum, ele_map, label='map data', color='r')
plt.xlabel('distance (m)')
plt.ylabel('elevation (m)')
plt.legend(loc=0)
plt.show()
ぱっと見た感じ、ほとんど差が無いように見えます。両者の差の分布をヒストグラムにしてみます。
ele_diff = np.subtract(ele, ele_map)
med = np.median(ele_diff)
std = np.std(ele_diff)
plt.hist(ele_diff, bins=20, color='lightblue')
plt.xlabel('difference of elevation (m)')
plt.axvline(med, linestyle='-', color='b', linewidth=0.5)
plt.axvline(med-2*std, linestyle='--', color='b', linewidth=0.5)
plt.axvline(med+2*std, linestyle='--', color='b', linewidth=0.5)
plt.show()
print('中央値 =', round(med,1), '(m)')
print('標準偏差 =', round(std, 1), '(m)')
標準偏差 = 8.9 (m)
差分の中央値は5.4mでした。40m以上差がある点もいくつかあります。標準偏差σが8.9mですので、今回GPSで直接取得した標高と、GPSで取得した経緯度に基づく地形図データの標高とでは、概ね-12.4mから23.2mの差があるようです(中央値±2σ)。40m弱ですね。GPS端末に表示される標高は、これくらいの誤差を含んでいると知っておくと良さそうです。
せっかく標高を取得したので、保存しておきます。
df_gpx['ele_map'] = ele_map
df_gpx.to_csv('output.csv', index=False)References
