【ggplot2】札幌圏におけるセイコーマートの強さをボロノイ分割で把握(したかった)
Good morning.
This is Silver weeeeeeeek!!!
さて。
進学で北海道に移住してからかれこれ4年とちょっと。
デビューしたての頃はずっと思っていました。
「セイコーマートってなんだ？」
と。
セイコーマート
それは北海道最強を誇るコンビニです。
ローソンよりもセブンよりも、道内における圧倒的地位と人口カバー率を誇るコンビニです。
今回は
コンビニの出店位置を母点にしたボロノイ分割によって、
他のコンビニとの出店戦略の違いを把握出来たらいいな、みたいなのをやります。
首都圏だともうやられている方もたくさんいらっしゃるみたいです。
方向性
例のごとく、作業前の自らへの課題とか仮説とか。
・店舗住所を{rvest}で抜いてきてGoogleAPIに投げて緯度経度を取得したいよ
・最後はggplot2で書きたいよ
・【仮説】セイコーマートは中心部よりも郊外出店で広い商圏を稼いでる気がする
前処理
まず{rvest}でそれぞれ札幌圏のコンビニ店舗情報を取得します。
今回ここで文字コード問題でつっかかったりしていましたが、かの有名なr-wakarangで助けて頂きました...
Windows環境で文字化けてしまう方、encode="UTF8"の引数指定でだめな時はEncoding(df)<-"UTF-8"で直接変換するといけるみたいです
url<-"http://www.mapion.co.jp/phonebook/M02005007C/" page = read_html(url,encoding = "UTF-8") #ここでのencode指定は無視された accessx<-page %>% html_nodes(xpath='//td') %>% html_text() Encoding(accessz)<- "UTF-8" #強制Encode accessx #1ページ目 accessy #2ページ目 accessz #3ページ目 access_s<-c(accessx,accessy,accessz)#各ページの抽出結果を結合(セイコーマート) #結果参照 > head(access_s) [1] "北海道札幌市中央区円山西町４丁目７−２５" "011-616-7911" [3] "北海道札幌市中央区北１条西１９丁目２−１" "011-615-2626" [5] "北海道札幌市中央区南７条西５丁目２８９−５４" "011-532-8610" #電話番号が邪魔そうなので正規表現で抜きます。 #札幌圏の市外局番は011から始まっているので、前方一致「^011」で引っかかった列以外を抽出 Index<-grep("^(011)",access) access[-Index] #あとはゴミ行取ったり、ベクトルにしたりなんだりで完成したのがこんな感じ adress_s<-as.vector(as.matrix(adress_s)) > str(adress_s) chr [1:299] "北海道札幌市厚別区もみじ台東４丁目２-１２" "北海道札幌市厚別区厚別中央１条７丁目１７-３７" ... > tail(adress_s) [1] "北海道札幌市北区北２２条西４丁目２-１２" "北海道札幌市北区北２４条西１４丁目３-３" "北海道札幌市北区北２７条西２丁目１-２５" [4] "北海道札幌市北区北３１条西４丁目１-１" "北海道札幌市北区北３６条西２丁目２-１" "北海道札幌市北区北８条西５丁目１" #ローソンとセブンもそれぞれ作りました > str(adress_rr) #ローソン chr [1:242] "北海道札幌市厚別区厚別中央１条４丁目１-２５" "北海道札幌市厚別区厚別中央３条６丁目１-１" ... > str(adress_7) #セブン chr [1:296] "北海道札幌市厚別区厚別西４条２丁目４-１５" "北海道札幌市厚別区厚別中央２条１丁目２-１７" ...
セイコーマート、ローソン、セブンイレブン
札幌市内の店舗数はそれぞれ240～300店舗ほどと、それほど数に差があるようには見受けられません。
ちなみに全然関係無いんですけど、
> str(adress) chr [1:299] "北海道札幌市厚別区もみじ台東４丁目２－１２""| __truncated__ "北海道札幌市厚別区もみじ台北５丁目１－６""| __truncated__ ...
こういう"|_truncated"とかいう邪魔者をスマートに消す方法があったら知りたいです...そもそもだれなのこれ...
GoogleAPIから緯度経度を取得
文字列ベクトル(3店舗合計900件)の座標を取得します。
GoogleAPIはJSONで返してくれるので、一件数投げて、緯度と経度だけを取得するならこんな感じで
{Rcurl}と{rjson}を使う感じになるのでしょうか。
address <- "住所文字列" url <- paste("http://maps.googleapis.com/maps/api/geocode/json?address=", address,"&sensor=false®ion=JP&language=ja",sep="") # Geocoding APIからJSON形式で色々返してもらう json<-RCurl::getURL(url) # JSONからRに緯度・経度の情報を抽出してオブジェクトに格納 rjson <- rjson::fromJSON(json) lat <- rjson$results[[1]]$geometry$location$lat lon <- rjson$results[[1]]$geometry$location$lng
しっかし900件あるしな～forか関数作るかかなあ～という感じですが、
@teramonagi先生が複数対応の関数を作ってくださっていたので、一部お借りすることに。
#緯度経度情報へ変換(teramonagi先生関数) location_s <- GetLatitudeAndLongitude(adress_s) location_r <- GetLatitudeAndLongitude(adress_r) location_7 <- GetLatitudeAndLongitude(adress_7) #3店舗それぞれに店の種別「store」を足して別オブジェクトへ格納 > location_s %>>% + (dplyr::mutate(.,store="セイコーマート")->location_ss) %>>% + head(3) lat lon store 1 43.03022 141.4996 セイコーマート 2 43.03790 141.4800 セイコーマート 3 43.03768 141.4663 セイコーマート > location_r %>>% + (dplyr::mutate(.,store="ローソン")->location_rr) %>>% + head(3) lat lon store 1 43.03524 141.4687 ローソン 2 43.04238 141.4731 ローソン 3 43.03987 141.4886 ローソン > location_7 %>>% + (dplyr::mutate(.,store="セブンイレブン")->location_77) %>% + head(3) lat lon store 1 43.04993 141.4590 セブンイレブン 2 43.03077 141.4536 セブンイレブン 3 43.03966 141.4577 セブンイレブン #3店舗をマージ > location_ss%>>% + (rbind(.,location_rr)) %>>% + (rbind(.,location_77)->location) %>>% + tail(3) lat lon store 835 43.10258 141.3402 セブンイレブン 836 43.10580 141.3370 セブンイレブン 837 43.06910 141.3480 セブンイレブン > location$store<-as.factor(location$store) > str(location) #マージ完成 'data.frame': 837 obs. of 3 variables: $ lat : num 43 43 43 43 43 ... $ lon : num 141 141 141 141 141 ... $ store: Factor w/ 3 levels "セイコーマート",..: 1 1 1 1 1 1 1 1 1 1 ...
できた。計837件だそうです。
storeは図にするときに色分けのため、factor化してます。
ちょっと可視化
長き戦いを終えたので、ここでいったん地図に落とし込みます。
#札幌市の地図オブジェクトを取得(Zoom引数でちょっと調整) sapporo_center <- geocode("Sapporo Station", source = "google") #札幌駅の位置を取得 sapporo_map <- ggmap::get_googlemap(center = c(lon =sapporo_center[1,1] , lat = sapporo_center[1,2]), zoom = 12, scale = 2,maptype = "roadmap",color ="color") #satellite', 'roadmap' #コンビニマップを描画 conmap2<-ggmap::ggmap(sapporo_map)+geom_point(data=location,aes(x=lon,y=lat,colour=store),size=2)
できたぜ。
これだけでも中心部にローソンとセブンがわーーーっと集まってて、
セイコーマートは店舗を集中させずにまんべんなく出店させてそう...？みたいなのが見える気がする。
ボロノイ分割
さて。メインテーマであるボロノイ分割をやります。
真剣に説明しようとおもったけど、図を改めて作るのがアレだったので...
分かりやすい記事がありましたのでこちらを...
すっごい端的に言うと、
A点-------------------------B点
ここの間に垂直二等分線を引きます。
C点とD点、E点と...みたいにドロネー線に対してどんどん垂直二等分線(ボロノイ線)を引いていきます。
最後にドロネー線を消したら、各点を中心とした「領域」が出来るよね！って感じです。
やります。
library(deldir) #deldir::分割 deldir <- deldir::deldir(x=location$lon,y=location$lat, rw=c(min(location$lon),max(location$lon),min(location$lat),max(location$lat))) #すごいいっぱいのリストが出るよ
rw引数は分割するはじっこの線を決める感じで、この枠内を分割してねみたいな書き方をしてみました。
結果を表示したかったけれどもブログが長くなりすぎそうなので主要なものを以下に。
delsgs:ドロネー三角形。エッジの座標点（x1,y1,x2,y2）の後ろに母点(コンビニ店舗)のID
dirsgs:ディリクレ領域（ボロノイ領域）。エッジの座標点（x1,y1,x2,y2）の後ろに母点(コンビニ店舗)のID
n.data:分割されたデータ点(重複なし)
n.dum:分割されたデータ点(重複点なし)
del.area:ドロネー三角に分割された点の面積。
dir.area:ボロノイ境界領域によって囲まれている面積。
rw:正方形領域の枠。
....etc(もっといろいろあります)
今回はボロノイ領域を書きたいので、dirsgsの結果を使います。
geom_segment()でごにょごにょと。
conmap3<-ggmap::ggmap(sapporo_map) + geom_point(aes(x=lon, y=lat,colour=factor(store)),size=2, data = location) + geom_segment(aes(x = x1, y = y1, xend = x2, yend = y2), size = .25,data = deldir$dirsgs, linetype = 1);conmap3#境界が書ける
見えない。見えにくいよパトラッシュ。
元の地図テーマを変えてみた。
なんかこわい。まあいいか。
ボロノイ境界を引いてみてわかったことは、特に北東部、南西部でセイコーマートのボロノイ領域が「広い」
西側はなかなかセブンも頑張っていますが、やはり市内の端にいけばいくほどセイコーマートの存在感が増しています。
ローソンの店舗は地下鉄南北線、東西線に沿ってる感じですがね、人が多いところに出店しまくってる感じがします。
セブンは人が多いところを中心に、たまに穴場的なところへの出店。
セイコーマートは中心部はやや避け、とにかくまんべんなく、人口＆土地カバー率といったところでしょうか。
見にくいから色を塗りたい
ここからはうまくいっていない編です。
geom_polygon()を使って母点(コンビニ)の色で領域を色分けしたら見やすいよなあって感じなんですが。。
どうやってtileと母点のFactorを結びつけるか方針が決まりませんでした...
#tile情報だけを切り出してきて、その中身を色塗り tilelist <- unclass(tile.list(deldir)) tilelist <- lapply(tilelist, function(l){ data.frame(x = l$x, y = l$y) }) library(reshape2) tiledf <- reshape2::melt(tilelist, id.vars = c("x")) tiledf <- tiledf[,c(4,1,3)] names(tiledf) <- c('tile','x','y') tiledf$tile<-as.factor(tiledf$tile) #各tileの頂点の座標が取得出来てます。ここに母点情報を足せればいいのかなあ... > head(tiledf) tile x y 1 1 142.1004 43.17976 2 1 142.0103 43.27621 3 1 141.5015 43.03878 4 1 141.4951 43.03427 5 1 141.4888 43.02622 6 1 141.4889 43.02588 #描画 conmap3<-conmap3+geom_polygon(aes(x = x, y = y, fill = tile),data = tiledf, colour = 'black', alpha = .3)
ちがうそうじゃない....ってなる....
現状邪魔なラベルを消す元気もないです。
うーーーむ。
以上、現場からでした。
