読者です 読者をやめる 読者になる 読者になる

ill-identified diary

所属組織の業務や見解とは一切無関係なアフィリエイト付きメモ帳。所属とは関係ないけどここを見て所属先にも興味を持っていただけると喜びます。

[GIS] [R] 日本国内の鉄道網を可視化してみる

R GIS ggplot2 鉄オタ データハンドリング/データ加工 視覚化

今回やること

国土交通省国土数値情報ダウンロードサービスから鉄道の時系列データをダウンロードし, 国内の鉄道路線網がどう変わっていったかを, R を用いた処理方法を解説しつつ, 可視化してみる.

今回も R を使ってグラフを作成する. 以前の[R] Rで学ぶ都知事選のデータ可視化【地理データ編】で言及したように, 地理データのプロット専用のパッケージは使いづらいため, やはり ggplot2 を使用する. 国土数値情報の鉄道時系列データから, 平成25年度版のデータをダウンロードした.

使用する R のパッケージは maptools, ggplot2 の 2 つだけ*1である.

シェイプファイルの読み込み

まずは鉄道網の時系列データを読み込む方法を解説する. ダウンロードしてきたファイルは, 「シェイプファイル」とよばれる種類のファイルで .shp という拡張子がファイル名に付いている. readShapeLines() 関数は, 線分形式の地理データを読み込むときに使う. 今回は鉄道路線なので, この関数を用いるのが適当だろう. ダウンロードしたファイルを解凍すると複数のファイルが現れるが, 基本的に .shp ファイルを読みこめばいい.

proj4string=... という引数は, 読み込む地理データの測地系を指定するのに使う. これを正しく指定しないと, 実際の位置とずれる可能性がある. 国土数値情報のダウンロードページで座標系が JGD2000 であると明記されている. しかし, ここで指定できる引数は EPSG コードである. インフォマティックス社の 空間情報クラブ でこの辺の話が解説されており, また同ページで紹介されている EPSG Geodetic Parameter Registry で対応する EPSG コードを検索することができる. 結論から言うと, JGD2000 に対応する EPSG コードは 4612 であるので, 下記のようなコードで路線の地理データを読み込む (ファイルパスは任意で変更).

library(maptools)
library(ggplot2)

railHist <- readShapeLines('GIS/mlit/rail_hist_2013/N05-13_RailroadSection2.shp', proj4string=CRS("+init=epsg:4612"))

また, ダウンロードしたファイルの中には N05-13_Station2.shp というファイルも存在する. これは駅の名称や位置情報のデータであるが, 今回は使用しない.

国土数値情報ダウロードサービスのデータは全て Shift-JIS (CP932) なので, ウィンドウズ以外の環境では文字化けする可能性がある. そこで, 自分は以下の関数 sp.conv.932,utf8() を作成して, ファイルを全て UTF-8 に予め変換してある. ウィンドウズであれば, この処理は必要ない.

sp.conv.932.utf8 <- function(sp.data){
  for (i in colnames(sp.data@data)){
    sp.data@data[,i] <- iconv(sp.data@data[,i], 'CP932','UTF-8')
  }
  return(sp.data)
}
railHist <- sp.conv.932.utf8(railHist)

読み込んだオブジェクトは, sp というクラスになる. このクラスには, datalines といったスロットが存在する.

f:id:ill-identified:20140803233028p:plain

lines には線分の位置データが*2, data には, 位置情報を記録された地物の属性情報, つまり路線の名称や路線を所有する企業名などが記録されている. この data スロットはおなじみの data.frame 形式である. さきほどの方法で読み込んだ場合, この data の列名が N03_XXX というふうになっており, 一見何なのかわかりづらい. ダウンロードページに詳細が書かれているので, 以下のコードでわかりやすく改名している.

# 列名のベクトル
column.names.new <- c('type','rail.name','company','begin.service', 'begin','end', 'rail.id', 'transition.id', 'change.value', 'note', 'others', 'station.name')
colnames(railHist@data) <- column.names.new[1:11]

