次元の海で溺れる

Rとデータ解析と統計手法たちとわたし

【sf,mapview】ラーメン屋さんになるにはまだ商圏人口についての理解が足りない

そういえば1年ぐらい前に試される大地を離れて関東に引っ越しました。
文明の利器「えあこん」は素晴らしい。もうあなた無しでは生きられない。

前段

`sf`パッケージなるものがあるらしい、という話を聞いたのはいつのことだったか。
詳しく書いてくださるみなさんの記事を見ているとなるほど気になる....

uribo.hatenablog.com

notchained.hatenablog.com

ksmzn.hatenablog.com

テーマ

思い返せば私はラーメン屋さんになりたかった気がする。
日々のドタバタにかまけてそんなことも忘れていた。

wafdata.hatenablog.com

しかしながら明らかに過去の検討は甘い。
人口で重みづけした重心を求めるだけで立地を決めてよいのか。
もっと集客可能な商圏内についての考察が必要なのではなかろうか。
真剣みが足りない。

今日のゴール

  1. 店舗の付近円商圏(5km,10km)の地域を特定し、各商圏人口を概算する
  2. sfに慣れる
  3. mapviewに慣れる

関連データ

  1. 千葉県船橋市の小地域シェープファイル(世界測地系・緯度経度 shape形式)

 ※ 地図で見る統計(統計GIS) | 政府統計の総合窓口 平成22年国勢調査


船橋にしたのは偶然です。決して最近船橋のだし茶漬け屋さんにハマったとかではありません断じて。

データの読み込み

いつもは`CRS`諸々指定して読み込むところですが、今回は`sf`なので恐る恐る。

library(dplyr)
library(ggplot2)
library(sf)
library(mapview)
library(units)
raw <-sf::read_sf('./data/A002005212015DDSWC12204/h27ka12204.shp',
            options="encoding=cp932")
raw %>% head()

>Simple feature collection with 6 features and 35 fields
geometry type:  POLYGON
dimension:      XY
bbox:           xmin: 139.9887 ymin: 35.689 xmax: 140.0012 ymax: 35.70095
epsg (SRID):    NA
proj4string:    +proj=longlat +ellps=GRS80 +no_defs
     KEY_CODE PREF CITY S_AREA PREF_NAME CITY_NAME     S_NAME KIGO_E HCODE      AREA PERIMETER H27KAxx_
1 12204001001   12  204 001001    千葉県    船橋市 宮本1丁目   <NA>  8101  95071.29  1386.969     2856
2 12204001002   12  204 001002    千葉県    船橋市 宮本2丁目   <NA>  8101  98592.73  1453.452     3810
3 12204001003   12  204 001003    千葉県    船橋市 宮本3丁目   <NA>  8101  88038.99  1357.873     1045
4 12204001004   12  204 001004    千葉県    船橋市 宮本4丁目   <NA>  8101 103077.04  1623.495     3336
5 12204001005   12  204 001005    千葉県    船橋市 宮本5丁目   <NA>  8101  91841.85  1284.460     2095
6 12204001006   12  204 001006    千葉県    船橋市 宮本6丁目   <NA>  8101 222729.01  2134.591     5262
  H27KAxx_ID KEN KEN_NAME SITYO_NAME GST_NAME CSS_NAME KIHON1 DUMMY1 KIHON2  KEYCODE1  KEYCODE2 AREA_MAX_F
1       2855  12   千葉県       <NA>   船橋市     <NA>   0010      -     01 204001001 204001001          M
2       3809  12   千葉県       <NA>   船橋市     <NA>   0010      -     02 204001002 204001002          M
3       1044  12   千葉県       <NA>   船橋市     <NA>   0010      -     03 204001003 204001003          M
4       3335  12   千葉県       <NA>   船橋市     <NA>   0010      -     04 204001004 204001004          M
5       2094  12   千葉県       <NA>   船橋市     <NA>   0010      -     05 204001005 204001005          M
6       5261  12   千葉県       <NA>   船橋市     <NA>   0010      -     06 204001006 204001006          M
  KIGO_D N_KEN N_CITY KIGO_I       MOJI KBSUM JINKO SETAI   X_CODE   Y_CODE  KCODE1
