Terria (5) GPXをCZMLに変換して読み込んでみる
概要
- GPXファイルをCZMLに変換してTerriaMapに読み込む。
- 変換はpython(特にgpxpy)を使用する
参考
- ルートヒストリー
- Python3.7によるGPXファイルの読み込みと移動距離の算出
- GIS実習オープン教材(CZML入門)
- CZML仕様
GPS記録について
gpxファイルは自分で過去に記録したものを使用しました(さすがにこれだけで個人を特定はできないと思いますが・・・自分のGPSデータを使用するのは気を付けないといけないですね・・・)。探せばネットにもたくさん落ちていると思います。
GPRロガーとしては「ルートヒストリー」を使用しました。シンプルでよいアプリだと思います。
入力に使用したGPXファイル(抜粋)
ルートヒストリーでGPX形式で出力すると以下のようになります。
<?xml version="1.0" encoding="UTF-8"?> <gpx version="1.1" creator="RouteHistory - https://www.ateow.com" xmlns="http://www.topografix.com/GPX/1/1"> <trk> <trkseg> <trkpt lat="35.77923847359902" lon="137.25219315313424"> <time>2022-05-04T06:25:02Z</time> <ele>353.5262222290039</ele> </trkpt> <trkpt lat="35.77924467425494" lon="137.25207535075572"> <time>2022-05-04T06:25:02Z</time> <ele>353.4094467163086</ele> </trkpt> <trkpt lat="35.77933668128351" lon="137.25101184773823"> <time>2022-05-04T06:25:06Z</time> <ele>350.3656341670015</ele> </trkpt> <trkpt lat="35.77934171424567" lon="137.2509021258634"> <time>2022-05-04T06:25:07Z</time> <ele>349.1141156737026</ele> </trkpt> <trkpt lat="35.779352662004264" lon="137.25073872997493"> <time>2022-05-04T06:25:08Z</time> <ele>348.62688992877827</ele> </trkpt> </trkseg> </trk> </gpx>
pythonでGPXをCZMLに変換
jupyter notebookで作業してます。GPXの読み込みはQiitaの記事を参考にpygpxで行い、czmlは「GIS実習オープン教材」を参考に作成。GPXは高さが標高で格納されていますが、czmlでは楕円体高での表現とする必要があるので標高⇒楕円体高の変換が必要です(czmlの仕様にまだ詳しくないので断言はできないですが、おそらく楕円体高しかだめっぽい?)。今回はそこまで精度を求めていないのでざっくり40mを足して変換しています。
czmlの中では時間に依存するポイントと軌跡を表現するためのポリラインの二つを作成しています。
変換に使用したnotebook
import json import os import glob import gpxpy import gpxpy.gpx from pytz import timezone import pandas as pd
def make_df_from_gpx(gpx_path): # https://qiita.com/toran/items/6b4bcd103292164e8e28 # タイムゾーン dt_tz = timezone('Asia/Tokyo') # 日時文字列形式 dt_fmt = '%Y-%m-%dT%H:%M:%SZ' # 配列の初期化 DateTime = [] Lat = [] Lng = [] Alt = [] # GPXファイルの読み込み with open(gpx_path,'r') as gpx_file_r: # GPXファイルのパース gpx = gpxpy.parse(gpx_file_r) # GPXデータの読み込み for track in gpx.tracks: for segment in track.segments: # ポイントデータリストの読み込み points = segment.points # ポイントデータの長さ N = len(points) # ポイントデータの読み込み for i in range(N): # ポイントデータ point = points[i] # データ抽出 datetime = point.time.astimezone(dt_tz).strftime(dt_fmt) lat = point.latitude lng = point.longitude alt = point.elevation # データ代入 DateTime.append(datetime) Lat.append(lat) Lng.append(lng) Alt.append(alt) df = pd.DataFrame(columns=["datetime", "lat", "lon", "alt"]) df["datetime"] = DateTime df["lat"] = Lat df["lon"] = Lng df["alt"] = Alt return df
def make_czml(timestamp, lat, lon, alt, diff_geoid=0): czml = [{ "id" : "document", "name" : "name", "version" : "1.0" }] elem = { "id": "1", "name": "gpx_point_with_timestamp", "description": "GPXの時刻付きポイントです", f"availability":f"{timestamp[0]}/{timestamp[-1]}", "billboard" : { "image" : "http://maps.google.com/mapfiles/kml/pushpin/ylw-pushpin.png", "scale" : 0.5 }, "position":{ "epoch" : f"{time_stamp[0]}", "cartographicDegrees":[] } } coord = elem["position"]["cartographicDegrees"] for i in range(len(timestamp)): coord += [timestamp[i], float(lon[i]), float(lat[i]), float(alt[i])+diff_geoid] czml.append(elem) elem_2 = { "id" : "2", "name" : "gpx_line", "description" : "GPX全体のポリラインです", "polyline" : { "positions" : { "cartographicDegrees" : [] }, "material" : { "solidColor" : { "color" : { "rgba" : [0, 0, 255, 100] } } }, "width" : 2.5 } } coord_2 = elem_2["polyline"]["positions"]["cartographicDegrees"] for i in range(len(timestamp)): coord_2 += [float(lon[i]), float(lat[i]), float(alt[i])+diff_geoid] czml.append(elem_2) return czml def save_czml(czml, out_dir, out_name): out_path = os.path.join(out_dir, f"{out_name}.czml") with open(out_path, 'w') as f: json.dump(czml, f, indent=4)
gpx_path = "./Log20220504-152502.gpx" out_dir = "./" df = make_df_from_gpx(gpx_path) out_name = os.path.basename(gpx_path).split(".")[0] lat = df["lat"].values lon = df["lon"].values alt = df["alt"].values time_stamp = df["datetime"].values # ざっくり40m程度足して標高⇒楕円体高の変換をする czml = make_czml(time_stamp, lat, lon, alt, diff_geoid=40) save_czml(czml, out_dir, out_name)
Terriaでの表示
Terriamapにドラッグ&ドロップでczmlを読み込むと以下のような感じ。下に出てくるバーで時間を動かして表示できます。 (背景が寂しかったので地理院の空中写真を読み込んでいます。)
- 地理院タイル一覧
- https://maps.gsi.go.jp/development/ichiran.html
- 「簡易空中写真(2004年~)」を読み込んでいます。たまたま画像がある箇所でよかった・・・。