<ちょっと間が開いてしまいましたが、自作のAndroidアプリにログ機能を搭載したという話を書きました。
alasixosaka.hatenablog.com
ところが、そのログは、緯度、経度、高度、時刻を順番にカンマ（,）で区切ったテキストをだらだらと書き出しただけのものなので、辿ったルートを地図上に表示させようとするとGPXなどの形式に変換する必要があります。
今回は、その変換用のプログラムの話です。
変換はPC上で行い、プログラムはPythonで記述しました。
GPX形式とは
まず、GPX形式に変換するといっても、それがどんな形式わからないと変換のしようもありません。
といっても、ほぼ素人の自分がここであれこれ解説するよりも、詳しく説明されているサイトを見てもらったほうが早いので、詳しくは下記のリンクを見てください。
tancro.e-central.tv
www.kunimiyasoft.com
要するに、XML形式の一種で、位置情報に特化したファイル形式ということが言えます。ですので、一つの区切りが”＜”で始まり、”/＞”で終わる形式になっています。
例えば、あるポイントは、
<trkpt lat="36.41274926854364" lon="139.42386455894342"><ele>310</ele><time>2010-09-18T23:14:56Z</time></trkpt>
のように記述されています。lat=の後が緯度、lon=の後が経度、eleで囲まれている部分が高度、timeで囲まれている部分が時刻です。
ですので、カンマで区切られた緯度、経度、高度、時刻を順番にこの形式で出力してやれば基本的にはGPXファイルに変換できることになります。
それ以外には、ヘッダーとか細かい部分がありますが、それらは、上の2つのリンクで見てもわかるように、適当に別のGPXファイルを見てコピペしても問題ないようです。
ただし、問題が一つあります。
高度は変換が必要
アプリの開発のところでも書きましたが、Android端末から得られるGPSの高度というのは、準拠楕円体高というもので、実際の標高とは少し違うものということです。これを実標高であるジオイド高に変換する必要があります。
楕円体高、ジオイド高についての詳細は下記のサイトを見てください。
blog.geolonia.com
ja.wikipedia.org
私はあんまりよく理解していません。とにかく何か変換が必要というくらいの認識です。
そして、日本で標高を計算するには、国土地理院が提供しているデータを用いる必要があるとのことです。
qiita.com
どんな計算をしているのかを理解し、自分でプログラムを組むのは結構大変そうですが、幸いなことにPythonでプログラムを作って公開されている方がおられますので、変換の計算はこのプログラムを利用させていただくことにしました。
github.com
上記のサイトのREAD.MEにも書いてありますが、国土地理院のデータはあらかじめ各自でダウンロードしておく必要があります。
fgd.gsi.go.jp
最新のデータはverが2_2になっています。参考サイトのプログラムはver2_1を使用していますのでその部分も書き換える必要があります。最初エラーになってしばらく何のことかわかりませんでした。
プログラムの詳細
定義など
まず、定数、グローバル変数を定義します。なお、importは省略しています。
GSI_GEOID_FILE_NAME = 'gsigeo2011_ver2_2.asc' # ジオイドデータはglobalで使う global geoid global geoidparam global glamn global glomn global dgla global dglo global nla global nlo
ここはオリジナルのプログラムとほぼ同じですが、ジオイドのファイル名が上にも書いたようにダウンロードしたバージョンが2_2でしたのでそれに合わせてあります。
ジオイドデータ読み込み
ここはオリジナルと全く同じです。
# ジオイドデータを読込みます。事前作成したデータファイルモジュール（geoidData.py）が存在しない時は
# 国土地理院のジオイドデータファイルから作成します。（次回以降の処理が若干早くなる）
def getGeoidData():
global geoid
global glamn
global glomn
global dgla
global dglo
global nla
global nlo
if (os.path.exists(os.path.join(os.getcwd(), 'geoidData.py')) == False):
print('国土地理院のジオイドファイルからデータファイルを作成します(初回のみ)')
if (os.path.exists(os.path.join(os.getcwd(), GSI_GEOID_FILE_NAME)) == False):
print('エラー:データファイルがありません')
print('同じフォルダに国土地理院のジオイドファイル「' + GSI_GEOID_FILE_NAME + '」を置いてください')
return False
createGeoidData()
else:
print('データファイルを読込中')
module = importlib.import_module('geoidData')
geoid = module.setGeoid()
misc = module.setMiscData()
glamn = misc["glamn"]
glomn = misc["glomn"]
dgla = misc["dgla"]
dglo = misc["dglo"]
nla = misc["nla"]
nlo = misc["nlo"]
# ヘッダのdglaの有効桁数違いで地理院プログラムと微小な差異があったので、より計算値に近づける
dgla = math.floor(dgla * (nla - 1)) / (nla - 1)
dglo = math.floor(dglo * (nlo - 1)) / (nlo - 1)
ジオイドデータの書き出し
ここもオリジナルそのままです。初回だけこの関数が実行され、実行されると、geoidData.pyというファイルができます。
# 国土地理院のジオイドデータからPythonのデータを書きだして、次回から使えるようにします。
def createGeoidData():
global geoid
global glamn
global glomn
global dgla
global dglo
global nla
global nlo
f = open(os.path.join(os.getcwd(), GSI_GEOID_FILE_NAME), 'r')
# f = open('gsigeo2011_ver2_1.asc', 'r')
line = f.readline()
# 20.00000 120.00000 0.016667 0.025000 1801 1201 1 ver2.1
linestr = '' + line
linestr = linestr.strip()
header = linestr.split(" ")
glamn = float(header[0])
glomn = float(header[1])
dgla = float(header[2])
dglo = float(header[3])
nla = int(header[4])
nlo = int(header[5])
geoid = {}
la = 0
lo = 0
while line:
line = f.readline()
linestr = '' + line
linestr = linestr.strip()
g = linestr.split(" ")
glen = len(g)
for i in range(0, glen):
if (g[i].strip() == ""):
continue
if (g[i] != '999.0000'):
geoid[str(la) + "_" + str(lo)] = float(g[i])
lo += 1
if (lo == nlo):
lo = 0
la += 1
f.close()
fw = open('geoidData.py', 'w') # 書き込みモードで開く
fw.writelines('def setGeoid():\n')
# geoid = [[999] * 1201 for i in range(1801)]
fw.writelines('\tgeoid = {}\n')
for la in range(0, nla):
for lo in range(0, nlo):
if (str(la) + "_" + str(lo) in geoid.keys()):
fw.writelines(
'\tgeoid["' + str(la) + '_' + str(lo) + '"] = ' + str(geoid[str(la) + "_" + str(lo)]) + '\n')
fw.writelines('\treturn geoid\n')
fw.writelines('def setMiscData():\n')
fw.writelines('\tmisc = {}\n')
fw.writelines('\tmisc["glamn"] = ' + str(glamn) + '\n')
fw.writelines('\tmisc["glomn"] = ' + str(glomn) + '\n')
fw.writelines('\tmisc["dgla"] = ' + str(dgla) + '\n')
fw.writelines('\tmisc["dglo"] = ' + str(dglo) + '\n')
fw.writelines('\tmisc["nla"] = ' + str(nla) + '\n')
fw.writelines('\tmisc["nlo"] = ' + str(nlo) + '\n')
fw.writelines('\treturn misc\n')
fw.close()
ジオイド高の計算
ここが、計算の主体となる部分です。ここもオリジナルそのままです。
緯度、経度を引数にして関数を呼び出すと、ジオイド高を返すというルーチンです。
# 緯度経度からジオイド値を求める
def getGeoidValue(lat, lon):
global geoid
global glamn
global glomn
global dgla
global dglo
global nla
global nlo
# 囲う矩形を求める
j = int(math.floor((lon - glomn) / dglo))
i = int(math.floor((lat - glamn) / dgla))
if (i < 0 or i >= nla or j < 0 or j >= nlo):
# print('エラー：緯度経度が範囲外です')
return 999.00
if ((not (str(i) + "_" + str(j) in geoid.keys())) or (not (str(i) + "_" + str(j + 1) in geoid.keys())) or (
not (str(i + 1) + "_" + str(j) in geoid.keys())) or (not (str(i + 1) + "_" + str(j + 1) in geoid.keys()))):
return 999.00
wlon = glomn + j * dglo
elon = glomn + (j + 1) * dglo
slat = glamn + i * dgla
nlat = glamn + (i + 1) * dgla
t = (lat - slat) / (nlat - slat)
u = (lon - wlon) / (elon - wlon)
Z = (1 - t) * (1 - u) * geoid[str(i) + "_" + str(j)] + (1 - t) * u * geoid[str(i) + "_" + str(j + 1)] + t * (
1 - u) * geoid[str(i + 1) + "_" + str(j)] + t * u * geoid[str(i + 1) + "_" + str(j + 1)]
Z = Z * 100000
Z = math.floor(Z + 0.5)
Z = Z / 100000
return Z
高度を変換してGPXファイルに書き出す。
元のプログラムはCSVファイルを書き出すようになっていましたが、ここをGPXファイルを書き出すように変更しました。
関数名はconvertCSVとそのままですが。
まず、ログファイルを読み込みます。読み込んだ値はdataという配列に格納されます。配列は1次元配列で、緯度、経度、高度、時刻の順に並んでいます。
ですので、配列要素の４の倍数の時が緯度、４の倍数+1の時が経度、４の倍数+２の時が高度、４の倍数+3の時が時刻となっています。
GPSデータの緯度はlat、経度はlan、楕円体高はEllapsoidHeightという変数に格納されます。そこで緯度、経度からジオイド高を読みだして、楕円体高-ジオイド高より標高を計算してelavationという変数に格納します。
そうしたら、 with open(outfile, "w") as f: 以下で、GPXファイルとして書き出しています。最初の部分は上に書いたおまじないで、適当なほかのGPXファイルからとってきた部分をコピペして書き出しています。
そのあとは、GPXファイルの形式に従って、緯度、経度、高度、時刻を書き出すようになっています。
# CSVファイル中の高さ(楕円体高)カラムを標高値に置き換えて出力します
def convertCSV(infile, outfile):
global geoid
# header check
hasInvalid = False
with open(infile) as f:
lines = f.readlines()
for s in lines:
data = s.split(",")
length = len(data)//4
for row in range(length):
temp = data[row*4+2]
EllapsoidHeight = float(data[row*4+2].replace(',', ''))
lat = float(data[row*4].replace(',', ''))
lon = float(data[row*4+1].replace(',', ''))
geoidval = getGeoidValue(lat, lon)
elevation = EllapsoidHeight - geoidval # 標高＝楕円体高-ジオイド高
if (geoidval == 999.00):
# ジオイド値が正常に取得できない
elevation = 999999
hasInvalid = True
data[row*4+2] = "{0:.3f}".format(elevation) # 高さカラムを差し替え
#fout.writelines(delimChar.join(row) + '\n') # 出力
print(data[row*4+2])
if (hasInvalid):
print('ジオイドの取得に失敗したデータがあります。確認ください。該当するデータは標高を 999999 にしています。')
else:
print('処理を正常に終了しました')
with open(outfile, "w") as f:
f.write('<?xml version="1.0" encoding="utf-8" standalone="no"?>"\n')
f.write('<gpx xmlns="http://www.topografix.com/GPX/1/0" xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance" '
'xsi:schemaLocation="http://www.topografix.com/GPX/1/0 http://www.topografix.com/GPX/1/0/gpx.xsd" '
'creator="Bike Route Toaster http://bikeroutetoaster.com/" version="1.0">\n')
f.write('<trk>\n')
f.write('<Name>'+infile+'</Name>\n')
f.write('<trkseg>\n')
leng = len(data)//4
for i in range(leng):
Time = convertTime(data[i*4+3])
f.write('<trkpt lat="'+data[i*4]+'" lon="'+data[i*4+1]+'">\n')
f.write('<ele>'+data[i*4+2]+'</ele>\n')
f.write('<time>'+Time+'</time>\n')
f.write('</trkpt>\n')
f.write('</trkseg>\n')
f.write('</trk>\n')
f.write('</gpx>\n')
時刻の変換
時刻については、元のログファイルとGPXファイルの形式が異なるので、これも変換が必要です。
convertTimeという変数を使って変換しています。ログファイルの時刻は、年が４桁、月、日がそれぞれ２桁、そして、時、分、秒もそれぞれ２桁の合計14桁の整数として表現されています。
そこで、元の数値を10の10乗で割った４桁が年、そして残りの数字を10の8乗で割った値が月、というように順番に割り算を繰り返し、時刻を求めてそれらをGPX形式に合わせた文字列として返しています。
def convertTime(timeSS):
time = float(timeSS)
year = time//1e10
yearS = str(year)
time = time - year*1e10
month = time//1e8
monthS = str(month)
time = time - month*1e8
day = time//1e6
dayS = str(day)
time = time - day*1e6
hour = time//1e4
hourS = str(hour)
time = time - hour*1e4
minute = time//100
minuteS = str(minute)
sec = time - minute*100
secS = str(sec)
timeS = yearS+"-"+monthS+"-"+dayS+"T"+hourS+":"+minuteS+":"+secS+"Z"
return timeS
メインルーチン
最後がメインルーチンです。
ファイルダイアログを表示して、変換するログファイルを指定し、出力ファイルはファイル名をそのままにして、拡張子を".log"から".gpx"に変換したファイルにしています。
getGeoidDataは元のプログラムのままで、ジオイドデータを一度だけ変換するものです。convertCSVが変換の主体で、ジオイド高から実際の標高を求めてGPXファイルの出力します。
# main
def main():
#args = sys.argv
# 引数とファイルの存在チェック
root = tk.Tk()
root.withdraw()
typ = [("ログファイル",".log")]
dir = 'h:\gps\logfile'
filein = filedialog.askopenfilename(filetypes=typ, initialdir=dir)
fileout = filein.replace(".log","raw.gpx",1)
# ジオイド情報の取得
getGeoidData()
# CSVファイルの高さ情報からジオイド高を引いたものを複製する
convertCSV(filein, fileout)
# 実行
main()
実行結果
実際に自転車で近所を走った時のデータを変換してみました。
軌跡については下の図のように問題ありませんでした。
ところが標高については、下のグラフのように変な値になっています。
ところどころでフラットな部分があって、正しく標高が取得できていないようです。GPSの精度の問題でしょうか？ 変換前のグラフと変換後のグラフを比べると20-30mほど標高に差があり、変換自体はうまくできているようです。
yamapのアプリなどでは、あとで標高を見てもこんなに変なことになっていないので、理由はよくわかりません。まあ、最悪軌跡（緯度と経度）がちゃんととれていれば、標高は後から計算することができるので何とかなるので、一応これで良しとしようかと思います。
下のグラフはカシミール３Dで標高を直した図です。こんな感じでできるのでちょっと手間がかかりますが、そもそもがバックアップが目的なので、まあこれでいいのでは、と思っています。
ちなみにCOROSのデータはこんな感じです。COROSの方がデータが滑らかな感じがしますけど。ほぼほぼあっている感じですね。
ちなみに、軌跡のスタート地点はここです。北摂のラルプデュエズということで最近ちょっと話題のポイントです。ちゃんとグーグルマップにも載っています。道路の先に見えるのは最近完成した安威川ダムです。
本物のラルプデュエズには比べるまでもないですが、それでもここを一気に登るのは結構きつく、いい練習になります。
本当は走り始めがここというわけではなく、自宅からここまで自転車でやってきて、やっとアプリのことを思い出して、ここからスタートしたというのが実際のところです。この日は、忍頂寺まで行って帰ってきました。まだ、コロナから完全に回復せず、咳は出るは、体は思いはで、どうも今一です。