# 以下の3列は数値なので型変換
int.list <- c('begin.service','begin','end')
for ( i in int.list){
  railHist@data[,i] <- as.integer(railHist@data[,i])
}

下準備が終わったので, 路線データを年代別に分ける. 属性情報の begin に設置開始年が, end に開始終了年が4桁で与えられている. ただし, 1950年から 2013年までの情報しかなく, それ以前から開始している場合でも begin=1950 になっており, また, 不明な場合は 999 が与えられている. 一方, 現在も存続している場合は end=9999 になっている. よって, 簡単な条件式である時期に存在した路線を抽出することができる:

railHist.1960 <- railHist[railHist$begin <= 1960 & 1960 <= railHist$end, ]
railHist.1970 <- railHist[railHist$begin <= 1970 & 1970 <= railHist$end, ]
railHist.1980 <- railHist[railHist$begin <= 1980 & 1980 <= railHist$end, ]
railHist.1990 <- railHist[railHist$begin <= 1990 & 1990 <= railHist$end, ]
railHist.2000 <- railHist[railHist$begin <= 2000 & 2000 <= railHist$end, ]
railHist.2010 <- railHist[railHist$begin <= 2010 & 2010 <= railHist$end, ]
railHist.2013 <- railHist[railHist$begin <= 2013 & 2013 <= railHist$end, ]

これで, 1960, 1970, ... 2010 までの10年毎と, 2013年での国内の全路線図のデータがそれぞれ作成される. 開始・終了年には実は不明を意味する 999 が存在したが, 数件だけなので今回は無視することにした.

次に, ggplot でグラフ化できるようにデータを加工する. [R] Rで学ぶ都知事選のデータ可視化【地理データ編】 で言及したように, fortify() 関数を使って, 位置情報をデータフレームに変換してくれる. 今回は, データを抽出する目的以外で属性情報を使用しないが, fortify() 関数だと属性情報が切り捨てられるため, 位置と属性を紐付けたい場合, 別途に結合する必要があるので注意が必要である. 以下のコードでは, 7種類の条件で取り出したデータをリストでまとめている.

rail.df.ls <- list()
years <- seq(from = 1960, to=2010, by = 10)
for (i in 1:length(years) ){
   rail.df.ls[[i]] <- fortify(get(paste('railHist',years[i],sep='.')))
   
}
names(rail.df.ls) <- as.character(years)
rail.df.ls[[7]] <- fortify(railHist.2013)
names(rail.df.ls)[7] <- '2013'

データの活用方法

ここで, 他の抽出方法についても簡単に触れておく. rail.name(N05_002) には路線名が, company (N05_003) には運営会社が含まれており, これらを利用して条件付けをすることもできる. 例えば, 阪急線の全路線を取得する方法を考えてみよう. このデータに記載されている会社名がわからないため, grep() で部分一致検索をする.

unique(railHist@data$company[grep('阪急', railHist@data$company)])

結果はこうだった.

[1] "阪急電鉄"       "北大阪急行電鉄"

阪急以外にも, 北大"阪急"行電鉄がヒットしてしまったが, このデータでは 阪急は "阪急電鉄" と書かれていることがわかったので十分だ.

unique(railHist@data$rail.name[railHist@data$company=='阪急電鉄'])
 [1] "京都線"     "嵐山線"     "千里線"     "宝塚線"     "箕面線"     "神戸線"     "今津線"     "伊丹線"    
 [9] "甲陽線"     "神戸高速線"

となった. しかし, これは年に関する条件を置いてないので, 過去から現在までの全ての路線になる. よって, 今年のデータだけを出すように修正する.

unique(railHist@data$rail.name[railHist@data$company=='阪急電鉄' & railHist@data$begin <= 2013 & 2013 <= railHist@data$end ])
 [1] "京都線"     "嵐山線"     "千里線"     "宝塚線"     "箕面線"     "神戸線"     "今津線"     "伊丹線"    
 [9] "甲陽線"     "神戸高速線"

