【orloca,leaflet】ラーメン屋さんになりたい ~SappoRoR#6反省会~
先日SappoRoRへ行ってきて、
おもむろに「ラーメン屋さんになりたい」と言ってきました。
田舎の両親、ラーメン屋さんのみなさま、アホな娘でごめんなさい。
あまり後悔はしてません。
資料各種
- 改行がダサかったので上げ直し。
- コード類 github.com
概要
テーマ
- ラーメン屋さんになりたい
- 札幌市中央区で最も立地に適した場所を探す
データ準備
e-stat(国勢調査平成22年/小地域)より、札幌市中央区のシェープファイルを取得。
いつもの流れ。
getinfo.shape(shape.file) ## Shapefile type: Polygon, (5), # of Shapes: 990 pj<-CRS("+proj=longlat +datum=WGS84") Chuo<- readShapePoly(shape.file,proj4string = pj) pj<-CRS("+init=epsg:2454") CRSargs(pj) ## [1] "+init=epsg:2454 +proj=tmerc +lat_0=44 +lon_0=142.25 +k=0.9999 +x_0=0 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs" Chuo2<- readShapePoly(shape.file,proj4string = pj) str(Chuo@data) ## 'data.frame': 990 obs. of 38 variables: ## $ AREA : num 34241 23602 65333 450562 13707 ... ## $ PERIMETER : num 755 649 1215 2664 536 ... ## $ H22KA01_ : int 13586 13677 13726 13830 13848 13951 14018 14067 14102 14140 ... ## $ H22KA01_ID: int 13585 13676 13725 13829 13847 13950 14017 14066 14101 14139 ... ## $ KEN : Factor w/ 1 level "01": 1 1 1 1 1 1 1 1 1 1 ... ## $ CITY : Factor w/ 1 level "101": 1 1 1 1 1 1 1 1 1 1 ... ## $ KEN_NAME : Factor w/ 1 level "北海道": 1 1 1 1 1 1 1 1 1 1 ... ## $ SITYO_NAME: Factor w/ 1 level "石狩振興局": 1 1 1 1 1 1 1 1 1 1 ... ## $ GST_NAME : Factor w/ 1 level "札幌市": 1 1 1 1 1 1 1 1 1 1 ... ## $ CSS_NAME : Factor w/ 1 level "中央区": 1 1 1 1 1 1 1 1 1 1 ... ## $ HCODE : int 8101 8101 8101 8101 8101 8101 8101 8101 8101 8101 ... ## $ KIHON1 : Factor w/ 81 levels "0200","0300",..: 81 80 79 76 78 77 76 75 76 74 ... ## $ DUMMY1 : Factor w/ 1 level "-": 1 1 1 1 1 1 1 1 1 1 ... ## $ KIHON2 : Factor w/ 31 levels "00","01","02",..: 16 16 16 17 16 16 16 16 22 16 ... ## $ KEYCODE1 : Factor w/ 987 levels "101020000","101030000",..: 987 986 985 980 984 983 979 977 982 974 ... ## $ KEYCODE2 : Factor w/ 987 levels "1010200","1010300",..: 987 986 985 980 984 983 979 977 982 974 ... ## $ AREA_MAX_F: Factor w/ 1 level "M": 1 1 1 1 1 1 1 1 1 1 ... ## $ KIGO_D : Factor w/ 0 levels: NA NA NA NA NA NA NA NA NA NA ... ## $ N_KEN : Factor w/ 0 levels: NA NA NA NA NA NA NA NA NA NA ... ## $ N_CITY : Factor w/ 0 levels: NA NA NA NA NA NA NA NA NA NA ... ## $ N_C1 : int 0 0 0 0 0 0 0 0 0 0 ... ## $ KIGO_E : Factor w/ 3 levels "E1","E2","E3": NA NA NA NA NA NA NA NA NA NA ... ## $ KIGO_I : Factor w/ 0 levels: NA NA NA NA NA NA NA NA NA NA ... ## $ TATE : int 0 0 0 0 0 0 0 0 0 0 ... ## $ DIR : int 0 0 0 0 0 0 0 0 0 0 ... ## $ HIGHT : int 50 50 50 50 50 50 50 50 50 50 ... ## $ JIKAKU : int 10 10 10 10 10 10 10 10 10 10 ... ## $ NMOJI : int 6 6 5 9 5 5 9 9 9 9 ... ## $ MOJI : Factor w/ 987 levels "旭ケ丘1丁目",..: 905 903 904 900 898 879 899 869 902 876 ... ## $ SEQ_NO2 : int 13585 13676 13725 13829 13847 13950 14017 14066 14101 14139 ... ## $ KSUM : int 2 1 8 1 3 3 3 7 2 9 ... ## $ CSUM : int 2 5 8 1 3 3 3 7 2 12 ... ## $ JINKO : int 29 518 206 0 136 149 142 285 41 699 ... ## $ SETAI : int 1 304 108 0 72 72 80 197 20 307 ... ## $ X_CODE : num 141 141 141 141 141 ... ## $ Y_CODE : num 43.1 43.1 43.1 43.1 43.1 ... ## $ KCODE1 : Factor w/ 987 levels "0200-00","0300-00",..: 987 986 985 980 984 983 979 977 982 974 ... ## $ KEY_CODE : Factor w/ 987 levels "011010200","011010300",..: 987 986 985 980 984 983 979 977 982 974 ... ## - attr(*, "data_types")= chr "N" "N" "N" "N" ... Chuo.dc<- data.frame(MOJI=Chuo@data$MOJI, X_CODE=Chuo@data$X_CODE, Y_CODE=Chuo@data$Y_CODE, KEY_CODE=Chuo@data$KEY_CODE, JINKO=Chuo@data$JINKO, SETAI=Chuo@data$SETAI, X_CEN=coordinates(Chuo)[,1], Y_CEN=coordinates(Chuo)[,2]) knitr::kable(head(Chuo.dc))
MOJI | X_CODE | Y_CODE | KEY_CODE | JINKO | SETAI | X_CEN | Y_CEN | |
---|---|---|---|---|---|---|---|---|
0 | 北二十二条西 | 141.3284 | 43.08521 | 011017922 | 29 | 1 | 141.3284 | 43.08521 |
1 | 北二十一条西 | 141.3285 | 43.08400 | 011017921 | 518 | 304 | 141.3285 | 43.08400 |
2 | 北二十条西 | 141.3280 | 43.08234 | 011017920 | 206 | 108 | 141.3280 | 43.08234 |
3 | 北十六条西16丁目 | 141.3260 | 43.07778 | 01101791616 | 0 | 0 | 141.3260 | 43.07778 |
4 | 北十八条西 | 141.3298 | 43.08034 | 011017918 | 136 | 72 | 141.3299 | 43.08034 |
5 | 北十七条西 | 141.3302 | 43.07896 | 011017917 | 149 | 72 | 141.3301 | 43.07896 |
人口コロプレス図で全体を把握
編集したデータを人口コロプレス図にして描画。
ここまではいつもの流れ。
pal <- colorNumeric( palette = "YlOrRd", domain = Chuo.dc$JINKO ) leaflet::leaflet(Chuo) %>% addProviderTiles(provider = "OpenStreetMap.BlackAndWhite") %>% addPolygons(stroke = FALSE, fillOpacity = 0.9, color = ~pal(Chuo.dc$JINKO))
各ポリゴンの代表点を抽出して「世帯」と見立てる
いつもの如く、ワタシのような一般ピーポーには個別世帯のデータなど手に入らないので
各ポリゴンの代表点を取得し、「世帯」とする。
コロプレス図に描画した各小地域の人口を「各世帯の世帯人口」としてあてはめ
leaflet::leaflet(Chuo) %>% addProviderTiles(provider = "OpenStreetMap.BlackAndWhite") %>% addPolygons(stroke = FALSE, fillOpacity = 0.9, color = ~pal(Chuo.dc$JINKO)) %>% addCircles(lng = coordinates(Chuo)[,1],lat=coordinates(Chuo)[,2],color = "gray80")
立地分析(Minsum問題)
先ほど作った「世帯」を「世帯人口」で重み付けし、
移動コストが最小となる地点を探す。
今回のような施設数=1(ラーメン屋)となる算出については
orlocaパッケージのzsummin関数がよさそう。
loca<-loca.p(x=Chuo.dc$X_CEN,y=Chuo.dc$Y_CEN,w=Chuo.dc$JINKO) (point<-zsummin(loca)) ## [1] 141.33469 43.05066 x<- as.data.frame(loca@x) y<- as.data.frame(loca@y) w<- as.data.frame(loca@w) knitr::kable(head(data.frame(x=x,y=y,w=w)))
loca.x | loca.y | loca.w |
---|---|---|
141.3284 | 43.08521 | 29 |
141.3285 | 43.08400 | 518 |
141.3280 | 43.08234 | 206 |
141.3260 | 43.07778 | 0 |
141.3299 | 43.08034 | 136 |
141.3301 | 43.07896 | 149 |
算出されたラーメン屋さんの位置を可視化
おもむろにallows関数で描いてみる。
plot(Chuo) for(i in 1:nrow(x)){ arrows(x[i,],y[i,],point[1],point[2],col="red") } arrows(x[1,],y[1,],point[1],point[2],col="red",lwd=0.2)
こわすぎ。
leafletでも可視化したい
小地域ごとの線が見栄え上邪魔なので、ポリゴンを結合して
札幌市中央区の「外側の線」を作る
#ポリゴンの結合------------------------- maptools::gpclibPermit() ## [1] FALSE Chuo.union<- maptools::unionSpatialPolygons(Chuo,IDs=rep(1,times=nrow(Chuo@data))) str(Chuo.union) ## Formal class 'SpatialPolygons' [package "sp"] with 4 slots ## ..@ polygons :List of 1 ## .. ..$ :Formal class 'Polygons' [package "sp"] with 5 slots ## .. .. .. ..@ Polygons :List of 1 ## .. .. .. .. ..$ :Formal class 'Polygon' [package "sp"] with 5 slots ## .. .. .. .. .. .. ..@ labpt : num [1:2] 141 43 ## .. .. .. .. .. .. ..@ area : num 0.00513 ## .. .. .. .. .. .. ..@ hole : logi FALSE ## .. .. .. .. .. .. ..@ ringDir: int 1 ## .. .. .. .. .. .. ..@ coords : num [1:1267, 1:2] 141 141 141 141 141 ... ## .. .. .. .. .. .. .. ..- attr(*, "dimnames")=List of 2 ## .. .. .. .. .. .. .. .. ..$ : NULL ## .. .. .. .. .. .. .. .. ..$ : chr [1:2] "x" "y" ## .. .. .. ..@ plotOrder: int 1 ## .. .. .. ..@ labpt : num [1:2] 141 43 ## .. .. .. ..@ ID : chr "1" ## .. .. .. ..@ area : num 0.00513 ## ..@ plotOrder : int 1 ## ..@ bbox : num [1:2, 1:2] 141.2 43 141.4 43.1 ## .. ..- attr(*, "dimnames")=List of 2 ## .. .. ..$ : chr [1:2] "x" "y" ## .. .. ..$ : chr [1:2] "min" "max" ## ..@ proj4string:Formal class 'CRS' [package "sp"] with 1 slot ## .. .. ..@ projargs: chr "+proj=longlat +datum=WGS84 +ellps=WGS84 +towgs84=0,0,0" #描画ワクだけ (result<-leaflet::leaflet(Chuo.union) %>% addProviderTiles(provider = "OpenStreetMap.Mapnik") %>% addPolygons(stroke = TRUE,color="gray20"))
改めて結果を可視化
ここがあまりきれいなコードじゃなくて納得いってない..... 各世帯の座標→ラーメン屋 の線をそれぞれ世帯数分引いている。なんかいいやり方はないものか。
#-立地分析leaflet-------------------------------------------------- loca<-loca.p(x=Chuo.dc$X_CEN,y=Chuo.dc$Y_CEN,w=Chuo.dc$JINKO) point<-zsummin(loca) x<- as.data.frame(loca@x) y<- as.data.frame(loca@y) w<- as.data.frame(loca@w) resultdata<- data.frame(x=x[1:990,], y=y[1:990,], pointx=rep(point[1],990), pointy=rep(point[2],990)) knitr::kable(tail(resultdata))
x | y | pointx | pointy | |
---|---|---|---|---|
985 | 141.3443 | 43.02096 | 141.3347 | 43.05066 |
986 | 141.3471 | 43.02086 | 141.3347 | 43.05066 |
987 | 141.3488 | 43.02091 | 141.3347 | 43.05066 |
988 | 141.3486 | 43.01963 | 141.3347 | 43.05066 |
989 | 141.3473 | 43.01900 | 141.3347 | 43.05066 |
990 | 141.3455 | 43.01864 | 141.3347 | 43.05066 |
test<-result for(i in 1:nrow(resultdata)){ testd<-data.frame(x = as.numeric(resultdata[i, c(1, 3)]), y = as.numeric(resultdata[i, c(2, 4)])) test <- test %>% addPolylines(data=testd, lng = ~x, lat = ~y, color="#f30", weight="2") } test %>% setView(lng = point[1], lat = point[2], zoom = 13) %>% addPopups(lng=point[1],lat=point[2],popup="here!!!")
まとめ
算出した「here!」に居た先客とは!!
...というくだりには触れないでおく...札幌ネタ....
LTにはちょうどいいくらいの内容だったはずなのに、
他のみなさまがまさかのrevealjsスライド祭りで大変な汗をかきました。
次は10月なのかな。お待ちしてます。