鉄道路線の線形を可視化する
思いついたときに地図ネタを置いておく。
概要
よく「線形が良い・悪い」という言葉を聞くが、小さな縮尺の地図でなんとなく眺めていても分からない。そこで、急カーブの部分だけ色をわけて描画できないかと思い、トライしてみた。
データの取得・変換
まずはデータが無いと始まらない。
鉄道路線の座標データは、国土地理院の国土数値情報ダウンロードサービス（要登録：無料）から落とせる。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乗くらいの値を用いると急カーブから緩いカーブまでうまくグラデーションにできそうな感触。パラメータが大きすぎると急カーブばかり濃い色になる。
