自転車で北海道一周してきたのでGPSログをLaTeXで表示する

この記事はTeX & LaTeX Advent Calendar 2018 18日目の記事です. LuaTeX関係ありません.遅くなってしまって申し訳ありません.

はじめに

私はサイクリングが趣味で,休みの日はロングライドに出かけることもあります. そして,タイトルにある通り学生最後の夏休みに17日かけて室蘭市から北海道を反時計回りに出発し函館市までツーリングをしてきました.

f:id:muscle_keisuke:20181218162317j:plainf:id:muscle_keisuke:20181218162322j:plain
f:id:muscle_keisuke:20181218162328j:plainf:id:muscle_keisuke:20181218162342j:plain

私は自転車でロングライドする時には必ずGPSログを取ります.北海道一周のときにもこのログは取り続けていました. 自分が保存していたGPSログはGPXというフォーマットです. GPXファイルには位置,速度,標高などの情報がXMLスキーマベースで保存されています. 全てテキストデータになっています.

そこで思いつきました.

LaTeXGPSログを読み出せそう

ということで,GPXファイルをLaTeX上で読み込んで,そのルートをTikZで描画することにしました.

GPXファイルの中身をLaTeXで読む

GPXファイルについて

GPXファイルは以下のようなXMLスキーマベースで記録されています*1

<?xml version="1.0" encoding="UTF-8"?>
<gpx creator="StravaGPX" version="1.1" xmlns="http://www.topografix.com/GPX/1/1" xmlns:xsi="http://www.w3.org/2001/XMLSchema-instance" xsi:schemaLocation="http://www.topografix.com/GPX/1/1 http://www.topografix.com/GPX/1/1/gpx.xsd">
 <metadata>
  <time>2017-09-04T00:57:35Z</time>
 </metadata>
 <trk>
  <name>Morning Ride</name>
  <trkseg>
   <trkpt lat="42.3847010" lon="141.0319980">
    <ele>65.4</ele>
    <time>2017-09-04T00:57:35Z</time>
   </trkpt>
   <trkpt lat="42.3846550" lon="141.0323940">
    <ele>64.9</ele>
    <time>2017-09-04T00:57:39Z</time>
   </trkpt>
   ...
  </trkseg>
 </trk>
</gpx>

<trk>というタグがトラックの始まりを表します.トラックには時間や緯度経度,標高などが記録されています. それぞれの数値はタグ内に記録されており,対応するタグは以下の通りです.

  • <trkpt>(緯度経度)
  • <time>(時間)
  • <ele>(標高)

今回はPDF上にルートを表示したいので緯度経度情報のみ取り出すことにします.

LaTeXでファイルの中身を読み出す

ファイルの読み込み

LaTeXにはファイルを読み出すコマンドがあります. 以下のソースはファイルを読み込み1文ずつ\filelineに格納していくプログラムです.

\newread\file
 \openin\file=filename
 \loop\unless\ifeof\file
  \read\file to\fileline
 \repeat
\closein\file

これを使ってGPXファイルから緯度経度情報の一覧を取り出します.

ある文字列を含む文のみ取り出す

LaTeXでファイルから1文ずつ取り出すことはできました. 次は必要な文だけ取り出すことを考えます. 緯度経度情報はGPXファイルの<trkpt>タグ内に記録されています. 次のようなフォーマットです.

<trkpt lat="緯度" lon="経度">...</trkpt>

ここから緯度,経度の数値のみ取り出します. まずは<trkpt ... を含む文章を取り出します. xstringパッケージは文字列処理に関するマクロが定義されています. その一つの\IfSubStrコマンドは文中に特定の文字列が入っているかどうか判定することができます. これを使って条件分岐を行います.

\IfSubStr{\filename}{<trkpt}{
 ...
}{}

第一引数に文,第二引数に含まれるか判定する文字列を渡します. 第一引数に今回は取り出した文が入った\filenameを指定し,<trkptを第二引数に指定します. これで<trkptから始まる文のみに分岐して処理を進めることができます.

次に取り出した文から緯度経度の数値のみ取り出します. 同じくxstringパッケージの\StrBetweenコマンドを使えば 特定の文字列に挟まれた文字列を取り出すことができます.

\StrBetween{}{取り出したい文字列の前の文字列}{取り出したい文字列の後の文字列}[取り出した後に格納する変数]

まずは緯度を取り出します.以下のフォーマットなので,

<trkpt lat="緯度" lon="経度">...</trkpt>
\StrBetween{\fileline}{lat="}{" lon}[\gpxlat]%

