anorithのブログ

主にGIS系の記事を書く。

Terria (5) GPXをCZMLに変換して読み込んでみる

概要

  • GPXファイルをCZMLに変換してTerriaMapに読み込む。
  • 変換はpython(特にgpxpy)を使用する

参考

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を読み込むと以下のような感じ。下に出てくるバーで時間を動かして表示できます。 (背景が寂しかったので地理院の空中写真を読み込んでいます。)

youtu.be