オープンデータを使ってみよう！-流山市の桜を最適なルートで回ろう！-

今回のテーマはオープンデータ！

最近、何かと話題のオープンデータ。今回は、オープンデータで遊んでみたよという記事です。

f:id:tomoshige_n:20140801221248p:plain

オープンデータと言えば、「慶應義塾大学 SFC研究所データビジネス創造・ラボ」と「アクセンチュア」が第2回データビジネス創造コンテスト | 慶應義塾大学 × アクセンチュアを開催するらしいので学生のみんなは応募してみては？＆一緒に応募してみたい人募集しています。

今回の解析で使うデータ

今回は「流山市」のデータを使わさせてもらいました。一応クリエイティブ・コモンズでデータを使うと必ずリンクを載せてくださいということなので、最初に載せておきます。
流山市オープンデータ一覧｜流山市

あと、一応参考までに上のデータ解析コンペに掲載されている「課題」もせっかくなので書いておきましょう。

流山市の課題

流山市は、「都心から一番近い森のまち」、「母になるなら流山市」というキャッチフレーズのとおり、緑豊かな住環境の中で、子育てに優しいまちを目指しており、若い世代の方々が続々とお越しになっていますが、全国の自治体にみられるように、若年層の投票率が比較的低い状況です。平成27年4月に予定される統一地方選挙を控え、この投票率を向上させるためには、どのような方策が考えられるでしょうか。若い皆さんのユニークなアイデアを是非お聞かせください。（第2回データビジネス創造コンテストより抜粋）

今回の解析

今回の解析は、上の課題と関係あるかわかりませんが、4月と言えば桜の季節！そして、流山市は桜が有名らしい（確証はないですが、オープンデータに名所が記録されてたし。。。）。気候も良く運動もしたくなる季節ということで、流山市にある桜スポットを全部巡る最短ルートを書いてみましょう！ということにチャレンジします。

解析の方針

解析の方針としては、とりあえずデータを取得して、Googleさんに頼ってスポット同士の距離（徒歩）を算出、その後に巡回セールスマン問題だ！と考えて最適ルートを計算（パッケージ頼り...）、最後に図示とします。

お世話になるRパッケージや、資料など

Rパッケージ

ggmap（地図を作成する）
TSP（巡回セールスマン問題を解く）
ggplot2（美しい描画を可能に！）

資料（完全に頼りっぱなし...）

R advent calendar201320140102分

早速解析スタート

まず最初に、macで解析を行っているので、urlから直接取得はできません（文字コードが違う！！！ってエラーが出ます）。なので、一度流山市のHPの桜の名所（オープンデータ）｜流山市からデータをダウンロードしてsakura.csvとして保存します。

さて、ここからはRで処理していきましょう。まず、データの確認ですね。（後ろにアクセスやら、その場所の説明などが入っていますが、今回は「徒歩」で全て回るので関係ありません。削除！

#データの読み込み
sakura=read.csv("sakura.csv")
#変数名の変更
colnames(sakura)=c("No","Name","Place","Latitude","Longitude","Access","Description")

  No         Name                    Place Latitude Longitude
1  1     利根運河        流山市西深井829他 35.91525  139.9029
2  2 におどり公園 流山市西深井字大堺1028-1 35.91320  139.8885
3  3       清瀧院         流山市名都借1024 35.84851  139.9375
4  4       東福寺         流山市鰭ヶ崎1303 35.84177  139.9081
5  5 総合運動公園       流山市野々下1-29-4 35.85803  139.9193
6  6     大宮公園         流山市平和台5-44 35.85480  139.9077

さて、地図を書くために必要なのは"Latitude","Longitude"です。さて、書いてきましょう。といっても、パッケージに頼ります。コマンドの説明は全て中に書いています。

#必要なライブラリ
#RJSONIO,rjson,RgoogleMaps,ggmap
library(ggmap)
#緯度（縦側）
lat=sakura$Latitude
#軽度（横側）
lon=sakura$Longitude
#図示したい範囲
loc=c(min(lon),min(lat),max(lon),max(lat))

#googleから地図データを取得（zoomの大きさを間違えるとこの後エラーがでます）
M=ggmap(get_map(location=loc,zoom=13,source="google"))+xlab("")+ylab("")
#桜があるスポットをgoogle mapの上に書きましょう
MP=M+geom_point(data=sakura,mapping=aes(x=Longitude,y=Latitude),size=7,colour="red")

ということで、場所にポイントをつけることができました。図はこんな感じになります。赤い点がある場所が桜のある場所です。ここを最適に回るルートを考えましょう。

f:id:tomoshige_n:20140802023808p:plain

