【sf,mapview】ラーメン屋さんになるにはまだ商圏人口についての理解が足りない
そういえば1年ぐらい前に試される大地を離れて関東に引っ越しました。
文明の利器「えあこん」は素晴らしい。もうあなた無しでは生きられない。
前段
`sf`パッケージなるものがあるらしい、という話を聞いたのはいつのことだったか。
詳しく書いてくださるみなさんの記事を見ているとなるほど気になる....
テーマ
思い返せば私はラーメン屋さんになりたかった気がする。
日々のドタバタにかまけてそんなことも忘れていた。
しかしながら明らかに過去の検討は甘い。
人口で重みづけした重心を求めるだけで立地を決めてよいのか。
もっと集客可能な商圏内についての考察が必要なのではなかろうか。
真剣みが足りない。
今日のゴール
- 店舗の付近円商圏(5km,10km)の地域を特定し、各商圏人口を概算する
- sfに慣れる
- mapviewに慣れる
データの読み込み
いつもは`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 千葉県 船橋市 宮本１丁目 <NA> 8101 95071.29 1386.969 2856 2 12204001002 12 204 001002 千葉県 船橋市 宮本２丁目 <NA> 8101 98592.73 1453.452 3810 3 12204001003 12 204 001003 千葉県 船橋市 宮本３丁目 <NA> 8101 88038.99 1357.873 1045 4 12204001004 12 204 001004 千葉県 船橋市 宮本４丁目 <NA> 8101 103077.04 1623.495 3336 5 12204001005 12 204 001005 千葉県 船橋市 宮本５丁目 <NA> 8101 91841.85 1284.460 2095 6 12204001006 12 204 001006 千葉県 船橋市 宮本６丁目 <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> 宮本１丁目 34 1849 1057 139.9920 35.69872 0010-01 2 <NA> <NA> <NA> <NA> 宮本２丁目 23 1717 963 139.9906 35.69405 0010-02 3 <NA> <NA> <NA> <NA> 宮本３丁目 16 1365 738 139.9934 35.69077 0010-03 4 <NA> <NA> <NA> <NA> 宮本４丁目 21 1216 632 139.9936 35.69238 0010-04 5 <NA> <NA> <NA> <NA> 宮本５丁目 13 805 388 139.9936 35.69508 0010-05 6 <NA> <NA> <NA> <NA> 宮本６丁目 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(.))
何このロンギヌスの槍形状。かっこいい。
とはいえこの鋭利な刃物は海。我がラーメン屋の商圏としてはいかがか、という感じなので
要素指定で削り、ベースマップとして整理する。
base_map <- raw[-337,] base_map %>% mapview::mapview(label=row.names(.))
いい感じ。
店舗位置を中心とした円商圏の設定
まずラーメン屋の座標を設定する。
ここをキャンプ地とする、と私が言ったらそこがラーメン屋なのだ。
船橋駅の座標とぴったり一致してる？...そんなことは知らない。
store_point <- st_point(c(139.985338,35.701557)) mapview::mapview(base_map,label=row.names(base_map))+ mapview::mapview(store_point,color="red")
建設完了。
あとはここを中心に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商圏ライン")
よさげ。
この各円(ポリゴン)の中に含まれる小地域は「商圏内」としたい。
円商圏に含まれる小地域を抽出する
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商圏ライン")
なるほど円商圏の端の部分の小地域は"切り取られる"。
これはこれでとっても使い勝手がありそう。
でも今回は人口を集計しちゃいたいので、円商圏が少しでもかぶっている地域は商圏内としてしまいたい。
次に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商圏ライン")
アップ
わーい、きれいに取れた
多すぎるレイヤーを整理しつつ商圏人口を反映
あとは商圏人口を計算して、見た目整えていく。
最近は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)
できた。
フォーカスすると各商圏の人口の総計が見られる。
そしてmapviewの何が素敵って、描画後にレイヤーを任意に切り替えできること。
(↓タイルを切り替えてレイヤーを10km商圏のみに)
レイヤーやハイライトの切り替えこそが可視化の醍醐味、だって視界の深度は可変だもの。
焦点を合わせるものは少なければ少ない方が美しいのは日常もそうでしょう、といった気持ち。
まとめ
円の取り方、もっとうまく取りたいし最終的にはポリゴンじゃなくてlinestringみたいにしたい...
閃かなかったので詳しいかたどうか。。
あとsfはパイプを使えるのが便利ポイントなはずなのにあまりパイプらなかった気がする。
st_intersectsした後の抽出とかもっと綺麗なコードにならないかな....あとで考える。
追伸
1億年ぶりに書いたから文体も分からなければラーメン屋さんギャグも思い浮かばず。無念でならない。
【leaflet,jsonlite】家を探しのススメ～路線が分からない編～
あけましておめでとうございます。
新年早々水道を止められかけるなどしましたが元気です。
今回のテーマ
土地勘が無いところの家探しをすると仮定します。
世の中には様々な不動産サイトがありますが、
土地勘が無いともうなんだか全部どうでもよくなってきます。
そもそも住みたいエリアも分からんそもそも知らん。となってきます。
せめてエリアを決めるために以下を仮定して、なんとなく便利げなエリアを把握しよう、というのが
今回の目標です。
- 目的地は市ヶ谷 （特に意味はない）
- 乗り換えたくない
- 駅前に住みたい （300M以内)
ではやってみる。
パッケージ
library(leaflet) library(dplyr) library(magrittr) library(jsonlite) library(knitr)
APIを叩いて路線情報を取得
#総武線 11313 #東京メトロ有楽町線 28006 #都営新宿線 99304 #東京メトロ南北線 28009 url<-c("http://www.ekidata.jp/api/l/11313.json", "http://www.ekidata.jp/api/l/28006.json", "http://www.ekidata.jp/api/l/99304.json", "http://www.ekidata.jp/api/l/28009.json") for(i in 1:length(url)){ readLines(url[i],encoding = "UTF-8")[3] %>% gsub("xml.data = ","",.,fixed=TRUE) %>% fromJSON()->JSON switch(i, "1"= JSON->soubu_line, "2"= JSON->yurakutyo_line, "3"= JSON->sinzyuku_line, "4"= JSON->nanboku_line) }
4路線一気にいきたかったけどそうもいかなかった...
いったん路線別に取得したものは、fromJSON()でこんな感じのリストになっています。
> soubu_line $line_cd [1] 11313 $line_name [1] "JR中央・総武線" $line_lon [1] 139.8372 $line_lat [1] 35.70164 $line_zoom [1] 10 $station_l station_cd station_g_cd station_name lon 1 1131301 1131105 三鷹 139.5603 2 1131302 1131104 吉祥寺 139.5798 3 1131303 1131218 西荻窪 139.5994 4 1131304 1131217 荻窪 139.6201 5 1131305 1131216 阿佐ケ谷 139.6359 6 1131306 1131215 高円寺 139.6497 7 1131307 1131214 中野 139.6658 [省略]
「station_l」に入っているdata.frameが目的の路線情報なので
ここを抽出して4路線分くっつけておきます。
前処理
station_lには路線自体の情報が含まれていないので、
後から使うためにstation_cdの前5桁をline_cdとして切り出し、それをもとに路線名称をくっつけます。
これで後からFilter出来るはず。
#市ヶ谷駅 ichigaya <- c("139.736642","35.691295") #4路線を統合、路線CDが失われるので、station_cdの前5桁を取得 line<-dplyr::bind_rows(soubu_line$station_l,yurakutyo_line$station_l,sinzyuku_line$station_l,nanboku_line$station_l) %>% dplyr::mutate(line_cd=as.character(substr(station_cd,1,5))) #路線名マスタ line_cd_master<-data.frame(line_cd=c("11313","28006","99304","28009"), line_name=c("総武線","有楽町線","新宿線","南北線"),stringsAsFactors=FALSE) #路線名をJOIN line<-dplyr::left_join(line,line_cd_master,by="line_cd") line$line_name<-as.factor(line$line_name) #総武線 11313 #東京メトロ有楽町線 28006 #都営新宿線 99304 #東京メトロ南北線 28009 #色パレット関数作成 colpal<-colorFactor(c("tomato","RoyalBlue","MediumSeaGreen","Firebrick"),domain=c("総武線","有楽町線","新宿線","南北線"))
完成後はこんな感じ。
> kable(head(line),format="markdown")
| station_cd | station_g_cd | station_name | lon | lat | line_cd | line_name |
| ----------: | ------------: | :------------ | --------: | --------: | :------- | :--------- |
| 1131301 | 1131105 | 三鷹 | 139.5603 | 35.70268 | 11313 | 総武線 |
| 1131302 | 1131104 | 吉祥寺 | 139.5798 | 35.70312 | 11313 | 総武線 |
| 1131303 | 1131218 | 西荻窪 | 139.5994 | 35.70384 | 11313 | 総武線 |
| 1131304 | 1131217 | 荻窪 | 139.6201 | 35.70452 | 11313 | 総武線 |
| 1131305 | 1131216 | 阿佐ケ谷 | 139.6359 | 35.70482 | 11313 | 総武線 |
| 1131306 | 1131215 | 高円寺 | 139.6497 | 35.70533 | 11313 | 総武線 |
描いちゃう
leaflet(line) %>% addTiles() %>% setView(lng=ichigaya[1],lat=ichigaya[2],zoom=12) %>% addMarkers(lng=ichigaya[1],lat=ichigaya[2]) %>% addPolylines(lng=~lon,lat=~lat,color=~colpal(line_name),weight=3,data=dplyr::filter(.data=line,line_name=="総武線"),group="総武線") %>% addCircles(lng=~lon,lat=~lat,radius=300,color=~colpal(line_name),weight=5,data=dplyr::filter(.data=line,line_name=="総武線"),group="総武線") %>% addPolylines(lng=~lon,lat=~lat,color=~colpal(line_name),weight=3,data=dplyr::filter(.data=line,line_name=="有楽町線"),group="有楽町線") %>% addCircles(lng=~lon,lat=~lat,radius=300,color=~colpal(line_name),weight=5,data=dplyr::filter(.data=line,line_name=="有楽町線"),group="有楽町線") %>% addPolylines(lng=~lon,lat=~lat,color=~colpal(line_name),weight=3,data=dplyr::filter(.data=line,line_name=="新宿線"),group="新宿線") %>% addCircles(lng=~lon,lat=~lat,radius=300,color=~colpal(line_name),weight=5,data=dplyr::filter(.data=line,line_name=="新宿線"),group="新宿線") %>% addPolylines(lng=~lon,lat=~lat,color=~colpal(line_name),weight=3,data=dplyr::filter(.data=line,line_name=="南北線"),group="南北線") %>% addCircles(lng=~lon,lat=~lat,radius=300,color=~colpal(line_name),weight=5,data=dplyr::filter(.data=line,line_name=="南北線"),group="南北線") %>% addLegend(position = "topright", pal = colpal, value = ~line_name, title = "Line") %>% addLayersControl(overlayGroups = c("総武線","有楽町線","新宿線","南北線"),options = layersControlOptions(collapsed = FALSE))
事前に、colorFactor()で各路線をdomainとして定義しておき、
leaflet内では路線ごとにFilterして、それぞれ線を引いたり、駅に円を置いたりする。
addCircles()内のradius引数はm単位なので、今回の「300m」指定に合わせる。
これで駅から300m範囲に円が作れる。
後はaddLayersControl()で各路線を表示ON/OFF出来るようにすればOK。
総武線と有楽町線～とか
チェックをポチポチすればいいので便利。便利だ！！！
おわり
あとはこれに各駅の乗降者数とかを乗っければ混む駅とかが分かる。はず。
平均家賃はスクレイピング要りそうでめんどくさいからやめよう。
後はもうすこし情報載せて駅候補決めて、
不動産へ行って「ぼくこの円の中に住みたいんだよね！」って言えばOK。
.......なんで目的地市ヶ谷にしたんだろう.....謎....
参考
addLayersContoroll()を使ってみたくて書いた記事ですが
冷静に考えたらすでにぞうさん先生が大分以前に台風プロットで使ってらっしゃいました。
分かりやすい.....
P.S
年末色々ありすぎて数々の不義理を重ねました。
謝罪参りしなきゃいけないことが多すぎる。水道局とか水道局とか。
Japan.Rに行けなかったのを今でも引きずっているので
誰か居酒屋.RでLT聞いてやってください。
【maptools】小さいPolygonを結合するだけじゃ足りない、王国を作りたい
「たまにはシムシティじゃなくてCivilizationがやりたい」
そんな自分のためのただのメモ。
今回のテーマ
1.小さな都市をつくる(うそ)
2.小さな都市をさらに繋げて大きな国にする(うそ)
データ
小さな都市を作るための小さな町を集めます。
例のごとく以下のe-statから、シェープファイルをダウンロード。
政府統計の総合窓口 GL01010101
「地図で見る統計(統計GIS)」→「データダウンロード」→「平成22年度国勢調査(小地域)」
1. 札幌市中央区
2. 札幌市北区
3. 札幌市東区
4. 札幌市白石区
5. 札幌市豊平区
6. 札幌市南区
7. 札幌市西区
8. 札幌市厚別区
9. 札幌市手稲区
10.札幌市清田区
全て世界測地系緯度経度･shape形式です。
データの読み込み
とりあえず10ファイル全てを読み込みます。なんか汚い。
#shapefile pj<-CRS("+proj=longlat +datum=WGS84") Chuo_row<- readShapePoly("data/A002005212010DDSWC01101/h22ka01101",proj4string = pj) kita_row<-readShapePoly("data/A002005212010DDSWC01102/h22ka01102",proj4string = pj) higashi_row<-readShapePoly("data/A002005212010DDSWC01103/h22ka01103",proj4string = pj) shiroishi_row<-readShapePoly("data/A002005212010DDSWC01104/h22ka01104",proj4string = pj) toyohira_row<-readShapePoly("data/A002005212010DDSWC01105/h22ka01105",proj4string = pj) minami_row<-readShapePoly("data/A002005212010DDSWC01106/h22ka01106",proj4string = pj) nishi_row<-readShapePoly("data/A002005212010DDSWC01107/h22ka01107",proj4string = pj) atsubetsu_row<-readShapePoly("data/A002005212010DDSWC01108/h22ka01108",proj4string = pj) teine_row<-readShapePoly("data/A002005212010DDSWC01109/h22ka01109",proj4string = pj) kiyota_row<-readShapePoly("data/A002005212010DDSWC01110/h22ka01110",proj4string = pj)
元のシェープファイルには各小地域の人口等、属性テーブルが含まれるのでSpatialPolygonsDataFrameとして読み込まれます。
10個のオブジェクトは完全にバラバラなのでこんな感じ。
小地域データなので各区の中の番地ごとにPolygonが分かれ、ID付与されています。
小さな都市を作る
前も同じようなことをしたので割愛。
wafdata.hatenablog.com
何のことはない、これを
↓
こうするだけ。
maptools::gpclibPermit() Chuo_u<-maptools::unionSpatialPolygons(Chuo_row,IDs=rep(1,times=nrow(Chuo@data))) kita_u<-maptools::unionSpatialPolygons(kita_row,IDs=rep(2,times=nrow(kita_row@data))) higashi_u<-maptools::unionSpatialPolygons(higashi_row,IDs=rep(3,times=nrow(higashi_row@data))) shiroishi_u<-maptools::unionSpatialPolygons(shiroishi_row,IDs=rep(4,times=nrow(shiroishi_row@data))) toyohira_u<-maptools::unionSpatialPolygons(toyohira_row,IDs=rep(5,times=nrow(toyohira_row@data))) minami_u<-maptools::unionSpatialPolygons(minami_row,IDs=rep(6,times=nrow(minami_row@data))) nishi_u<-maptools::unionSpatialPolygons(nishi_row,IDs=rep(7,times=nrow(nishi_row@data))) atsubetsu_u<-maptools::unionSpatialPolygons(atsubetsu_row,IDs=rep(8,times=nrow(atsubetsu_row@data))) teine_u<-maptools::unionSpatialPolygons(teine_row,IDs=rep(9,times=nrow(teine_row@data))) kiyota_u<-maptools::unionSpatialPolygons(kiyota_row,IDs=rep(10,times=nrow(kiyota_row@data)))
王国を作るという野望のため、
結合後のIDは各区ごとにユニークになるように振っておく。
ここで各オブジェクトはSpatialPolygonsになる。
王国をつくる
いざバラバラのオブジェクトを一つのSpatialPolygonに。
これを
↓
こうする。
これをやるためには、Spatialオブジェクトのrbind、ことmaptools::spRbind()を使う。
今回はSpatialPolygonsとSpatialPolygonsの合わせ技だけど、関数としては以下の組み合わせに対応出来るらしい。
| 引数1 | 引数2 |
|---|---|
| SpatialPoints | SpatialPoints |
| SpatialPointsDataFrame | SpatialPointsDataFrame |
| SpatialLines | SpatialLines |
| SpatialLinesDataFrame | SpatialLinesDataFrame |
| SpatialPolygons | SpatialPolygons |
| SpatialPolygonsDataFrame | SpatialPolygonsDataFrame |
(ちなみに*DataFlame系の場合は同じ属性テーブルを有する場合に限る。ちょっと試せてはいない。)
でも....うっすら気付いてたけどこれって.....と思って関数ヘルプを眺めていたら
spRbind-methods can bind 2 objects, whereas rbind-methods can bind multiple object
10個一気に繋げたいんだけど....なあ.....
ごにゃごにゃやってみたけど結局以下で愚直に王国を作った。
Sapporo<-spRbind(Chuo_u,kita_u) %>% spRbind(.,higashi_u) %>% spRbind(.,shiroishi_u) %>% spRbind(.,toyohira_u) %>% spRbind(.,minami_u) %>% spRbind(.,nishi_u) %>% spRbind(.,atsubetsu_u) %>% spRbind(.,teine_u) %>% spRbind(.,kiyota_u)
最終的に、1つのSpatialPolygons(内部に区ごとの10個のPolygonを有する)が出来た。
str(Sapporo) >Formal class 'SpatialPolygons' [package "sp"] with 4 slots ..@ polygons :List of 10 .. ..$ :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 .. ..$ :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.4 43.1 .. .. .. .. .. .. ..@ area : num 0.00701 .. .. .. .. .. .. ..@ hole : logi FALSE .. .. .. .. .. .. ..@ ringDir: int 1 .. .. .. .. .. .. ..@ coords : num [1:1493, 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.4 43.1 .. .. .. ..@ ID : chr "2" .. .. .. ..@ area : num 0.00701 .. ..$ :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.4 43.1 .. .. .. .. .. .. ..@ area : num 0.0063 .. .. .. .. .. .. ..@ hole : logi FALSE .. .. .. .. .. .. ..@ ringDir: int 1 .. .. .. .. .. .. ..@ coords : num [1:744, 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.4 43.1 .. .. .. ..@ ID : chr "3" .. .. .. ..@ area : num 0.0063 #･･･長いので省略
これでどんどん領地を拡大していけばかのアッバース朝みたいな国を作れる。
【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 "旭ケ丘１丁目",..: 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 | 北十六条西１６丁目 | 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月なのかな。お待ちしてます。
【sp,deldir,leaflet】あなたの全てを知りたい編～Polygon内に含まれるポイントデータを集計～
GWなので家出をしました。
帰りたくないけどおうちに帰るまでが家出だからなあ...
本日のテーマ
掲題の通りです。
あなた(Polygon)の全部が知りたいのよ！という気持ちなので
Polygon内に含まれる複数の点の情報を集計するなどしたいです。
たとえばさ、
算出したポリゴン内に複数の世帯が含まれていたとして、
そのポリゴン内の人口を算出したい、みたいなの。
データ
そうはいっても不条理な世の中では私のような下々のパンピーが世帯別のデータを手に入れることは出来ません。
乱数作ってポリゴン内に配置して～とかが主流なのでしょうが、
妄想でもいいのでもう少しリアルに近いことがしたい。
ということで以下の通りに。
#------------------------
①札幌市中央区の小地域データからポリゴン内代表点を抽出。
→「世帯」だと思い込んで取扱い
②札幌市中央区内の小学校の座標を用いて中央区をボロノイ分割
→集計したい「ポリゴン」を作成
③②に含まれる①の点を集計してポリゴン内人口を算出する。
#------------------------
イメージは「小学校区(ボロノイ領域基準)内の人口を把握する」という分析をする感じです。
ざっくり同じくらいの人口になるんじゃないかなと仮説を立てつつ。
妄想の世界へ
①、②は過去にやったやつの焼き直し的なコードなので細かい記述は省きます。
wafdata.hatenablog.com
wafdata.hatenablog.com
#データの準備 pj<-CRS("+proj=longlat +datum=WGS84") #緯度経度なので+projにlonglatを指定 Chuo<- maptools::readShapePoly(shape.file,proj4string = pj) #SpatialPolygonDataframeとして読み込み #lng,latに代表点を格納(世帯のダミーデータ) Chuo.d<- 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, lng = coordinates(Chuo)[,1], lat=coordinates(Chuo)[,2]) knitr::kable(head(Chuo.d,10),format="markdown") | |MOJI | X_CODE| Y_CODE|KEY_CODE | JINKO| SETAI| lng| lat| |:--|:------------------|--------:|--------:|:-----------|-----:|-----:|--------:|--------:| |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 |北十六条西１６丁目 | 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| |6 |北十六条西１５丁目 | 141.3304| 43.07795|01101791615 | 142| 80| 141.3304| 43.07795| |7 |北十五条西１５丁目 | 141.3306| 43.07682|01101791515 | 285| 197| 141.3306| 43.07682| |8 |北十六条西２１丁目 | 141.3209| 43.07592|01101791621 | 41| 20| 141.3209| 43.07592| |9 |北十四条西１５丁目 | 141.3314| 43.07522|01101791415 | 699| 307| 141.3314| 43.07522| #描画してみる ##カラーパレットの生成(連続値) pal <- colorNumeric( palette = "YlOrRd", domain = Chuo@data$JINKO ) ##人口コロプレスに代表点を載せて描画 (m<-leaflet::leaflet(Chuo) %>% addProviderTiles(provider="OpenStreetMap.BlackAndWhite") %>% addPolygons(stroke = FALSE, fillOpacity = 0.9, color = ~pal(Chuo.d$JINKO)) %>% addCircles(lng=coordinates(Chuo)[,1],lat=coordinates(Chuo)[,2],color = "gray80"))
この黒いツブツブが代表点で、各ポリゴン一つずつ。
元データに小地域ごとの人口が入っているので、そこと紐づけて世帯＋世帯人口として取り扱います。
次にボロノイ分割で領域作り。
小学校の座標データがいまいち取れなかったので作りました。
#ボロノイ #データの読み込み school<-read.csv("school0501.csv",header = TRUE) #16校くらい。学校名＋lng+lat #描画の範囲を固める xrng <- scales::expand_range(range(school$lng), .05) yrng <- scales::expand_range(range(school$lat), .05) #ボロノイ分割(戻り値はlist) result_deldir <- deldir::deldir(x=school$lng,y=school$lat, rw = c(xrng, yrng)) tilelist <- deldir::tile.list(result_deldir) #分割したtileを取得 #内部構造を見てみる str(tilelist,10) List of 16 $ :List of 5 ..$ ptNum: int 1 ..$ pt : Named num [1:2] 141.4 43.1 .. ..- attr(*, "names")= chr [1:2] "x" "y" ..$ x : num [1:4] 141 141 141 141 ..$ y : num [1:4] 43.1 43.1 43 43 ..$ bp : logi [1:4] TRUE TRUE FALSE TRUE $ :List of 5 ..$ ptNum: int 2 ..$ pt : Named num [1:2] 141.3 43.1 .. ..- attr(*, "names")= chr [1:2] "x" "y" ..$ x : num [1:6] 141 141 141 141 141 ... ..$ y : num [1:6] 43.1 43.1 43.1 43 43 ... ..$ bp : logi [1:6] TRUE TRUE FALSE FALSE FALSE FALSE #(長いので省略) #ボロノイ分割の結果を見てみる ggplot2::qplot(lng, lat, data = school)+geom_segment(aes(x = x1, y = y1, xend = x2, yend = y2), size = .25,data = result_deldir$dirsgs)
#ボロノイ領域のポリゴン化 polys<-vector(mode='list', length=length(tilelist)) for (i in seq(along=polys)) { voronoi<-cbind(tilelist[[i]]$x, tilelist[[i]]$y) voronoi<-rbind(voronoi, voronoi[1,]) polys[[i]]<-Polygons(list(Polygon(voronoi)), ID=as.character(i)) } SP<-SpatialPolygons(polys) slotNames(SP) #viewの準備 sapporo <- leaflet() %>% setView(lng = 141.3508, lat = 43.06862,zoom = 10) (school_voronoi<-sapporo %>% addProviderTiles("OpenStreetMap.Mapnik") %>% addPolygons(stroke=FALSE,data = SP) %>% addPolylines(data = SP,color = "black",weight = 2.5)) #小地域代表点を世帯と思い込んで school_voronoi %>% addCircles(data=Chuo.d,lng=coordinates(Chuo)[,1],lat=coordinates(Chuo)[,2],color = "gray80")
ggplot2で描いたやつのポリゴン内点は各小学校(母点)、作成されたポリゴンがボロノイ領域です。
leafletの上図は結果をSpatialPolygonにして格納・描画したもの。まだ色はつけてません。
下図はさっき用意した世帯をその上に乗せたものです。
鳥肌が立ちますね。
これで戦の準備が整いました。
本題～ポリゴン内集計
どうやらsp::over関数を使えばいけそうである。
引数にfunctionを指定出来るようなので、sumを指定してJINKO(人口)列を集計する。
また、over関数は使用前に代表点を(lng,lat)の形で格納しておく必要がありそう。
coordinates(df)<-~x+yの形でやったらうまくできた。
#代表点を格納 coordinates(Chuo.d)<- ~lng+lat #集計値算出 result_table<-as.data.frame(over(SP,Chuo.d[,5:6],fn=sum)) >knitr::kable(result_table) | JINKO| SETAI| |-----:|-----:| | 15917| 10583| | 18987| 14341| | 19397| 12922| | 23965| 14162| | 14638| 7978| | 14459| 7186| | 12187| 6455| | 22489| 12870| | 14146| 6102| | 9593| 5084| | 15983| 6439| | 6866| 3036| | 11778| 5494| | 7905| 2553| | 2255| 926| | NA| NA|
ポリゴンが16あるので、各ポリゴンごとに16セットできた。
ポリゴン内に世帯が一つもないところはNAになっている。
この集計結果をさきほどのボロノイポリゴンに色付けして完成にする。
#結果を@dataに格納し、SPDFに変換 result_polygon<-SpatialPolygonsDataFrame(SP,result_table) #NAを0に置換 result_polygon@data$JINKO[is.na(result_polygon@data$JINKO)]<-0 #描画 pal2<-colorNumeric("Oranges", domain = range(result_polygon@data$JINKO)) (result_voronoi<-sapporo %>% addProviderTiles("OpenStreetMap.BlackAndWhite") %>% addPolygons(stroke=FALSE, data = result_polygon, color = ~pal2(result_polygon@data$JINKO), fillOpacity = 0.9) %>% addPolylines(data = result_polygon,color = "black",weight = 2.5))
できた。満足。
感想
山間(左側)は人口が他に比べて少なめ。学校のボロノイ領域が広いのも元々人口少なめなのが関係してるのかも。
やはり中心部は人口多いですねー小学校も生徒人数多かったりするのかしら。
詳しく調べるには世帯の年齢とかも加味しないとダメそうなのでこれ以上のコメントはやめておく.....
おわり
先日久々のTokyoRへ行ってきました。次はいつ来られるかなあ
ちなみにSappoRoRは7月9日です。ぜひ。
本州は暑くてね.......
【maptools,leaflet,rgdal】 (愛と希望で)世界を一つにまとめたい
何の予定も無いけど、なんだか突然LTがしたくなったので
自宅.Rをやるべく資料作りに勤しむ遊びにハマっています。
待てど暮らせど北国に春が来ない。
データセット
平成22年 国勢調査(小地域)境界データ
(地図で見る統計（統計ＧＩＳ） 政府統計の総合窓口 GJ01000001)
「データダウンロード」より進み、
「札幌市中央区-男女別人口総数及び世帯総数」から「世界測地系・緯度経度.shape形式」を取得します。
データの読み込み(世界測地系･緯度経度)
早速データを読み込み。今回は世界測地系(緯度経度)なんですが、e-statでは平面直角座標のshape形式も提供されているので
そっちの読み込み方法も少しメモ。
ライブラリ
library(maptools) library(leaflet) library(rgdal) library(magrittr)
読み込み
shape.file #ダウンロードした.shapeのフルパス格納 > getinfo.shape(shape.file) Shapefile type: Polygon, (5), # of Shapes: 990 pj<-CRS("+proj=longlat +datum=WGS84") #緯度経度なので+projにlonglatを指定 Chuo<- maptools::readShapePoly(shape.file,proj4string = pj) #SpatialPolygonDataframeとして読み込み #確認 > 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 "旭ケ丘１丁目",..: 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" ... #変数が多すぎなので、e-stat上の定義書を見ながら使いそうなものを確認 Chuo.d<- 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) |MOJI | X_CODE| Y_CODE|KEY_CODE | JINKO| SETAI| |:------------------|--------:|--------:|:-----------|-----:|-----:| |北二十二条西 | 141.3284| 43.08521|011017922 | 29| 1| |北二十一条西 | 141.3285| 43.08400|011017921 | 518| 304| |北二十条西 | 141.3280| 43.08234|011017920 | 206| 108| |北十六条西１６丁目 | 141.3260| 43.07778|01101791616 | 0| 0| |北十八条西 | 141.3298| 43.08034|011017918 | 136| 72| |北十七条西 | 141.3302| 43.07896|011017917 | 149| 72| |北十六条西１５丁目 | 141.3304| 43.07795|01101791615 | 142| 80| |北十五条西１５丁目 | 141.3306| 43.07682|01101791515 | 285| 197| |北十六条西２１丁目 | 141.3209| 43.07592|01101791621 | 41| 20| |北十四条西１５丁目 | 141.3314| 43.07522|01101791415 | 699| 307|
小地域は市町村よりもさらに細かいデータなので、〇丁目とか、〇番地とかそういう単位で構成されています。
定義書曰く、
X_CODE 小地域の図形中心座標
Y_CODE 小地域の図形中心座標
KEYCODE 丁目・番地ユニークコード(都道府県コード＋市町村コード+小地域)
JINKO 該当地域人口
SETAI 該当地域世帯総数
データの読み込み(平面直角座標系だった場合)
ちなみに、平面直角座標系の場合はCRS()に直接+projを指定出来ません。
なので、
平面直角座標系の系統番号→EPSG(European Petroleum Survey Group)の番号読み替え→+init引数に指定
の手順になります。
(参考:平面直角座標系(平成18年改定)国土地理院)
国勢調査の定義書曰く、「北海道はⅫ系だよ！」とのことなので、EPSGのコードだと「2454」に該当しそう。
>pj<-CRS("+init=epsg:2454") #読み替えたEPSGで指定 >rgdal::CRSargs(pj) #中身を見ると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) #読み込み
まあメモ程度に。
いったん描画
楽なのでleafletで描きます。
ついでなので@dataの人口を使ってコロプレス図にしてみる。
#カラーパレットの生成(連続値) pal <- colorNumeric( palette = "YlOrRd", domain = Chuo@data$JINKO ) #コロプレス図だけ描くver leaflet::leaflet(Chuo) %>% addProviderTiles(provider = "OpenStreetMap.BlackAndWhite") %>% addPolygons(stroke = FALSE, fillOpacity = 0.9, color = ~pal(Chuo.d$JINKO)) #代表点を求めて打つver leaflet::leaflet(Chuo) %>% addProviderTiles(provider = "OpenStreetMap.BlackAndWhite") %>% addPolygons(stroke = FALSE, fillOpacity = 0.9, color = ~pal(Chuo.d$JINKO)) %>% addCircles(lng = coordinates(Chuo)[,1],lat=coordinates(Chuo)[,2],color = "gray80")
下の図鳥肌が立つ。。。。。
ここでちょっとstrokeが気になり始める
leaflet::addPolygons()とかの描画系関数の中にあるstroke引数は、ポリゴンの外枠の線を制御します。
stroke=TRUE 外線あり
stroke=FALSE 外線なし
描いてみるとこんな感じ
leaflet::leaflet(Chuo) %>% addProviderTiles(provider = "OpenStreetMap.BlackAndWhite") %>% addPolygons(stroke = TRUE, fillOpacity = 0.9, color = ~pal(Chuo.d$JINKO)) %>% addPolylines(color="gray80") #外線の色変えのためにaddPolylineを指定
小地域だとなんか線がいっぱいになって見づらい。
でも外線一切なしだと儚げな感じになっちゃってそれもそれで見づらい。
→一番外側(丁目単位ではなく、札幌市中央区の外側)だけ色指定したい！！！！
そこでポンコツはかんがえた
SpatialPolygonは、たくさんのPolygonが集合して全体の形を構成している。ので、
他のPolygonと接していない部分を「一番外側」として定義しない限り、思い描いたようには出来なそう。
→ええい
→もうポリゴンを一つにしちゃったらいいじゃない。
→世界をひとつにしちゃえばいいのよ。
やってみる
#元のPolygon。小地域の数990個分ある。 > getinfo.shape(shape.file) Shapefile type: Polygon, (5), # of Shapes: 990 maptools::gpclibPermit() #こいつを先にやっとかないとなんか怒られる #ポリゴンの結合 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"
maptools::unionSpatialPolygons()でポリゴンの結合が出来るのだけど、
今回は990個をぜーーーんぶ一つにしてしまいたいので、IDsには全部1を指定する。
IDsは「新しいID」をベクトル指定出来るので、たとえば
元のID=c(1,2,3)なPolygonがあったら
IDs=c(1,1,2)みたいに指定すると、任意にくっつけるPolygonを指定出来る。（この場合ID=3だけ結合されない)
てことでちゃんと結合できたか描画してみる
leaflet::leaflet(Chuo.union) %>% addProviderTiles(provider = "OpenStreetMap.BlackAndWhite") %>% addPolygons(stroke = TRUE)
小地域の区切りが消えてすっきり。
重ねて描いてみる
最初に描いてたコロプレスと、さっき結合したものを重ねて外側だけ線を引く。
#コロプレスと重ねる leaflet::leaflet(Chuo) %>% addProviderTiles(provider = "OpenStreetMap.BlackAndWhite") %>% addPolygons(stroke = FALSE, fillOpacity = 0.9, color = ~pal(Chuo.d$JINKO)) %>% addPolygons(data = Chuo.union,stroke = TRUE,fill = FALSE,color="gray80") #代表点と重ねる leaflet::leaflet(Chuo) %>% addProviderTiles(provider = "OpenStreetMap.BlackAndWhite") %>% addCircles(lng = coordinates(Chuo)[,1],lat=coordinates(Chuo)[,2]) %>% addPolygons(data = Chuo.union,stroke = TRUE,fill = FALSE,color="blue")
できた。
閑話休題
世界をつなげる力を手に入れた。
そんなことより、@dataに含まれるX_CODE,Y_CODEに謎がある。
sp::coordinates()で代表点を求めずに、図形中心座標(X_CODE,Y_CODE)で点を打つとなんか微妙にズレる。
leaflet::leaflet(Chuo) %>% addProviderTiles(provider = "OpenStreetMap.BlackAndWhite") %>% addPolygons(stroke = FALSE, fillOpacity = 0.9, color = ~pal(Chuo.d$JINKO)) %>% addCircles(lng = Chuo@data$X_CODE,lat=Chuo@data$Y_CODE)
なんか左上にズレる.....。なんで....。どうしてなのパトラッシュ....
ここで言ってるX_CODEとY_CODEが何をどうやって求めてるのか分かればすっきり出来るのかしら。
おわり
ねむいのでおわります。
足りないものが足りないままなのすごくいらいらするなあ、というまいにち。
【sp,ggplot2,leaflet】日曜日なので豆について考えてみた（3）
待てど暮らせど春が来ない。
本日のテーマ
使い回し続けている大豆データで、
コロプレス図を色んな方法で書きまくります。
コロプレス図、色んな方法で描けるのだけど、
各方法の描画の仕方と見え方について比較などしてみたくなった、というのが動機です。
あとどさくさに紛れてWindowsfontの適用をやっと覚えたのでメモ。
データ
wafdata.hatenablog.com
wafdata.hatenablog.com
シリーズで使いまわしている大豆に関するデータ(農林水産省/大豆関連データ集)と、
Nipponパッケージをインストールすると入る日本のシェープファイルを使います。
前処理については前回(たぶん)やっているので端折りますが、
・NipponパッケージのシェープファイルをSpatialPolygonsDataFrameとして読み込み
・@data部に、大豆データの「都道府県別納豆消費量(円)」「都道府県別豆腐消費量(円)」「都道府県別大豆生産量(t)」をマージ
して、以下の状態になっています。
#データ部だけ抜き出すとこんな感じ > knitr::kable(map@data) |SP_ID |jiscode |name | population|region | natto| tofu| mame_t| |:-----|:-------|:---------|----------:|:----------------|-----:|----:|------:| |1 |01 |Hokkaido | 5506419|Hokkaido | 3847| 4362| 73600| |2 |02 |Aomori | 1373339|Tohoku | 4807| 5140| 5370| |3 |03 |Iwate | 1330147|Tohoku | 4654| 6945| 5470| |4 |04 |Miyagi | 2348165|Tohoku | 5485| 6264| 19300| |5 |05 |Akita | 1085997|Tohoku | 4477| 5123| 9640| |6 |06 |Yamagata | 1168924|Tohoku | 5569| 5906| 7720| |7 |07 |Fukushima | 2029064|Tohoku | 5298| 5254| 2220| |8 |08 |Ibaraki | 2969770|Kanto | 5916| 4839| 5410| |9 |09 |Tochigi | 2007683|Kanto | 3962| 5111| 4250| |10 |10 |Gunma | 2008068|Kanto | 4808| 5226| 444| |11 |11 |Saitama | 7194556|Kanto | 4406| 6328| 786| |12 |12 |Chiba | 6216289|Kanto | 4064| 6151| 1080| |13 |13 |Tokyo | 13159388|Kanto | 3618| 5969| 4| |14 |14 |Kanagawa | 9048331|Kanto | 3418| 5820| 64| |15 |15 |Niigata | 2374450|Chubu | 4586| 5070| 8840| |16 |16 |Toyama | 1093247|Chubu | 4329| 4927| 7630| |17 |17 |Ishikawa | 1169788|Chubu | 3505| 5176| 2310| |18 |18 |Fukui | 806314|Chubu | 3238| 4815| 2430| |19 |19 |Yamanashi | 863075|Chubu | 4501| 5274| 289| |20 |20 |Nagano | 2152449|Chubu | 4157| 4510| 3380| |21 |21 |Gifu | 2080773|Chubu | 3160| 5078| 3430| |22 |22 |Shizuoka | 3765007|Chubu | 3735| 6057| 317| |23 |23 |Aichi | 7410719|Chubu | 3348| 5350| 7180| |24 |24 |Mie | 1854724|Chubu | 2429| 5307| 3750| |25 |25 |Shiga | 1410777|Kinki | 2912| 5059| 9450| |26 |26 |Kyoto | 2636092|Kinki | 2680| 5982| 455| |27 |27 |Osaka | 8865245|Kinki | 1945| 5470| 20| |28 |28 |Hyogo | 5588133|Kinki | 2289| 5709| 3320| |29 |29 |Nara | 1400728|Kinki | 2421| 5408| 256| |30 |30 |Wakayama | 1002198|Kinki | 1994| 4864| 43| |31 |31 |Tottori | 588667|Chugoku | 2658| 6728| 1120| |32 |32 |Shimane | 717397|Chugoku | 2605| 6213| 1320| |33 |33 |Okayama | 1945276|Chugoku | 2428| 4534| 2370| |34 |34 |Hiroshima | 2860750|Chugoku | 2592| 5258| 667| |35 |35 |Yamaguchi | 1451338|Chugoku | 2744| 4971| 1020| |36 |36 |Tokushima | 785491|Shikoku | 2306| 6763| 45| |37 |37 |Kagawa | 995842|Shikoku | 2528| 5425| 97| |38 |38 |Ehime | 1431493|Shikoku | 2628| 6170| 457| |39 |39 |Kochi | 764456|Shikoku | 2017| 6009| 61| |40 |40 |Fukuoka | 5071968|Kyushu / Okinawa | 3125| 4563| 14300| |41 |41 |Saga | 849788|Kyushu / Okinawa | 3317| 6053| 15300| |42 |42 |Nagasaki | 1426779|Kyushu / Okinawa | 2672| 4824| 496| |43 |43 |Kumamoto | 1817426|Kyushu / Okinawa | 3861| 5417| 3710| |44 |44 |Oita | 1196529|Kyushu / Okinawa | 3604| 6288| 1690| |45 |45 |Miyazaki | 1135233|Kyushu / Okinawa | 2937| 5239| 317| |46 |46 |Kagoshima | 1706242|Kyushu / Okinawa | 3513| 5471| 301| |47 |47 |Okinawa | 1392818|Kyushu / Okinawa | 2389| 7478| 1|
早速描く～spplot編～
spパッケージの中に、spplot()なる関数があるので、これでSP系の地図描画が出来ます。
これの便利なとこは、
何も考えずにSPDFを突っ込んじゃってOKなところ、と、
zcolという引数にc()で束ねた列を与えることで、同じ凡例を使った複数のコロプレス図を並べて描けるところ。
(何を言っているのか分からない文章)
とりあえずまずはデフォルトの状態で描いてみる。
#納豆と豆腐の消費額コロプレス図を並べて描く > sp::spplot(map,zcol=c("natto","tofu"),names.attr=c("納豆","豆腐"))
簡単にすっきり並べられた！！
....けどこう...なんかこう...デフォルトの色味がこう...いかがわしい。
なんかとてもいかがわしい色味。
ということでもう少し手を入れてみる。
消費額を離散化して、at引数に入れる。（1000円ごとに8グループ）
かつ、RColorBrewerパッケージの色パレットを適用させてみる。
#8階層に分けて色を決める c(0,1000,2000,3000,4000,5000,6000,7000,8000) %>% sp::spplot(map,zcol=c("natto","tofu"),at=.,names.attr=c("納豆","豆腐"), col.regions=RColorBrewer::brewer.pal(8,"BuPu"))
清純な感じになりました。
親御さんもこれで一安心です。
結果については前回記事でもグダグダ言ってましたが、
納豆の地域差が東西ではっきり出ていて、豆腐はそんなに地域差は無さそう,,
早速描く～ggplot2編～
次はみんな大好きggplot2で描きます。
が、その前に、
「plot内のフォントがどうしても綺麗じゃなくてやる気が出ない」
症状が現れ始めたので、fontについて考えます。
ggplot2の描画にお好きなfontを指定するには、
theme()を使ってelement_textを変更して、(font)familyを指定してあげればOK！なんだけど、
その前にfontfamilyにお好きなfontを追加したい...メイリオとかメイリオとかメイリオとかね...
font変更周りはmacの情報しか見つからず、
挙げ句うまくいかなくて唸っておりました。
が、下記記事でWindowsfont()をかませないとダメなことを知る。
http://sudori.info/stat/stat_fig_font.htmlsudori.info
#フォントファミリーの指定 #メイリオを"MEI"というfamilyとして定義 windowsFonts(MEI=windowsFont("Meiryo"))
これでいける、はず。
気を取り直して描画。
ggplot2ではSPDFも受け入れてくれるらしいのだけれど、
いったんfortify()してデータ整形する練習をしてみる。
> map_df<-ggplot2::fortify(map) #dataframeに変換 > knitr::kable(head(map_df,10)) | long| lat| order|hole |piece |group |id | |--------:|-------:|-----:|:-----|:-----|:-----|:--| | 139.7707| 42.3018| 1|FALSE |1 |0.1 |0 | | 139.8711| 42.6623| 2|FALSE |1 |0.1 |0 | | 140.1895| 42.8243| 3|FALSE |1 |0.1 |0 | | 140.3062| 42.7741| 4|FALSE |1 |0.1 |0 | | 140.5259| 42.9885| 5|FALSE |1 |0.1 |0 | | 140.3279| 43.2196| 6|FALSE |1 |0.1 |0 | | 140.3500| 43.3315| 7|FALSE |1 |0.1 |0 | | 140.4959| 43.3704| 8|FALSE |1 |0.1 |0 | | 140.7945| 43.1928| 9|FALSE |1 |0.1 |0 | | 141.0210| 43.2348| 10|FALSE |1 |0.1 |0 |
SPDFをfortify()すると、
ポリゴン部の情報がdataframe化されるので、元々data部にあったregionや納豆&豆腐情報は載ってこない。
描画に使用するデータをマージする必要がありそう。
(変換前)SPDFのdata部・・・・・SP_ID:1～47
(変換後)fortify()後・・・・・id:0～46
ちょっとごちゃごちゃするけど、今回はSP_IDを1ズラして、idとJOINすることにした。
（後々考えたら行番号で一致させるのが正解。たぶん。)
>map_join<-data.frame(id=as.numeric(as.character(map@data$SP_ID)-1),map@data)#ポリゴンIDとのJOINKeyを強引に生成 >map_join$id<-as.character(map_join$id) >map_df<- dplyr::left_join(map_df,map_join,by="id") #JOIN > knitr::kable(head(map_df,10)) | long| lat| order|hole |piece |group |id |SP_ID |jiscode |name | population|region | natto| tofu| mame_t| |--------:|-------:|-----:|:-----|:-----|:-----|:--|:-----|:-------|:--------|----------:|:--------|-----:|----:|------:| | 139.7707| 42.3018| 1|FALSE |1 |0.1 |0 |1 |01 |Hokkaido | 5506419|Hokkaido | 3847| 4362| 73600| | 139.8711| 42.6623| 2|FALSE |1 |0.1 |0 |1 |01 |Hokkaido | 5506419|Hokkaido | 3847| 4362| 73600| | 140.1895| 42.8243| 3|FALSE |1 |0.1 |0 |1 |01 |Hokkaido | 5506419|Hokkaido | 3847| 4362| 73600| | 140.3062| 42.7741| 4|FALSE |1 |0.1 |0 |1 |01 |Hokkaido | 5506419|Hokkaido | 3847| 4362| 73600| | 140.5259| 42.9885| 5|FALSE |1 |0.1 |0 |1 |01 |Hokkaido | 5506419|Hokkaido | 3847| 4362| 73600| | 140.3279| 43.2196| 6|FALSE |1 |0.1 |0 |1 |01 |Hokkaido | 5506419|Hokkaido | 3847| 4362| 73600| | 140.3500| 43.3315| 7|FALSE |1 |0.1 |0 |1 |01 |Hokkaido | 5506419|Hokkaido | 3847| 4362| 73600| | 140.4959| 43.3704| 8|FALSE |1 |0.1 |0 |1 |01 |Hokkaido | 5506419|Hokkaido | 3847| 4362| 73600| | 140.7945| 43.1928| 9|FALSE |1 |0.1 |0 |1 |01 |Hokkaido | 5506419|Hokkaido | 3847| 4362| 73600| | 141.0210| 43.2348| 10|FALSE |1 |0.1 |0 |1 |01 |Hokkaido | 5506419|Hokkaido | 3847| 4362| 73600|
やり方きれいじゃないですが、とりあえずJOINは成功。
ついでに離散化周りをやってから、描画。
#離散化 cut.line<-c(0,1000,2000,3000,4000,5000,6000,7000,8000) cut.label<-c("0-1000円","1000-2000円","2000-3000円","3000-4000円","4000-5000円","5000-6000円","6000-7000円","7000-8000円") map_df$natto.c<-cut(map_df$natto,cut.line,labels=cut.label) map_df$tofu.c<-cut(map_df$tofu,cut.line,labels=cut.label) #描画 (nattomap_2<-ggplot(map_df,aes(x=long,y=lat,map_id=region,fill=natto.c))+ geom_map(map=map_df,colour="gray",size=0.8)+ scale_fill_manual(values=RColorBrewer::brewer.pal(8,"Reds"))+ ggtitle("都道府県別納豆消費量")+ theme(plot.title=element_text(size = rel(2),family = "MEI"))+ coord_map(projection="mercator")) (tofumap_2<-ggplot(map_df,aes(x=long,y=lat,map_id=region,fill=tofu.c))+ geom_map(map=map_df,colour="gray",size=0.8)+ scale_fill_manual(values=RColorBrewer::brewer.pal(8,"YlOrRd"))+ ggtitle("都道府県別豆腐消費量")+ theme(plot.title=element_text(size = rel(2),family = "MEI"))+ coord_map(projection="mercator")) #並べてみる grid.arrange(nattomap_2,tofumap_2,ncol=2,nrow=1)
描画のところのコードが汚くなってしまった。猛省。。
geom_map()の方のmapの指定は省略したいんだけどなあ..
さっき定義したメイリオフォントはうまく適用できた。幸せ。
しかし、
せっかく同じ基準で離散化してるのに、パレットがなんか5色と4色になってしまっているのが
どうにも納得がいかない。
外出しでちゃんと8レベル分色を決めて、凡例にも8レベル分出したいなあ...む。
spplot()の結果のように、同じ凡例を使えるようにしたいのだけれど。
なにはともあれ
ggplot2で描くコロプレス図の特徴は、
自由度が高いことと綺麗なこと、かな。もっと細かい指定も出来るしちゃんと描く時はこれがメインの方法になりそう。
早速描いてみる～leaflet編～
なんだか眠くなってきたのでさらっと。
SPDFをleafletにそのまま渡すと、
x,yの指定なしの場合はdata内のlong,latを見つけてきてマッピングしてくれます。
leafletでは離散化せず、元の消費額を連続値のままcolorNumeric()を使って色決めしてみます。
map_l<-map #色決め colset.n<-colorNumeric("Reds", domain = range(map_l@data$natto)) colset.t<-colorNumeric("Oranges", domain = range(map_l@data$tofu)) map_l$natto.col<-colset.n(map_l@data$natto) #わざわざ列を作らなくていい気がする map_l$tofu.col<-colset.t(map_l@data$tofu) #わざわざ列を作らなくていい気がする #描画 leaflet() %>% addTiles() %>% addPolygons(data = map_l, color = ~natto.col, fillOpacity = 0.8, stroke = FALSE) leaflet() %>% addTiles() %>% addPolygons(data = map_l, color = ~tofu.col, fillOpacity = 0.8, stroke = FALSE)
※動きません
leafletは描くのも楽だし見た目も綺麗だし
少し透過させるとTileの地図情報も読み取れるのがいいですね。
ちからつきる
じつはcoloplethrパッケージで描くのもやりたかったんだけど、
ほんじつはここでちからつきた
ゆうしゃろと、ここにねむる。
P.S 書籍発売ラッシュ
- 作者: Hadley Wickham,石田基広,市川太祐,高柳慎一,福島真太朗
- 出版社/メーカー: 共立出版
- 発売日: 2016/02/10
- メディア: 単行本
- この商品を含むブログ (28件) を見る
かのHadley神の「AdvancedR」邦訳版がついに。
個人的には7章が好きです。
7章が好きです。
Rをはじめたての頃は、結果が返ればまあいいや、だったところ、
だんだん関数のソースを見てごにょごにょするようになってくると
ここのメソッドの呼び出しはこれ何をどう呼んでんだ...？
とか
結局クラスってなんだ（n回目)
とか
もやもやと悩んでいたので、
そういう「そもそもRって何がどうなっている言語なのか」について向き合える本だと思っています。
（目をそらしていたともいう）
のちのちゆっくりと書きます。
あとあと
岩波DSも2巻が発売されたようで、今回は自然言語処理周りとのこと。
- 作者: 岩波データサイエンス刊行委員会
- 出版社/メーカー: 岩波書店
- 発売日: 2016/02/17
- メディア: 単行本（ソフトカバー）
- この商品を含むブログ (2件) を見る
欲しいな..
積読したくないので来月かな･･･
