【Python】GPXファイルからXYZファイルへの変換

## はじめに

&emsp;山形大学大学院修士1年のHagianです。この記事は、[Python Advent Calendar 2022](https://qiita.com/advent-calendar/2022/python) の15日目となります。

## やりたかったこと

1. GPXファイル（GPSやGNSSの情報をXML形式で記録したもの）から緯度・経度・高さを読み込む
2. UTM（ユニバーサル横メルカトル）座標に変換する
3. 変換した座標データをXYZファイル（X, Y, Z方向の基準点からの距離をスペース区切りで記録したもの）で書き出す

## GPXとXYZについて

### ①GPX

&emsp;GPXは英語で、GPS Exchange Formatといい、アプリケーション間でGPSやGNSSの情報をやり取りするために用いられるフォーマットです。**初耳だ、何やら難しそう…** と思われる方もいらっしゃるかもしれませんが、**一言でまとめるとXML形式** です。GPXでは軌跡やルート、ポイントを記録することができます。

https://www.topografix.com/gpx.asp

実際に使われている例として登山者向けアプリの「YAMAP」、「ヤマレコ」や、フィットネスアプリとして有名な「Strava」等でアクティビティの軌跡を出力する際に、この形式を選択することができます。

https://yamap.com/

https://www.yamareco.com/

https://www.strava.com/

以下にGPXファイルの中身の例を示します。富士登山を想定して作ったデータであり、架空のものですので実在しません（笑）

```xml:Subashiri.gpx
<?xml version='1.0' encoding='UTF-8'?>
<gpx creator="YAMAP - https://yamap.com" xsi:schemaLocation="http://www.topografix.com/GPX/1/1 http://www.topografix.com/GPX/1/1/gpx.xsd" version="1.1" xmlns="http://www.topografix.com/GPX/1/1" xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance">
<metadata>
    <bounds maxlat="35.369134" maxlon="138.778500" minlat="35.366435" minlon="138.739171"/>
</metadata>
<trk>
    <name>須走口5合目〜8合目</name>
    <number>1</number>
    <trkseg>
        <trkpt lat="35.366435" lon="138.778500"><ele>1966.2</ele><time>2022-09-10T03:48:17Z</time></trkpt>
        <trkpt lat="35.366872" lon="138.778276"><ele>1969.2</ele><time>2022-09-10T03:50:50Z</time></trkpt>
        <trkpt lat="35.367560" lon="138.777829"><ele>1975.7</ele><time>2022-09-10T03:55:16Z</time></trkpt>
        .
        .
        .
        <trkpt lat="35.369134" lon="138.739171"><ele>3362.7</ele><time>2022-09-10T07:45:39Z</time></trkpt>
    </trkseg>
</trk>
</gpx>
```

### ②XYZ

&emsp;XYZはアルファベット最後の3文字…ではなくて、XYZの座標を記録したファイルのことです。単位はメートルで、ある位置（基準点）からの距離を表します。XYZファイルはあまり一般的ではない形式ですが、中身はいたって簡単でCSVの区切り文字がカンマから半角スペースになっているだけです。主にGISソフトで用いられるほか、シミュレーションなどをする際の3次元プロットの元データとしても用いられます。
以下に、XYZファイルの一例（先述したGPXファイルをXYZに変換したもの）を示します。

```
298175.47862468986 3915945.848919105 1966.2
298156.21326954296 3915994.7826505005 1969.2
298117.31068591564 3916072.0151253897 1975.7
294608.64727569465 3916326.204215371 3362.7
```

左からX、Y、Zの順に並んでいます。

## 実装方法

### 開発環境

* Python 3.10.8
* Visual Studio Code 1.74.0

### 使用ライブラリ

* ```gpxpy```

&emsp;簡単に使うことのできるGPXファイルパーサーです。

https://pypi.org/project/gpxpy/

* ```pyproj```

&emsp;地理的情報を扱う際に、座標系・測地系の変換や方位角・距離の計算を簡単に行えるライブラリです。

https://pypi.org/project/pyproj/

### 緯度・経度・高さをUTM座標に変換

&emsp;まず先述した```gpxpy```を使って、GPXファイルから緯度・経度・高さのみを読み込みます。下記リンクを参考にさせていただきました。

https://aotoshiro.jpn.org/2021/2376

以下の例は読み込んだ緯度・経度・高さをCSVとして書き出すプログラムです。

```python:gpx2utm.py
import csv
import gpxpy

with (open("Subashiri.gpx",mode="r") as gpx_file, open("Subashiri.csv",mode="w") as csv_file):
    coordinate_list = []
    gpx = gpxpy.parse(gpx_file)
    for trk in gpx.tracks:
        for seg in trk.segments:
            for pt in seg.points:
                coordinate_list.append([pt.latitude, pt.longitude, pt.elevation])
```

&emsp;次に、読み込んだ緯度・経度・高さをUTM座標に変換します。UTM（ユニバーサル横メルカトル）については下記リンクを参照ください。

https://www.esrij.com/gis-guide/coordinate-and-spatial/coordinate-system-japan/#:~:text=UTM%20%E5%BA%A7%E6%A8%99%E7%B3%BB,%E3%81%A7%E6%8E%A1%E7%94%A8%E3%81%95%E3%82%8C%E3%81%A6%E3%81%84%E3%81%BE%E3%81%99%E3%80%82

&emsp;UTM座標系では経度によって「UTMゾーン」というものが決められています。計算式で求めることもできるのですが、ここではUTMゾーンがわかっている前提で実装します。下記リンクを参考にさせていただきました。

https://qiita.com/spicy_HotChicken/items/4cdf303493d73e24dc14

```python:utm2csv.py
import csv
from pyproj import Proj

utm_zone = input("UTMゾーンの値を入力（半角数字）：")
# WGS84楕円体に準拠する場合、ellps="WGS84"とすればよい。
utm_conv = Proj(proj="utm", zone=utm_zone, ellps="GRS80")
out_list = []

with (open("Subashiri.csv", mode="r") as input_file, open("Subashiri_utm.csv", mode="w") as out_file):
    reader = csv.reader(input_file)
    for row in reader:
        coordinate_list = [float(i) for i in row]
        lat = coordinate_list[0]
        lon = coordinate_list[1]
        ele = coordinate_list[2]
        # 経度は東西方向なのでX方向、緯度は南北方向なのでY方向のデータになる
        utmx, utmy = utm_conv(lon, lat)

        # 南半球ならオフセットする
        if lat < 0:
            utmy = utmy + 10000000
        else:
            pass

        out_list.append([utmx, utmy, ele])
    writer = csv.writer(out_file)
    writer.writerows(out_list)
```

&emsp;ここでは準拠する楕円体をGRS80楕円体にしましたが、WGS84楕円体にすることも可能です。

### XYZファイルへの書き出し

&emsp;XYZファイルは拡張子こそ特殊ですが、区切り文字がスペースになったCSVのようなものです。そこで標準ライブラリのCSVを用いて書き出しを行い、保存するファイルの拡張子を```.xyz```にする裏技（？）を使います。

## 実装例

&emsp;完成したコードを以下に示します。

```python:gpx2xyz.py
# GPXファイルから緯度経度標高を抽出してXYZファイルに書き出すプログラム
import csv
import gpxpy
from pyproj import Proj

def from_gpx_to_xyz():
    file_name = input("GPXファイルの名称を入力（拡張子除く）：")
    utm_zone = input("UTMゾーンの値を入力（半角数字）：")
    utm_conv = Proj(proj="utm", zone=utm_zone, ellps="GRS80")
    xyz_list = []

    with (open(f"{file_name}.gpx",mode="r") as gpx_file, open(f"{file_name}_from_gpx.xyz",mode="w") as xyz_file):
        gpx = gpxpy.parse(gpx_file)
        for trk in gpx.tracks:
            for seg in trk.segments:
                for pt in seg.points:
        # 経度は東西方向なのでX方向、緯度は南北方向なのでY方向のデータになる
                    utmx, utmy = utm_conv(pt.longitude, pt.latitude)
                    # 南半球ならオフセットする
                    if pt.latitude < 0:
                        utmy = utmy + 10000000
                    else:
                        pass
                    xyz_list.append([utmx, utmy, pt.elevation])
        writer = csv.writer(xyz_file, delimiter=" ")
        writer.writerows(xyz_list)

if __name__ == "__main__":
    from_gpx_to_xyz()
```

&emsp;ファイルの名称とUTMゾーンを指定するだけで、GPXファイルからXYZファイルに一発変換するツールを作ることができました。

## おわりに

&emsp;実装例のコード、およびテストに用いたGPXファイルとXYZファイルはGithubで公開しています。

https://github.com/PyRadiolarus/for-qiita/tree/master/GPX2XYZ