最適なルートを見つけるためには、まず各地点間を歩いてどれぐらいの時間がかかるのかという情報が必要です。この情報を取るためには、google先生に頼ります。library(ggmap)の中には、「出発地点の緯度・経度」と「ゴール地点の緯度・経度」を与えると最適なルートを探してくれるrouteという関数があります。これを用いて、各地点の距離を計算しましょう。コマンドはこんな感じになります。

#距離・時間行列の作成
#どんな移動手段をつかうのか
travelmode="walking"
#地点数（データの行数に対応）
nr=nrow(sakura)
#各地点間の距離を保存する行列
time=matrix(0,nr,nr)
#各地点間の歩行距離を保存する行列
distance=matrix(0,nr,nr)

#実際の計算
for(i in 1:nr){
	x_i=sakura[i,]
	loc_i=c(x_i$Longitude,x_i$Latitude)
	for(j in (i+1):nr){
		x_j=sakura[j,]
		loc_j=c(x_j$Longitude,x_j$Latitude)
		route.ij=route(loc_i,loc_j,mode=travelmode,structure="route")
		time[i,j]=time[j,i]<-sum(route.ij$minutes,na.rm=T)
		distance[i,j]=distance[j,i]<-sum(route.ij$km,na.rm=T)
	}
}

これで、実際の地点間の距離がわかりました。結果を見てみましょう。単位は分です。

         [,1]    [,2]    [,3]    [,4]   [,5]   [,6]   [,7]    [,8]   [,9]  [,10]
 [1,]   0.000  20.650 116.233 114.400 94.233 96.267 92.433 112.633 81.083 38.167
 [2,]  20.650   0.000 121.667 107.283 94.250 88.133 83.950 109.033 88.233 26.383
 [3,] 116.233 121.667   0.000  42.550 32.917 42.650 48.017  33.900 48.817 95.633
 [4,] 114.400 107.283  42.550   0.000 38.933 23.517 24.167  11.267 69.667 83.767
 [5,]  94.233  94.250  32.917  38.933  0.000 23.017 23.700  36.383 41.500 68.900
 [6,]  96.267  88.133  42.650  23.517 23.017  0.000  6.100  23.100 53.967 64.733
 [7,]  92.433  83.950  48.017  24.167 23.700  6.100  0.000  26.100 54.283 60.600
 [8,] 112.633 109.033  33.900  11.267 36.383 23.100 26.100   0.000 66.300 83.333
 [9,]  81.083  88.233  48.817  69.667 41.500 53.967 54.283  66.300  0.000 61.717
[10,]  38.167  26.383  95.633  83.767 68.900 64.733 60.600  83.333 61.717  0.000

次に、地点間の距離も確認しましょう。単位はkmです。近い場所は、散歩コースに良さそうです。

       [,1]  [,2]  [,3]  [,4]  [,5]  [,6]  [,7]  [,8]  [,9] [,10]
 [1,] 0.000 1.740 9.577 9.403 7.764 7.940 7.627 9.294 6.649 3.184
 [2,] 1.740 0.000 9.986 8.802 7.702 7.242 6.909 8.978 7.204 2.191
 [3,] 9.577 9.986 0.000 3.453 2.645 3.449 3.851 2.764 3.898 7.850
 [4,] 9.403 8.802 3.453 0.000 3.164 1.927 1.959 0.962 5.663 6.857
 [5,] 7.764 7.702 2.645 3.164 0.000 1.886 1.914 3.008 3.377 5.646
 [6,] 7.940 7.242 3.449 1.927 1.886 0.000 0.464 1.876 4.386 5.297
 [7,] 7.627 6.909 3.851 1.959 1.914 0.464 0.000 2.168 4.414 4.984
 [8,] 9.294 8.978 2.764 0.962 3.008 1.876 2.168 0.000 5.330 6.767
 [9,] 6.649 7.204 3.898 5.663 3.377 4.386 4.414 5.330 0.000 5.114
[10,] 3.184 2.191 7.850 6.857 5.646 5.297 4.984 6.767 5.114 0.000

最適なルートを探す（巡回セールスマン問題）

さて、ここからは最適なルートを探すことになります。これは巡回セールスマン問題と呼ばれる問題で計算がとても難しいことで知られています。詳しく知りたい方はwikipediaを参照してください。最近はアルゴリズムも開発されているようなので、計算機の発展とともに解ける幅は広がりそうです。巡回セールスマン問題 - Wikipedia

さてRで巡回セールスマン問題を扱いたい場合には、library(TSP)を用います。TSP=traveling salesman problemですね。これを用いて、先ほど作成したtimeという行列を入れるだけで結果が出ます。。。便利ですね...中身知らないのに使っていいのかなという不満・不安もありますが、ここは頼ることにしましょう。

