RでTSPの練習

導入

このスライドの内容を参考に, RでTSPします. http://www.slideshare.net/hiroki84/r-advent-calendar201320140102

千葉県流山市のオープンデータトライアルにて, 公共施設の座標が公開されています. これを利用して, 流山市にある"平和台"を含む公園が10箇所について, この10箇所を巡回する最短路を, TSPを解くことで調べます.

データ

library(dplyr)
library(data.table)
## 流山市オープンデータトライアルからデータを取得
## http://www.city.nagareyama.chiba.jp/10763/index.html
parks = fread("shisetsu_kouen.csv") %>% as.data.frame
## 平和台をgrep
parksSelected = parks[grep("平和台", parks$名称), 1:7]
parksSelected %>% head
##     大分類 小分類          名称                所在地  緯度  経度
## 269  公園   公園 平和台1号公園      流山市平和台5-40 35.86 139.9
## 270  公園   公園 平和台2号公園      流山市平和台2-12 35.85 139.9
## 271  公園   公園 平和台3号公園       流山市平和台1-4 35.86 139.9
## 272  公園   公園 平和台4号公園      流山市平和台5-60 35.86 139.9
## 273  公園   公園 平和台5号公園     流山市平和台5-733 35.86 139.9
## 274  公園   公園 平和台6号公園 流山市平和台4-1640-37 35.85 139.9
##         電話番号
## 269 04-7150-6092
## 270 04-7150-6092
## 271 04-7150-6092
## 272 04-7150-6092
## 273 04-7150-6092
## 274 04-7150-6092
parksSelected %>% dim
## [1] 10  7

プロット

GoogleMap上にプロットします. ggmapパッケージを利用します.

library(ggmap)
## 公園の経度緯度情報
lat = parksSelected$緯度
lon = parksSelected$経度
## 地図を描く範囲
loc = c(min(lon), min(lat), max(lon), max(lat))
## 結果をプロット
## ggmapでなんやかんやする
library(ggplot2)
map = ggmap(get_map(location = loc, zoom=15, source="google"))
map = map + geom_point(data = parksSelected, aes(x=経度, y=緯度), size=7)
map

plot of chunk unnamed-chunk-2

点の位置と, 公園の位置は, 合っていないそうです. 緯度経度を与えても, その場所に点を打ってくれないみたいです. 緯度経度の値が正しいみたいですが, よく分かりません.

距離行列の作成

ggmapパッケージの, route関数を使います. 2点の緯度, 経度を与えることで, 移動の道のりを計算してくれます. 中身はよく分かりません. 最短距離?

## 格納する行列を作る
nr = parksSelected %>% nrow
time = matrix(0, nr, nr)
distance = matrix(0, nr, nr)

## GoogleMapsApiからなんやかんやで距離を計算する
## 時間がかかるので, 1回計算して結果をとっておく
# for(i in 1:nr){
#   parksSelected_i = parksSelected[i,]
#   loc_i = c(parksSelected_i$経度, parksSelected_i$緯度)
#   for(j in (i+1):nr){
#     parksSelected_j = parksSelected[j,]
#     loc_j = c(parksSelected_j$経度, parksSelected_j$緯度)
#     route_ij = route(loc_i, loc_j, mode="walking", structure="route")
#     time[i,j] = time[j,i] = 
#       sum(route_ij$minutes, na.rm=TRUE)
#     distance[i,j] = distance[j,i] = 
#       sum(route_ij$km, na.rm=TRUE)
#   }
# }

## 結果を保存したファイルを読み込む
time = fread("timeTable.csv") %>% as.matrix
distance = fread("distanceTable.csv") %>% as.matrix

TSPを実行する

距離行列ができたので, TSPを回します. TSPパッケージに, 距離行列を流すだけです.

library(TSP)
ans = solve_TSP(ATSP(time))
## 答えです
ans[1:10]
##  V9  V5  V4  V1  V3 V10  V2  V6  V7  V8 
##   9   5   4   1   3  10   2   6   7   8

結果を出力

結果を, 地図上にプロットします route関数で道順を調べて, geom_pathするだけです.

library(ggplot2)
map = ggmap(get_map(location = loc, zoom=15, source="google"))
#map = map + geom_point(data = parksSelected, aes(x=経度, y=緯度), size=7)

## 順番に道を描く
for(k in 1:(nr-1)){
  h1 = ans[k]
  h2 = ans[k+1]
  route = route(c(lon[h1], lat[h1]), c(lon[h2], lat[h2]), mode="walking", structure="route")
  map = map + geom_path(data=route, aes(x=lon, y=lat), size=1.5)
## 一本道の場合は, geom_pathが使えない
## 2点のgeom_pathをすればいい? 
## よく分からないので保留です
#   if(nrow(route<2)){
#     lonlat = data.frame(lon = c(lon[h1], lon[h2]), 
#                         lat = c(lat[h1], lat[h2]))
#     map = map + geom_path(data = lonlat, 
#                           aes(x=lon, y=lat), size=1.5, 
#                           colour = "red")  
#     print(h1)
#   }
}
## 結果
map