1   <NA>  <NA>   <NA>   <NA> 宮本1丁目    34  1849  1057 139.9920 35.69872 0010-01
2   <NA>  <NA>   <NA>   <NA> 宮本2丁目    23  1717   963 139.9906 35.69405 0010-02
3   <NA>  <NA>   <NA>   <NA> 宮本3丁目    16  1365   738 139.9934 35.69077 0010-03
4   <NA>  <NA>   <NA>   <NA> 宮本4丁目    21  1216   632 139.9936 35.69238 0010-04
5   <NA>  <NA>   <NA>   <NA> 宮本5丁目    13   805   388 139.9936 35.69508 0010-05
6   <NA>  <NA>   <NA>   <NA> 宮本6丁目    40  2849  1330 139.9963 35.69858 0010-06
                        geometry
1 POLYGON ((139.9926 35.69725...
2 POLYGON ((139.9912 35.69277...
3 POLYGON ((139.9955 35.69058...
4 POLYGON ((139.9971 35.69015...
5 POLYGON ((139.9947 35.69422...
6 POLYGON ((139.996 35.6969, ...

geometeryのところに各ポリゴンの座標が入ってる形。
ざっとmapviewで可視化してみる。

raw %>% 
mapview::mapview(label=row.names(.))

f:id:WAFkw:20180724234421p:plain

何このロンギヌスの槍形状。かっこいい。

とはいえこの鋭利な刃物は海。我がラーメン屋の商圏としてはいかがか、という感じなので
要素指定で削り、ベースマップとして整理する。

base_map <- raw[-337,]
base_map %>% 
mapview::mapview(label=row.names(.))

f:id:WAFkw:20180724234830p:plain

いい感じ。

店舗位置を中心とした円商圏の設定

まずラーメン屋の座標を設定する。

ここをキャンプ地とする、と私が言ったらそこがラーメン屋なのだ。
船橋駅の座標とぴったり一致してる?...そんなことは知らない。

store_point <- st_point(c(139.985338,35.701557))

mapview::mapview(base_map,label=row.names(base_map))+
  mapview::mapview(store_point,color="red")

f:id:WAFkw:20180724235319p:plain

建設完了。

あとはここを中心に5km,10kmの円商圏を設定する。
ちょっと悩んだけどsf::st_buffer()で円ポリゴンを作ることにした。

# kmの単位で取り扱いできるよう設定
r1 <- set_units(5,km)
r2 <- set_units(10,km)

# そのままでもst_bufferを使えるが、kmがうまく使えなかったり楕円になってしまう。
# 一度平面直角座標系に変換してからbufferを取る→戻す、という戦略で進めてみる

store_point<-st_sfc(store_point,crs=4612)
store_point<-st_transform(store_point,2451)

market_5km <- sf::st_buffer(store_point,dist=r1)
market_10km <- sf::st_buffer(store_point,dist=r2)

store_point<-st_transform(store_point,4612)
market_5km<-st_transform(market_5km,4612)
market_10km<-st_transform(market_10km,4612)

mapview(base_map)+
  mapview(store_point,color="red")+
  mapview(market_10km,lwd=3,alpha.regions=0.0,layer.name="10km商圏ライン")+
  mapview(market_5km,lwd=3,alpha.regions=0.0,layer.name="5km商圏ライン")

f:id:WAFkw:20180724235824p:plain

よさげ。

この各円(ポリゴン)の中に含まれる小地域は「商圏内」としたい。

円商圏に含まれる小地域を抽出する

sfパッケージにはsf::st_intersects()なる関数とsf::st_intersection()なる似たような名前の関数がある。
いったんその違いについてみてみる。

sf::st_intersection()。戻り値のクラスはsfcなのでそのまま描画できる。

mapview(sf::st_intersection(market_5km,base_map))+
  mapview(market_5km,lwd=3,alpha.regions=0.0,layer.name="5km商圏ライン")

f:id:WAFkw:20180725000540p:plain

なるほど円商圏の端の部分の小地域は"切り取られる"。
これはこれでとっても使い勝手がありそう。
でも今回は人口を集計しちゃいたいので、円商圏が少しでもかぶっている地域は商圏内としてしまいたい。

次にsf::st_intersects()を使ってみる。
これは戻り値がリストなので、
「かぶっているポリゴンのリスト」を使ってベースマップから対象ポリゴンを抽出するのがよさそう...

#intersectionでやると端が切れちゃうのでintersectsで取り直し
market_area_5km <- base_map[sf::st_intersects(market_5km,base_map)[[1]],]
market_area_10km <- base_map[sf::st_intersects(market_10km,base_map)[[1]],]

mapview(base_map,col.regions="gray")+
  mapview(store_point,color="red")+
  mapview(market_area_10km,alpha.regions=0.5,col.regions="yellow")+
  mapview(market_10km,lwd=3,alpha.regions=0.0,layer.name="10km商圏ライン")+
  mapview(market_area_5km,alpha.regions=0.5,col.regions="red")+
  mapview(market_5km,lwd=3,alpha.regions=0.0,layer.name="5km商圏ライン")

f:id:WAFkw:20180725001240p:plain

アップ

f:id:WAFkw:20180725001359p:plain


わーい、きれいに取れた

多すぎるレイヤーを整理しつつ商圏人口を反映

あとは商圏人口を計算して、見た目整えていく。
最近はglueパッケージでラベル作っちゃうことが多い...

#商圏ラインはまとめてしまう
market_line <- sf::st_sf(c(market_5km,market_10km))

# label
label_5km <- glue::glue("5km商圏人口 : {JINKO_5km} 人",JINKO_5km = market_area_5km$JINKO %>% sum() %>% as.character())
label_10km <-  glue::glue("10km商圏人口 : {JINKO_10km} 人",JINKO_10km = market_area_10km$JINKO %>% sum() %>% as.character())

# view
mapview(base_map,col.regions="gray",layer.name="船橋小地域")+
  mapview(store_point,color="black",col.regions="black",layer.name="店舗A")+
  mapview(market_line,lwd=4,alpha.regions=0.0,layer.name="円商圏ライン",label=NA)+
  mapview(market_area_10km,alpha.regions=0.5,col.regions="yellow",layer.name="10km商圏",label=label_10km)+
  mapview(market_area_5km,alpha.regions=0.7,col.regions="red",layer.name="5km商圏",label=label_5km)

f:id:WAFkw:20180725001815p:plain

できた。

フォーカスすると各商圏の人口の総計が見られる。

f:id:WAFkw:20180725002350p:plain

そしてmapviewの何が素敵って、描画後にレイヤーを任意に切り替えできること。
(↓タイルを切り替えてレイヤーを10km商圏のみに)

f:id:WAFkw:20180725002627p:plain

レイヤーやハイライトの切り替えこそが可視化の醍醐味、だって視界の深度は可変だもの。
焦点を合わせるものは少なければ少ない方が美しいのは日常もそうでしょう、といった気持ち。

まとめ

円の取り方、もっとうまく取りたいし最終的にはポリゴンじゃなくてlinestringみたいにしたい...
閃かなかったので詳しいかたどうか。。
あとsfはパイプを使えるのが便利ポイントなはずなのにあまりパイプらなかった気がする。
st_intersectsした後の抽出とかもっと綺麗なコードにならないかな....あとで考える。

追伸

1億年ぶりに書いたから文体も分からなければラーメン屋さんギャグも思い浮かばず。無念でならない。