#TSP（巡回セールスマン問題を解くアルゴリズム）
library(TSP)
#計算結果
ans=solve_TSP(TSP(time))
#結果の表示
ans[1:10]
[1]  7  5  1  2 10  9  3  8  4  6

ということで、7番から順番に上の順番で回るのが最適なルートということになります。最後にこれを図示しましょう。

#巡回経路の図示（日本語だと文字化けするけど...まぁ表示させる方法もあります...忘れた笑）
for(i in 1:(nr-1)){
	h1=ans[i]
	h2=ans[i+1]
	route1=na.omit(route(c(lon[h1],lat[h1]),c(lon[h2],lat[h2]),mode=travelmode,structure="route"))
	print(route1)
	#
	M=M+geom_point(data=route1[1,],aes(x=lon,y=lat),size=10,colour="red",alpha=0.75)
	M=M+geom_text(data=route1[1,],aes(x=lon,y=lat),label=i,colour="white")
	M=M+geom_text(data=route1[1,],aes(x=lon,y=lat+0.0005),label=sakura$Name[i],colour="red",size=4)
	M=M+geom_path(data=route1,aes(x=lon,y=lat),size=1.5,colour="red",alpha=0.75)
	M=M+geom_text(data=route1,aes(x=mean(lon),y=mean(lat)),label=paste(sprintf("%3.1f",sum(route1$minutes)),"min",sprintf("%3.1f",sum(route1$km)),"km",sep=""),colour="white",size=4)
	#if(nrow(route)<2){
	#	lonlat<-data.frame(lon=lon[c(h1,h2)],lat=lat[c(h1,h2)])
	#	M=M+geom_path(data=lonlat,aes(x=lon,y=lat),size=0.75,colour="red",alpha=0.75,type="dotted")
	#}
}
M=M+geom_point(data=route1[nrow(route1),],aes(x=lon,y=lat),size=10,colour="red",alpha=0.75)
M=M+geom_text(data=route1[nrow(route1),],aes(x=lon,y=lat),label=i+1,colour="white")
M=M+geom_text(data=route1[nrow(route1),],aes(x=lon,y=lat+0.0015),label=sakura$Name[h1],colour="red",size=4)

結果はこんな感じ。これが、google mapがはじき出した最適ルートです。
全体で28kmぐらいあることがわかります。と、ここまでが解析です。
f:id:tomoshige_n:20140802025510p:plain

感想とか

open dataを今回は可視化するということにチャレンジしてみました。いかがでしたか？
僕もあまりオープンデータをいじったことはなかったので、まだまだ知識がなく今回は上で紹介させてもらったslideに頼り切らせてもらうことになってしまいました。でも、やってみて結構勉強になったので、これからもちょこちょここういうことをやって行こうと思います。

そういえば、流山市は選挙の投票率（統一地方選挙）をしたいと言ってました。選挙、それ自体だけだと積極的にみんな行かないのでイベントと絡めるというのが良さそうですよね。統一地方選挙は４月ですから、例えば、桜の見物のついでに投票所に立ち寄れるような設計をすると投票率があがるかもしれませんね。そのためにはこの図の上に位置をプロットするのが効果的でしょう。それから人口の散布図なんかも更に書けると良さそうです。こういうのやってると、現地調査までやりたくなりますね、アンケート取ったり、インタビューした情報を突合するといい主張ができるかも。コンペでようかな...

それにしても、やはり、図示は強力な手法ですね。図示（ビジュアル）は統計との相性がやっぱり良くて、統計がわからない人にも、統計的なことを理解してもらうためのサポートになります。

告知的な

まだ時期はわからないんですけど、今年中にデータフェストっていうイベントを開催しようと思っています。データフェストはアメリカUCLA（カリフォルニア大学ロサンゼルス校）の統計学部が主体で開催している、学部生がデータ解析を勉強しつつ、実際の解析もやるっていうイベントです。

日本でもデータ解析のイベントたくさんありますが、ハードルが高いものが多いし、既にできる人向けなんですよね。そうじゃなくて、これから始める、始めたけどつまづいた人、やったことあるけど実際に使う機会がない人に向けて、データ解析のチュートリアル＋データ解析をやるイベントを作りたいなと考えています。データサイエンスできる人が必要だと言われて久しいですけど、あんまり初心者向けイベントないですし、僕らも専門家で閉じてても意味がなくて、知らない人に教える・伝えることを通して学ぶことがたくさんあるなって思うので。

詳しく決まってきたら、しっかり告知します！今回は、共有まで。

それでは、素敵な週末を。