今回やること

国土交通省の国土数値情報ダウンロードサービスから鉄道の時系列データをダウンロードし, 国内の鉄道路線網がどう変わっていったかを, 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 というクラスになる. このクラスには, data や lines といったスロットが存在する.

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年)')