で取り出し,\gpxlatに格納します. 経度も同じように

\StrBetween{\fileline}{lon="}{">}[\gpxlon]%

で取り出し,\gpxlonに格納します. これで,

<trkpt lat="42.3847010" lon="141.0319980">

から 42.3847010141.0319980のみ取り出すことができました.

地図を表示する

地図表現の手法について

地図を表示しなければルートを描画しても位置がわかりにくいのでルートの下に地図を表示します. 最初のアイデアは地図をベクタ画像で表示し,その上に緯度経度から平面座標に変換した点をプロットして点間を線で結び, ルートを描画する方法でした. しかし,後述する複雑な数式によって,単純な平面座標への変換を断念したため, 今回は正方形タイルによる地図表現*2という手法を使い地図の表示をします. この手法はGoogleにより考え出された手法であり,現在,様々な地図アプリに採用されています. 平面座標への変換がしやすいという点で今回採用しました.

正方形タイルによる地図表現

正方形タイルによる地図表現は地図を縦横方向に敷き詰めた正方形タイルで表現する手法です. ズームレベルというものを定めズームレベルが高くなるにつれ, 敷き詰める正方形タイルを多くしてより詳しく地図を表現します.正方形タイルのピクセル数は256 x 256で統一され, ズームレベル0を初期状態とし,タイル数1で世界を表現します.

http://www.trail-note.net/images/zoom0.png

© OpenStreetMap contributors

ズームレベル1で縦横それぞれ2倍のタイル数になり2 x 2になります.

f:id:muscle_keisuke:20181219000859p:plainf:id:muscle_keisuke:20181219001114p:plain
f:id:muscle_keisuke:20181219001130p:plainf:id:muscle_keisuke:20181219001139p:plain
© OpenStreetMap contributors

つまり,ズームレベルnでタイル数は縦横それぞれ 2^n x  2^nとなります.

タイルにはそれぞれタイル座標というものがあり,例えばズームレベル1であれば, 左上のタイルのタイル座標は(0,0)で,右下のタイルのタイル座標は(1,1)となります. タイル座標とは別にピクセル座標というものもあります.これは全タイルを含めた総ピクセル数の内, どこの座標を指し示すかを表す座標です.ズームレベル1の場合256 x 256 のタイルが2 x 2敷き詰められているので, 総ピクセル数は512 x 512となります.よって,タイル座標(0,1)の中心点のピクセル座標は(384,128)となります.

地図を入手する

この正方形タイルに対応した地図を無償で公開している国土地理院から画像ファイルを持ってきます.

地理院地図|地理院タイル一覧

様々な地図をズームレベル毎にダウンロードすることができます. 今回は標準地図を使わせてもらいました.また,ズームレベル毎にタイル画像を用意すると大量になってしまうので, 今回はズームレベル11で固定しました. ズームレベル11の北海道本土のタイルをあらかじめすべて用意します. 今回用意したタイルのタイル座標は(1818,731)から(1854,764)です. 37 x 33 = 1221枚のタイルをLaTeXで貼り付けることになります. 貼り付けは\foreachを使って行います. \spreadtileというマクロを定義し,第一引数,第二引数に左上のタイル座標を渡し,第三引数,第四引数に縦横方向のタイル数を渡します. 今回は(1818,731)から37 x 33のタイルを敷き詰めるので \spreadtile{1818}{731}{37}{33}のように引数を渡します.