結果は変わらなかった. では, 1960年まで遡るとどうなるだろうか.

unique(railHist@data$rail.name[railHist@data$company=='阪急電鉄' & railHist@data$begin <= 1960 & 1960 <= railHist@data$end ])
character(0)

ゼロ件になってしまった. 実は, 阪急電鉄は1973年以降の社名で, それ以前は 「京阪神急行電鉄」だからだ*3. というわけで, そのような社名がデータに存在しないか再び確認.

unique(railHist@data$company[grep('京阪神', railHist@data$company)])
[1] "京阪神急行電鉄"

ということで, 再び発見できたので, 「京阪神急行電鉄」の社名で1960年に存在した路線名を改めて取得する.

unique(railHist@data$rail.name[railHist@data$company=='京阪神急行電鉄' & railHist@data$begin <= 1960 & 1960 <= railHist@data$end ])
 [1] "京都線" "十三線" "嵐山線" "千里線" "宝塚線" "箕面線" "神戸線" "今津線" "伊丹線" "甲陽線"

となり, 神戸高速線の統合前の路線リストとなった. ただ, ウィキペディアの記述だと60年時点では千里線は「千里山線」という名称だったようだが, ここでは「千里線」となっている (これに関してはこれまでの情報ではどちらが正しいかは断定できない).

プロット

また, グラフにするにあたって, 背景として日本の白地図を用意した. 海岸線や都道府県の境界データも国交省のサイトでダウンロードできるが, 県単位なので使いづらい. ポリゴンデータを結合するのは結構時間がかかるからである. そこで, ここでは GADM からダウンロードすることにした. なお, このサイトでは .shp のような一般的な地理データの形式だけでなく, .Rdata 形式でもダウンロードできる.

load('GIS/GADM/JPN_adm0.RData')
jpn.map <- fortify(gadm)

最後に, ggplot を用いてプロットする. 線分データは geom_path() 関数を使う. geom_polygon() 関数とかだと正しく表示されない. 今回は 1960, 1980年, 2013年の3つのみ結果を用意した*4.

g.base <- ggplot(data=jpn.map) + geom_path(aes(long,lat,group=group), color='blue') + theme_bw() + labs(x='',y='') + xlim(127, 147) + ylim(30,47) #+ xlim(125,150) + ylim(25,46)
g.base + geom_path(data=rail.df.ls[['1960']], aes(long,lat, group=group)) + labs(title='国内鉄道路線 (1960年)')
g.base + geom_path(data=rail.df.ls[['1980']], aes(long,lat, group=group)) + labs(title='国内鉄道路線 (1980年)')
g.base + geom_path(data=rail.df.ls[['2013']], aes(long,lat, group=group)) + labs(title='国内鉄道路線 (2013年)')

f:id:ill-identified:20140803233126p:plain

f:id:ill-identified:20140803233134p:plain

f:id:ill-identified:20140803233142p:plain

北海道の路線が年々減ってきているのが分かる. これは, 旅客鉄道だけでなく, 鉱山鉄道の路線も含まれているからである.

しかし, 海岸線と重なっているところもあり, 任意に拡大して見ることが出来たほうが便利そうである. しかし, R でできるのはこの辺が限界である. 次回は他のソフトウェアとも連携し, より見やすくする方法を紹介する ---(すぐ更新するとは言ってない)----.

続きはこちら: [GIS] [R] 日本国内の鉄道網を可視化してみる (後編) - ill-identified diary

参考文献

*1:依存パッケージは自動で読み込まれるので数えていない

*2:今回は線分データを読み込んだので lines となった. データの種類によっては polygons という名前になる場合もある

*3:Wikipedia 見たただけでちゃんと裏を取ったわけではない

*4:出力に結構時間がかかる