鉄道路線の線形を可視化する
思いついたときに地図ネタを置いておく。
概要
よく「線形が良い・悪い」という言葉を聞くが、小さな縮尺の地図でなんとなく眺めていても分からない。そこで、急カーブの部分だけ色をわけて描画できないかと思い、トライしてみた。
データの取得・変換
まずはデータが無いと始まらない。
鉄道路線の座標データは、国土地理院の国土数値情報ダウンロードサービス（要登録：無料）から落とせる。shapeファイルとXMLファイルが選べるが、今回はQGISで変換するので扱いやすいshapeファイルを選択した。
全国の路線データを一気に描画しようとするとPCに負担がかかりすぎるので、一部の領域だけデータを抜き出す。
ダウンロードしたデータを、QGISで開いてみる。
路線データが表示できた。この中で、可視化したい場所を選択して、レイヤパネルの「N02-16_RailroadSection」部分を右クリックして名前をつけて保存。選択した地物だけを、GEOJSON形式で保存しておく。
GEOJSON形式の中身はこんな雰囲気。駅と駅間がそれぞれ1本ないし数本の曲線で表されている。
{ "type": "FeatureCollection", "name": "test", "crs": { "type": "name", "properties": { "name": "urn:ogc:def:crs:EPSG::4019" } }, "features": [ { "type": "Feature", "properties": { "N02_001": "11", "N02_002": "2", "N02_003": "飯田線", "N02_004": "東海旅客鉄道"}, "geometry": { "type": "LineString", "coordinates": [ [ 137.82087, 35.51948 ], [ 137.82, 35.51893 ], [ 137.81924, 35.51846 ], [ 137.81848, 35.518 ], [ 137.81822, 35.51784 ], [ 137.81774, 35.51756 ], ... { "type": "Feature", "properties": { "N02_001": "11", "N02_002": "2", "N02_003": "飯田線", "N02_004": "東海旅客鉄道"}, "geometry": { "type": "LineString", "coordinates": [ [ 137.90375, 35.57788 ], [ 137.90433, 35.578 ], [ 137.90492, 35.57813 ], [ 137.90549, 35.57827 ], [ 137.90593, 35.57836 ], [ 137.90635, 35.5784 ], ... { "type": "Feature", "properties": { "N02_001": "11", "N02_002": "2", "N02_003": "飯田線", "N02_004": "東海旅客鉄道"}, "geometry": { "type": "LineString", "coordinates": [ [ 137.82166, 35.51996 ], [ 137.82087, 35.51948 ] ] } }, { "type": "Feature", "properties": { "N02_001": "11", "N02_002": "2", "N02_003": "飯田線", "N02_004": "東海旅客鉄道"}, "geometry": { "type": "LineString", "coordinates": [ [ 137.84117, 35.49633 ], [ 137.84192, 35.49717 ] ] } }, { "type": "Feature", "properties": { "N02_001": "11", "N02_002": "2", "N02_003": "飯田線", "N02_004": "東海旅客鉄道"}, "geometry": { "type": "LineString", "coordinates": [ [ 137.91295, 35.5953 ], [ 137.91348, 35.5961 ] ] } }, { "type": "Feature", "properties": { "N02_001": "11", "N02_002": "2", "N02_003": "飯田線", "N02_004": "東海旅客鉄道"}, "geometry": { "type": "LineString", "coordinates": [ [ 137.83951, 35.48775 ], [ 137.83988, 35.48874 ], [ 137.83998, 35.48905 ], [ 137.84003, 35.48935 ], [ 137.83992, 35.49093 ], [ 137.83972, 35.49333 ],... ...
度単位→メートル単位に変換
座標が度単位なので、メートル単位に直す。投影法を気にし始めればきりがないが、狭い範囲なら地球を球と近似して中心点からの経緯度差をそのまま距離に変換する以下の式で十分、だと思う。
# 中心点は最初の曲線の開始点 cp = numpy.array(data[u'features'][0][u'geometry'][u'coordinates'][0]) R_e = 6378137.0 #地球半径 R_ns = R_e * math.pi / 180.0 # 南北方向角度1degあたり長さ R_ew = R_e * math.pi / 180.0 * math.cos(math.radians(centerPoint[1])) # 東西方向1degあたり長さ # i番目の曲線について：イテレータ省略 points = numpy.array(data[u'features'][i][u'geometry'][u'coordinates']) # 経緯度をm単位に points_m = numpy.array(points) * 0.0 points_m[:,0] = R_ew (points[:,0] - cp[0]) points_m[:,1] = R_ns (points[:,1] - cp[1])
曲率半径の解析
3点を通る円は一意に定まる。曲線上のある点について、その点と前後の点合わせて3点を通る円の半径を求めて、その点における線路の曲率半径を求めることが出来る。式としては以下の方程式を解いて中心点の座標を求め、中心から曲線上の点までの距離を求める感じ。（導出略）
pythonで実装してみる。1次方程式はnumpyを使って解かせる。3点が同一直線に近い場合は別に処理を書く必要がある。なお、隣接する点だと誤差が出やすいため、3つ隣の点で算出してみている。
#k番目の要素まわりの曲率半径 A = numpy.zeros((2,2)) v = numpy.zeros(2) # k-3番目の点がxp, k番目はx, k+3番目はxn xp = points_m[max(k-3,1),0] yp = points_m[max(k-3,1),1] x = points_m[k,0] y = points_m[k,1] xn = points_m[min(k+3,len(points_m)-1),0] yn = points_m[min(k+3,len(points_m)-1),1] A[0,0] = 2.0 * (x - xp) A[0,1] = 2.0 * (y - yp) v[0] = (x * x - xp * xp) + (y * y - yp * yp) A[1,0] = 2.0 * (xn - x) A[1,1] = 2.0 * (yn - y) v[1] = (xn * xn - x * x) + (yn * yn - y * y) # 中心座標 # 行列式がゼロ: 3点が同一直線上 便宜上曲率半径9999mとする if abs(numpy.linalg.det(A)) < 1e-10: p = [0.0, 0.0] r[k] = 9999.0 else: p = numpy.linalg.solve(A, v) r[k] = math.sqrt( (x - p[0])*(x - p[0]) + (y - p[1])*(y - p[1]) )
描画
えられた曲率半径を使って、matplotlibで描画してみた。
for k in range(0,len(points_m)-1): plt.plot([points_m[k,0],points_m[k+1,0]], [points_m[k,1],points_m[k+1,1]], color=cm.hot(0.9 * max(0.0, (r - 100.0)) / 4000.0 ))**0.2)) #100.0:最小半径 4000.0:最大半径 0.2:パラメータ(このへん調整)
うまく描画するためのパラメータ調整が難しいけど、半径の0.2乗くらいの値を用いると急カーブから緩いカーブまでうまくグラデーションにできそうな感触。パラメータが大きすぎると急カーブばかり濃い色になる。
Graphillionで最長片道切符のルート探索(3) 現時点での本州内のルートとか
久しぶりにGraphillionをいじっていたら、新しいバージョンが公開されていたので、以前挫折したWindowsPCへの導入を試みたところ、あっさりと成功してしまった。ということで、ルート探索でいくつか遊んでみた。
この記事の執筆時点(2015/01/26)における最長片道切符の本州内ルート(中小国→小倉の最長ルート)は、以下のようになる。(Yahoo!地図をもとに作成)
赤いほうが営業キロベース、青い方は運賃計算キロベースで計算した結果である。とくに東日本で大きくルートが異なっている。
運賃計算の特例や、災害で長期運休となっている区間の扱いについては、
- 岩国・櫛ケ浜間で山陽本線(柳井経由)を使用しない
- 新下関・小倉間で新幹線と在来線の同時使用不可
- 釜石・宮古間およびいわき・岩沼間は代行バスが無いため使用しない
- 気仙沼線BRT仮復旧区間は、営業キロは鉄道線のものを使用、運賃計算キロは0kmとする
とした。
常磐線(竜田-原ノ町)の代行バスが運行開始すれば、このルートは大きく変化して以下のようになる。
北陸新幹線の開業後も、在来線が新幹線に置き換わる形でのルート変更しか無いが、仙石線と東北本線の連絡線(仙石東北ライン)が開業すると、東北全体でルートが変化するようだ。ただ、こちらは運賃の扱いがまだ発表されていないようなので掲載はしない。
さらにおまけとして、常磐線代行バス運行開始後の、出発地に戻ってくるルートのうち最長のものを載せておく。これの列挙は、片道ルートでpaths(start, goal)としていた部分をcycles()にすれば良いだけなので簡単だ。
普通の最長片道切符では、自宅から稚内まで行き、肥前山口からまた帰るという手間がかかるが、このルートならどこから始めても出発駅に戻って来られるので、ある意味お得かもしれない。
続・GraphillionをWindowsにインストールしようとしている→できた
追記(2015/1/25)
最新バージョンのGraphillionを使い、https://github.com/takemaru/graphillion/wikiに従えば問題なく無事インストールに成功し、チュートリアルも動作した。めでたしめでたし。
以下、過去の出来事
Linux機が挙動不審なためメイン機(Win7 64bit)に何とかしてGraphillionの環境を整えようとすることに。
開発者さんの日本語解説https://github.com/takemaru/graphillion/wikiを読むと、WindowsではAnacondaを利用するのが推奨されていた。前回は何とか使わないでもいけないかと試行錯誤してみたが、思いきってPythonをアンインストールしてしまい、Anacondaを入れなおすことを決意した。
あとは解説に書いてあるとおり、パスを通してIPythonからコマンドを打ってビルドしてみた。
これで何とかなるか・・・と思ったら、ビルドは警告を出しつつも成功したのだが、テストしてみようとすると途中で「python.exeは動作を停止しました」の文字が・・・
その際のコンソールへの表示は以下のとおり。
In [12]: run setup.py test running test running egg_info writing Graphillion.egg-info\PKG-INFO writing top-level names to Graphillion.egg-info\top_level.txt writing dependency_links to Graphillion.egg-info\dependency_links.txt reading manifest file 'Graphillion.egg-info\SOURCES.txt' reading manifest template 'MANIFEST.in' writing manifest file 'Graphillion.egg-info\SOURCES.txt' running build_ext copying build\lib.win-amd64-2.7\_graphillion.pyd -> test_binary_operators (graphillion.test.graphset.TestGraphSet) ...
ここまで表示されて落ちてしまった。
テストには失敗するものの、これをスキップしてrun setup.py installでインストールするまではエラーも出ず成功した。
さあ、とにかくテストを兼ねてチュートリアルを実行してみようとした。
In [1]: from graphillion import GraphSet as gs In [2]: import graphillion.tutorial as tl In [3]: univ = tl.grid(8,8) In [4]: gs.set_universe(univ)
ここまでは問題なく行ったが、次の
In [5]: tl.draw(univ)
で、以下のエラーメッセージが。
---------------------------------------------------------------------------
KeyError Traceback (most recent call last)
<ipython-input-4-574648453ea0> in <module>()
----> 1 tl.draw(univ)
C:\Anaconda\lib\site-packages\graphillion-0.96-py2.7-win-amd64.egg\graphillion\t
utorial.pyc in draw(g, universe)
59 if not isinstance(universe, nx.Graph):
60 universe = nx.Graph(list(universe))
---> 61 n = sorted(universe[1].keys())[1] - 1
62 m = universe.number_of_nodes() / n
63 g.add_nodes_from(universe.nodes())
C:\Anaconda\lib\site-packages\networkx\classes\graph.pyc in __getitem__(self, n)
317 {1: {}}
318 """
--> 319 return self.adj[n]
320
321
KeyError: 1さらに、
In [6]: paths = gs.paths(1,81)
を実行しようとしたら、また「python.exeは動作を停止しました」・・・。
ああもう。
Graphillionで最長片道切符のルート探索(2) 本州内のルート探索
前回、Graphillionを用いて北海道内のルートを探索することに成功したが、次は本州内のルートを探索してみる。
とりあえず計算
本州内の最長ルートということで、前回と似たようにまずはデータを準備して、中小国→小倉のルートを数え上げた上で最長のものを表示すればよい。
しかしまず最初の問題として、本州の路線図は北海道のそれに比べて規模がはるかに大きい。そのため手打ちで路線図を作っていくのは現実的ではない。そこで、SWA氏製作のLOP toolkit(
http://www.swa.gr.jp/lop/lop_a032.htmlにて公開
)に同梱された路線データedges.csvを読み込んで、2004年3月時点での最長経路を計算してみた。
edges.csvを読み込むスクリプトは以下の通り。
filename = 'edges.csv' f = codecs.open(filename, 'r', 'utf-8') lines = f.readlines() f.close() univ = [] univ2 = [] for line in lines: dat = line.split(',') dat[1] = int(dat[1]) dat[2] = int(dat[2]) dat[3] = float(dat[3]) edge = tuple(dat[1:4]) edge2 = tuple(dat[1:3]) if edge2 in univ2: print "Duplicate edge:", line err = 1 else: univ2.append(edge2) univ.append(edge) if err: print "Please check the duplicate edges and try again." exit() gs.set_universe(univ)
前に書いたように、Graphillionでは2つの頂点を結ぶ辺は1つしか扱えないので、重複があった場合は処理を中断させている。
edges.csvを読み込むと、重複する辺が5つ(北上→一ノ関,横浜→大船,福山→三原,三原→海田市,岩国→櫛ケ浜)出てくるので、対策をしておく。
そして、前回と同様のスクリプトで、中小国から小倉までのルートを数え上げ、最長のものを表示してみる。
さすがにグラフが大規模なだけあって、北海道のように一瞬で終わったりはしない。しかし、それでも手元のPCで、pathsメソッドを実行するのに45秒程度、最長のものを表示するのに10秒程度で済む。
SWA氏のLOP toolkitを使用してGLPK(v4.52)を用いると数分程度かかった。Graphillionを使うと解き直しの手間もかからないので、現状ではこの方法が最速と言えるのではないだろうか。
なお、得られた経路集合に対して、len()メソッドを用いてその総数を調べると、実に471,754,961,437,722,794,251,171,200通り(!)という結果が出る。
得られたルートを、edges.csvの情報をもとにわかりやすく表示してみる。
for eline in lines: dat = eline.strip().split(',') linename = dat[-2] staname = dat[-1] stacode1 = int(dat[1]) stacode2 = int(dat[2]) edge = tuple([stacode1,stacode2]) if edge in longest: print u'{0},{1}'.format(linename,staname).encode("utf-8")
出力は以下の通り。なんとなく眺めれば経路が浮かぶ・・・だろうか？？
津軽,中小国→青森 東北４,青森→野辺地 東北４,野辺地→八戸 東北新幹線,八戸→盛岡 山田,盛岡→茂市 山田,茂市→釜石 釜石,釜石→新花巻 釜石,新花巻→花巻 東北,花巻→北上 北上,北上→横手 東北,一ノ関→小牛田 大船渡,一ノ関→気仙沼 陸羽東,小牛田→古川 奥羽,秋田→大曲 羽越,秋田→余目 (以下略)
特例への対応
これを眺めると、今回探索した最長経路はSWA氏の解とやや異なっている。異なる点は、
の2点である。ともに運賃計算上の特例に起因するものだが(詳細はhttp://www.swa.gr.jp/lop/lop_a014.htmlも参照)、せっかくなので、Graphillionで対処してみる。
まず、岩国-櫛ケ浜は、この区間で山陽本線を使用してはならないという制約を加える。Graphillionでは、グラフ集合のうち特定の頂点や辺などの要素を含まない部分集合を返すexcludingメソッドが存在するので、一行で済む。
paths2 = paths.excluding([(299, 300)]) # 山陽,岩国→櫛ケ浜
次に、九州への入り方だが、幡生-新下関と新下関-小倉の同時使用を禁止する。これには、「幡生-新下関と新下関-小倉を両方使うグラフの集合」を求め、元のグラフ集合から除外する。
hatabu = paths2.including([(314, 315)]) # 山陽,新下関→幡生 prohibit = hatabu.including([(314, 318)]) # 新幹線,新下関→小倉 paths3 = paths2 - prohibit
得られたpaths3を先ほどの表示スクリプトで眺めると、SWA氏の算出したルートに一致することが確認された。めでたしめでたし。
Graphillionで最長片道切符のルート探索(1) 稚内→中小国の最長ルートを求める
Windows版のPythonにGraphillionを入れられなくて四苦八苦しているのが嫌なので、Ubuntuの入った先代のノートPCを引っ張り出して本題である最長片道切符の話を始めることに。
基本的には以下の記事を参考にしていった。合わせて読んでいただきたい。
Home · takemaru/graphillion Wiki · GitHub
SWA氏流にいえば、片道切符のルートとしては、タイプL, タイプO, タイプPの3種類が考えられる。Grphillionでは、タイプLについてはGraphSet.pathsというメソッドによって、たとえば始点と終点を指定した経路については全パターンを数え上げることができる。タイプOは、GraphSet.cyclesを使用すれば同様にことが済む。厄介なのはタイプPで、これは個別の扱いが必要になってきそう。
まずは始めの一歩として、「稚内→中小国の最長ルート」を求めてみる。このルートは現在の最長片道切符における道内ルートと一致するはずである。厳密に言えば、始点は稚内だと決めつけることは出来ないのであるが、それはまた追って。
Graphillionでグラフを扱うには、まずuniverseを設定する。これは考えるべきグラフ全体の空間を教えてやる作業で、要は路線図をインプットするということ。
今回は北海道の路線図ということで、以下の様な(SWA氏のedges.csv風)ファイルを用意した。
1, 3, 125.6, 津軽海峡, 中小国→五稜郭 2, 3, 3.4, 函館, 函館→五稜郭 3, 4, 23.6, 函館, 五稜郭→大沼 4, 5, 22.5, 函館, 大沼→森 4, 32, 33.5, 函館支, 大沼→東森 32, 5, 1.8, 函館支, 東森→森 5, 6, 62.7, 函館, 森→長万部 ...
カンマ区切りで、路線の始点(のコード)、終点(のコード)、営業キロの順になっている。グラフの「重み付き辺」を列挙する形式である。それ以降はすべてコメントとなる。
なお、函館本線の大沼→森間は経路が2つあるが、Graphillionでは全く同一の始点・終点をもつ辺が複数あるとエラーが発生するため、わざと支線上の東森で辺を分割してある。
これを読み込み、universeを設定する。スクリプトは以下の通り。
from graphillion import GraphSet as gs filename = 'edges_hokkaido.txt' f = open(filename, 'r') lines = f.readlines() f.close() univ = [] # 路線図全体のリスト for line in lines: dat = line.split(',') dat[2] = float(dat[2]) edge = tuple(dat[0:3]) # tupleにしなければ辺として扱えない univ.append(edge) gs.set_universe(univ) # universeを設定
実は、始点と終点のコードについては整数である必要はなく、この場合も文字列としてそのまま扱っている。駅名をそのまま入れてしまったほうが、目で見て分かりやすいかもしれない。プログラム上の扱いは面倒だが。
ここまでいったら、いよいよGraphillionの本領発揮。「稚内から中小国までの経路をすべて数え上げ、営業キロの長い順に表示する」には、以下のようなスクリプトを使う。
start = 30 # 稚内 end = 1 # 中小国 paths = gs.paths(start, goal) len(paths) # 経路の総数を表示 for path in paths.max_iter(): path # 最長経路を表示 break # break しないと、複数のパスが降順で取り出される
以上。ねぇ簡単でしょ？ 実行も一瞬で済む。
これで、「あの」最長経路が表示される。ただし、辺の順序はたどる順ではないので一見読みづらいかも。また、2番めに長い経路や最短経路も、最後の3行付近をいじれば同じく一瞬で分かる。
タイプPの取り扱いとか、始点を稚内以外にする話はまた今度。
GraphillionをWindowsにインストールしようとしてる
GraphillionをWindowsのPython 2.7(amd64)にインストールしようとしたメモ。
> easy_install graphillion
とコマンドを打てばいいようだが、走らせてみたところ以下のようなエラーが出てビルドに失敗した。
File "C:\Python27\lib\distutils\msvc9compiler.py", line 299, in query_vcvarsall
raise ValueError(str(list(result.keys())))
ValueError: [u'path']ソースをダウンロードして、python setup.py buildのコマンドで自分でビルドしようとしても、やはり同じエラーで止まった。どうもVC9でのコンパイルに失敗しているようだ。
エラーの文で調べると、同じような場所でつまづいている先人がいた。
Windows 7 x64にMercurialからgitを操作できるhg-gitをインストールする - digital matter
msvc9compiler.pyにパッチを当てれば解決するらしい。実際に現在使っていたmsvc9compiler.pyをのぞくと、この修正は適用されていなかった。そこでパッチを当てる。さらに追加で、
KEY_BASE = r"Software\Microsoft\\"
と前方に追加しないと、
NameError: name 'KEY_BASE' is not defined
というエラーが出るので、KEY_BASEという変数が出てくるより前に入れる。
これでMSVC9を用いてビルドを試みると、以下のエラーが出る。
src\pygraphillion.cc(464) : error C3861: 'strtoll': identifier not found
src\pygraphillion.cc(482) : error C2039: 'data' : is not a member of 'std::vecto
r<_Ty>'
with
[
_Ty=char
]最初のエラーだが、これはMSVCのstdlib.hにstrtollという関数が無いためである。_strtoi64という関数が代替で使えるので、Graphillionのソースコード(src/pygraphillion.cc)を書き換える。
#if defined(_MSC_VER) #define strtoll _strtoi64 #endif
の3行を最初に加えれば大丈夫。
次のエラーだが、以下の部分が原因のようだ。
static PyObject* setset_len2(PySetsetObject* self, PyObject* args) { PyObject* obj = NULL; if (!PyArg_ParseTuple(args, "|O", &obj)) return NULL; if (obj == NULL || obj == Py_None) { string size = self->ss->size(); vector<char> buf; for (string::const_iterator c = size.begin(); c != size.end(); ++c) buf.push_back(*c); buf.push_back('\0'); return PyLong_FromString(buf.data(), NULL, 0); } else if (PyInt_Check(obj)) { int len = PyLong_AsLong(obj); RETURN_NEW_SETSET(self, self->ss->size(len)); } else { PyErr_SetString(PyExc_TypeError, "not int"); return NULL; } }
調べてみると、vector::dataという関数が、VC9では使えないのか・・・？
Windowsで環境構築するのはほんと骨折れます。Ubuntuならeasy_install一発で済むのに。でもOSを複数入れると、切り替えが面倒でなぁ・・・
[追記09/12] Anaconda使ったほうがよさげです。まだ動作してないけど。詳しくは続編へ。
WindowsでPythonにscipyを入れるとき
Pythonで数値計算を走らせるのに欠かせないScipyだが、Windowsに導入するとBLASのコンパイルあたりでつまづきやすい。最初Cygwinに何とかして入れようとしたが、コンフィグがどうにもうまくいかないため、Windows版のインストーラを使うことに。
Python Extension Packages for Windows - Christoph Gohlke
このサイトに、様々なパッケージのインストーラが非公式ではあるが用意されている。これらを実行していけばとくに何も考えなくてもよさそう。
[追記09/12]ここではhttps://www.python.org/downloads/で入手できるPythonを対象としているが、ここに書いてある面倒なことをしなくてもAnacondaを入れるのが近道そう。