\newcommand{\spreadtile}[4]{
  \foreach \x in {0,...,#3}{
    \foreach \y in {0,...,#4}{
      \pgfmathparse{int(#1 + \x)}
      \let\numTileX\pgfmathresult
      \pgfmathparse{int(#2 + \y)}
      \let\numTileY\pgfmathresult
      \node[outer sep=0pt, inner sep=0pt] (image) at (\x*90.31,-\y*90.31){\includegraphics[]{maps/11:\numTileX:\numTileY.png}};
    }
  }
}

タイル画像を貼る座標指定についてですが,画像サイズがピクセルになっているため, LaTeX(TikZ)内の座標系と座標に変換する必要があります. tikzpicture環境は座標の刻み幅をxyそれぞれ1mmに設定しています.

\begin{tikzpicture}[x=1mm, y=1mm]
...
\end{tikzpicture}

なので256 x 256ピクセルの画像をmm単位に変換します. LaTeXのデフォルトのdpi(dot per inch)は72なので,1インチあたり72ピクセルのドットが入ります. 1ピクセルあたり1/72インチなので,256/72インチ(=3.555...インチ\approx90.31mm)になります. つまり,PDF内に原寸大で貼り付けられたタイル画像は90.31mm x 90.31mmになります. 次に座標系の変換です.といってもy軸方向の正負の向きが逆なだけです. よって,(x,y)番目の画像を貼り付ける座標は(90.31x, -90.31y)となります.

これで,地図タイルを隙間なく敷き詰めることができました.

緯度経度を平面座標に変換する

緯度経度 -> ピクセル座標 -> mm座標に変換する方法

PDFに描画するために緯度経度を平面座標に変換する必要があります. 単純に緯度経度を平面座標に変換する式は以下のリンクにありました. https://vldb.gsi.go.jp/sokuchi/surveycalc/surveycalc/algorithm/bl2xy/bl2xy.htm しかし,数式が複雑過ぎてLaTeXで計算するのが困難です. 数式を書くのは簡単なんですけどね

上記の理由により,地図表現を正方形タイルによる方法に変更しました. 正方形タイルによる地図表現の場合,以下の情報が分かれば,平面座標を求めることができます.

  • \lambda : 緯度
  • \varphi : 緯度
  • z : ズームレベル

緯度経度からピクセル座標に変換する数式は以下の通りになります.


x = 2^{z+7}\left(\frac{\lambda}{180} + 1\right)


y = \frac{2^{z+7}}{\pi}\left(-\mathrm{tanh}^{-1}\left(\sin\left(\frac{\pi}{180}\varphi\right)\right) + \mathrm{tanh}^{-1}\left(\sin\left(\frac{\pi}{180}L\right)\right)\right)


L = 85.05112878

xを求める式を式1,yを求める式を式2とします.

引用:http://www.trail-note.net/tech/coordinate/

上記の式によって求まる座標は前節で説明したピクセル座標ですが,全世界における絶対座標になります. 今回は北海道のみしか地図に表示していないため,これをタイル内における座標に変換します. まず,求めた絶対座標を(x,y)としたときのタイル座標を求めます. タイルは256 x 256ピクセルなので x,yをそれぞれ256で割ると,割った商がタイル座標になります. 更に割った余りが求めたタイル座標内におけるピクセル座標となります. これを前節で説明した方法でmmに変換すると,PDFに表示された地図上に適切に点をプロットすることができます.

LaTeXによる実装

変換方法が確立したので,次はそれをLaTeXで実装します. しかし,LaTeXは複雑な計算をするとコンパイルに時間がかかってしまうので, 定数など計算可能な部分は先に計算します. 式1はこれ以上計算のしようがないのでこのままLaTeXで計算します. #2が経度,#3がズームレベルです.

\FPeval{\abpixcelX}{2^(#3+7) * (#2/180 + 1)}

問題は式2です.まず,式内の逆双曲線関数LaTeXに求める関数がないので,これを


\mathrm{tanh}^{-1}z = \frac{1}{2}\ln\left(\frac{1 + z}{1 - z}\right)

で変形します. すると,式2は


y = \frac{2^{z+6}}{\pi}\ln\frac{\left(1 - \sin\left(\frac{\pi}{180}\varphi\right)\right)\left(1 + \sin\left(\frac{\pi}{180}L\right)\right)}{\left(1 + \sin\left(\frac{\pi}{180}\varphi\right)\right)\left(1 - \sin\left(\frac{\pi}{180}L\right)\right)}

のように変形できる. L = 85.05112878より,


y = \frac{2^{z+6}}{\pi}\ln\frac{1.99627207622\left(1 - \sin\left(\frac{\pi}{180}\varphi\right)\right)}{0.00372792378\left(1 + \sin\left(\frac{\pi}{180}\varphi\right)\right)}


= \frac{2^{z+6}}{\pi}\ln\frac{\left(1 - \sin\left(\frac{\pi}{180}\varphi\right)\right)}{\left(1 + \sin\left(\frac{\pi}{180}\varphi\right)\right)} + \frac{2^{z+6}}{\pi}\ln\frac{1.99627207622}{0.00372792378}


= 2^{z+6}\left(2 + \frac{1}{\pi}\ln\frac{\left(1 - \sin\left(\frac{\pi}{180}\varphi\right)\right)}{\left(1 + \sin\left(\frac{\pi}{180}\varphi\right)\right)} \right)


= 2^{z+6}\left(2 + 0.31830988618\ln\frac{\left(1 - \sin\left(\frac{\pi}{180}\varphi\right)\right)}{\left(1 + \sin\left(\frac{\pi}{180}\varphi\right)\right)} \right)

となります.これをLaTeXで計算します.

  \FPeval{\abpixcelY}{2^(#3+6)/pi * (ln(1.99627207622*(1 - sin(pi*#1/180))/(0.00372792378*(1 + sin(pi*#1/180)))))}

これで絶対ピクセル座標が求まる.次にタイル座標を求め,タイル内ピクセル座標を求めます.

\FPeval{\tileX}{trunc(\abpixcelX / 256, 0)}
\FPeval{\tileY}{trunc(\abpixcelY / 256, 0)}

余りを求める式がFPパッケージにないため, 割った値の小数点以下を切り捨てた後,商を掛け,元の値との差を計算して余りを求めます. また,ピクセル座標は左上が原点なのに対し,LaTeX内に貼り付けるタイル画像は画像中心が原点なので, 求めたタイル内ピクセル座標に対し,座標系を変換した上で(-128, 128)の平行移動をさせます.

\FPeval{\localX}{trunc(\abpixcelX-(256*trunc(\abpixcelX/256,0)),0) - 128}
\FPeval{\localY}{128 - trunc(\abpixcelY-(256*trunc(\abpixcelY/256,0)),0)}

最後に敷き詰められた(x,y)番目の画像は(90.31x, -90.31y)なので, タイル内ピクセル画像を平行移動させ,ピクセル -> mmに変換して描画します. dpiが72で1インチ=25.4mmより,1mmあたりのピクセル数は72/25.4=2.83464566929で求まります. よって,

\newcommand{\convertmm}[2]{
  \pgfmathparse{\pixcelXonPDF/2.83468054}
  \let\mmXonPDF\pgfmathresult
  \pgfmathparse{\pixcelYonPDF/2.83468054}
  \let\mmYonPDF\pgfmathresult
}
\pgfmathparse{int(\tileX - \value{initTileX})*256 + \localX}
\let\pixcelXonPDF\pgfmathresult
\pgfmathparse{int(\localY - int((\tileY - \value{initTileY})*256})
\let\pixcelYonPDF\pgfmathresult
\convertmm{\pixcelXonPDF}{\pixcelYonPDF}

ついにプロットする点の真の座標(\mmXonPDF,\mmYonPDF)を求めることができました.

ルートの描画

LaTeXの限界

早速求めた点をプロットしようと思い,テストデータとして, 100km程度のロングライドのGPXデータを用意してコンパイルを走らせました. しかし,

! TeX capacity exceeded, ...

というエラーが起こりました. 作業領域として確保しているメモリーを使い果たしてしまったようです. 100km程度のロングライドでメモリーを使い果たしてしまうので, 2000km程度の北海道一周ツーリングをプロットするのは不可能です. なので,<trkpt...から始まる文を見つける度にカウンタをインクリメントし,100で割り切れる時のみ, プロットをすることにしました.

\newcounter{numPoint}
...
\IfSubStr{\fileline}{<trkpt}{%
 \pgfmathparse{int(mod(\value{numPoint},100))};
 \ifnum\pgfmathresult=0
  ....
  \filldraw [red, ultra thin] (\mmXonPDF, \mmYonPDF) circle (1);
  \ifnum\value{numPoint}=0\else
   \draw [red, line width = 7pt] (\mmXonPDF, \mmYonPDF) -- (previous);
  \fi
  \coordinate (previous) at (\mmXonPDF, \mmYonPDF);
 \fi
 \stepcounter{numPoint}

実際に描画した

北海道1周のGPXをすべて読み込んで,点のプロットをしてみました.

f:id:muscle_keisuke:20181219021343p:plain

100個に1つの割合で点をプロットしていますが,それでもかなりの精度でプロットできています.

f:id:muscle_keisuke:20181219133440p:plainf:id:muscle_keisuke:20181219133436p:plain
f:id:muscle_keisuke:20181219133441p:plainf:id:muscle_keisuke:20181219133445p:plain

PDFの実際のサイズは27000 x 240000になりました.

最後に

LaTeXでこれだけ作れたので,普通にWebアプリケーションとしてGPXビューワ作れそうだなと思いました. ちなみに今回の北海道1周のデータをコンパイルするのに30分近くかかります.遅すぎる.

実は北海道1周の後に日本縦断もしたので, そのデータも描画してみたかったです.今度やってみます.