次元の海で溺れる

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

【leaflet,jsonlite】家を探しのススメ~路線が分からない編~

あけましておめでとうございます。
新年早々水道を止められかけるなどしましたが元気です。

今回のテーマ

土地勘が無いところの家探しをすると仮定します。
世の中には様々な不動産サイトがありますが、
土地勘が無いともうなんだか全部どうでもよくなってきます。
そもそも住みたいエリアも分からんそもそも知らん。となってきます。

せめてエリアを決めるために以下を仮定して、なんとなく便利げなエリアを把握しよう、というのが
今回の目標です。

  • 目的地は市ヶ谷 (特に意味はない)
  • 乗り換えたくない
  • 駅前に住みたい (300M以内)

ではやってみる。

データ

各種路線情報取得のため、今回はこちらのAPIを使用させて頂きました。
データはこれだけ。
www.ekidata.jp

市ヶ谷駅に乗り入れている路線は

の4つらしい。なるほどわからん。

パッケージ

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))

f:id:WAFkw:20170110005037p:plain

事前に、colorFactor()で各路線をdomainとして定義しておき、
leaflet内では路線ごとにFilterして、それぞれ線を引いたり、駅に円を置いたりする。

addCircles()内のradius引数はm単位なので、今回の「300m」指定に合わせる。
これで駅から300m範囲に円が作れる。

後はaddLayersControl()で各路線を表示ON/OFF出来るようにすればOK。

都営新宿線南北線~とか
f:id:WAFkw:20170110005654p:plain

総武線と有楽町線~とか
f:id:WAFkw:20170110005803p:plain

チェックをポチポチすればいいので便利。便利だ!!!

おわり

あとはこれに各駅の乗降者数とかを乗っければ混む駅とかが分かる。はず。
平均家賃はスクレイピング要りそうでめんどくさいからやめよう。

後はもうすこし情報載せて駅候補決めて、
不動産へ行って「ぼくこの円の中に住みたいんだよね!」って言えばOK。

.......なんで目的地市ヶ谷にしたんだろう.....謎....

参考

addLayersContoroll()を使ってみたくて書いた記事ですが
冷静に考えたらすでにぞうさん先生が大分以前に台風プロットで使ってらっしゃいました。
分かりやすい.....

leafletではじめるRによる地図プロット

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付与されています。

f:id:WAFkw:20161025194358p:plain

小さな都市を作る

前も同じようなことをしたので割愛。
wafdata.hatenablog.com

何のことはない、これを
f:id:WAFkw:20161025194358p:plain

こうするだけ。
f:id:WAFkw:20161025194703p:plain

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に。
これを

f:id:WAFkw:20161025194703p:plain


こうする。

f:id:WAFkw:20161025195753p:plain

これをやるためには、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
#・・・長いので省略

これでどんどん領地を拡大していけばかのアッバース朝みたいな国を作れる。

おわり

  • 何に使うのかは秘密。
  • 出だし何も考えずに、各オブジェクト内結合で全部ID=1にしてたら"non-uniqueID"だよって怒られた。そりゃそうだ。
  • spCbind()もある。
  • 冷静に考えたらウマイヤ朝とかモンゴル帝国とかの方が大きいので目標として適切ではない。

札幌は雪が降ったので引きこもりには最適な環境です。

【orloca,leaflet】ラーメン屋さんになりたい ~SappoRoR#6反省会~

先日SappoRoRへ行ってきて、
おもむろに「ラーメン屋さんになりたい」と言ってきました。

田舎の両親、ラーメン屋さんのみなさま、アホな娘でごめんなさい。
あまり後悔はしてません。

資料各種

  • 改行がダサかったので上げ直し。

概要

テーマ


  • ラーメン屋さんになりたい
  • 札幌市中央区で最も立地に適した場所を探す

データ準備

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))

f:id:WAFkw:20160716215334p:plain

各ポリゴンの代表点を抽出して「世帯」と見立てる

いつもの如く、ワタシのような一般ピーポーには個別世帯のデータなど手に入らないので
各ポリゴンの代表点を取得し、「世帯」とする。
コロプレス図に描画した各小地域の人口を「各世帯の世帯人口」としてあてはめ

    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")

f:id:WAFkw:20160716215402p:plain

立地分析(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)

f:id:WAFkw:20160716215425p:plain

こわすぎ。

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"))

f:id:WAFkw:20160716215449p:plain

改めて結果を可視化

ここがあまりきれいなコードじゃなくて納得いってない..... 各世帯の座標→ラーメン屋 の線をそれぞれ世帯数分引いている。なんかいいやり方はないものか。

    #-立地分析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!!!")

f:id:WAFkw:20160716215510p:plain

まとめ

算出した「here!」に居た先客とは!!
...というくだりには触れないでおく...札幌ネタ....

LTにはちょうどいいくらいの内容だったはずなのに、
他のみなさまがまさかのrevealjsスライド祭りで大変な汗をかきました。
次は10月なのかな。お待ちしてます。

【sp,deldir,leaflet】あなたの全てを知りたい編~Polygon内に含まれるポイントデータを集計~

GWなので家出をしました。
帰りたくないけどおうちに帰るまでが家出だからなあ...

本日のテーマ

掲題の通りです。

あなた(Polygon)の全部が知りたいのよ!という気持ちなので
Polygon内に含まれる複数の点の情報を集計するなどしたいです。

たとえばさ、
算出したポリゴン内に複数の世帯が含まれていたとして、
そのポリゴン内の人口を算出したい、みたいなの。

f:id:WAFkw:20160503220019p:plain

データ

そうはいっても不条理な世の中では私のような下々のパンピーが世帯別のデータを手に入れることは出来ません。
乱数作ってポリゴン内に配置して~とかが主流なのでしょうが、
妄想でもいいのでもう少しリアルに近いことがしたい。

ということで以下の通りに。

#------------------------
①札幌市中央区の小地域データからポリゴン内代表点を抽出。
     →「世帯」だと思い込んで取扱い
②札幌市中央区内の小学校の座標を用いて中央区をボロノイ分割
     →集計したい「ポリゴン」を作成
③②に含まれる①の点を集計してポリゴン内人口を算出する。
#------------------------

イメージは「小学校区(ボロノイ領域基準)内の人口を把握する」という分析をする感じです。
ざっくり同じくらいの人口になるんじゃないかなと仮説を立てつつ。

妄想の世界へ

①、②は過去にやったやつの焼き直し的なコードなので細かい記述は省きます。
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  |北十六条西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|
|6  |北十六条西15丁目 | 141.3304| 43.07795|01101791615 |   142|    80| 141.3304| 43.07795|
|7  |北十五条西15丁目 | 141.3306| 43.07682|01101791515 |   285|   197| 141.3306| 43.07682|
|8  |北十六条西21丁目 | 141.3209| 43.07592|01101791621 |    41|    20| 141.3209| 43.07592|
|9  |北十四条西15丁目 | 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"))

f:id:WAFkw:20160503223958p:plain

この黒いツブツブが代表点で、各ポリゴン一つずつ。
元データに小地域ごとの人口が入っているので、そこと紐づけて世帯+世帯人口として取り扱います。

次にボロノイ分割で領域作り。
小学校の座標データがいまいち取れなかったので作りました。

#ボロノイ

#データの読み込み
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)

f:id:WAFkw:20160503222734p:plain

#ボロノイ領域のポリゴン化
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")

f:id:WAFkw:20160503223835p:plain
f:id:WAFkw:20160503223718p:plain

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))

f:id:WAFkw:20160503224917p:plain

できた。満足。

感想

山間(左側)は人口が他に比べて少なめ。学校のボロノイ領域が広いのも元々人口少なめなのが関係してるのかも。
やはり中心部は人口多いですねー小学校も生徒人数多かったりするのかしら。

詳しく調べるには世帯の年齢とかも加味しないとダメそうなのでこれ以上のコメントはやめておく.....

おわり

先日久々のTokyoRへ行ってきました。次はいつ来られるかなあ
ちなみにSappoRoRは7月9日です。ぜひ。

本州は暑くてね.......

【maptools,leaflet,rgdal】 (愛と希望で)世界を一つにまとめたい

何の予定も無いけど、なんだか突然LTがしたくなったので
自宅.Rをやるべく資料作りに勤しむ遊びにハマっています。

待てど暮らせど北国に春が来ない。

今回のテーマ

「世界を一つにまとめたい(意訳:複数ポリゴンを任意に結合したい)」

LT資料(自宅.R用)を作るべく
国勢調査の小地域データを手に入れたので、その可視化と取扱いについてまとめつつ。

データセット

平成22年 国勢調査(小地域)境界データ
(地図で見る統計(統計GIS) 政府統計の総合窓口 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 "旭ケ丘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" ...


#変数が多すぎなので、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|
|北十六条西16丁目 | 141.3260| 43.07778|01101791616 |     0|     0|
|北十八条西         | 141.3298| 43.08034|011017918   |   136|    72|
|北十七条西         | 141.3302| 43.07896|011017917   |   149|    72|
|北十六条西15丁目 | 141.3304| 43.07795|01101791615 |   142|    80|
|北十五条西15丁目 | 141.3306| 43.07682|01101791515 |   285|   197|
|北十六条西21丁目 | 141.3209| 43.07592|01101791621 |    41|    20|
|北十四条西15丁目 | 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")

f:id:WAFkw:20160313202014p:plain
f:id:WAFkw:20160313211024p:plain

下の図鳥肌が立つ。。。。。


ここでちょっと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を指定

f:id:WAFkw:20160313203017p:plain


小地域だとなんか線がいっぱいになって見づらい。
でも外線一切なしだと儚げな感じになっちゃってそれもそれで見づらい。

→一番外側(丁目単位ではなく、札幌市中央区の外側)だけ色指定したい!!!!

そこでポンコツはかんがえた

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)

f:id:WAFkw:20160313204516p:plain

小地域の区切りが消えてすっきり。


重ねて描いてみる

最初に描いてたコロプレスと、さっき結合したものを重ねて外側だけ線を引く。

f:id:WAFkw:20160313205549p:plain

#コロプレスと重ねる
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")

f:id:WAFkw:20160313205752p:plain
f:id:WAFkw:20160313210055p:plain

できた。

閑話休題

世界をつなげる力を手に入れた。

そんなことより、@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)

f:id:WAFkw:20160313210453p:plain

なんか左上にズレる.....。なんで....。どうしてなのパトラッシュ....
ここで言ってる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("納豆","豆腐"))

f:id:WAFkw:20160221215205p:plain


簡単にすっきり並べられた!!

....けどこう...なんかこう...デフォルトの色味がこう...いかがわしい。
なんかとてもいかがわしい色味。

ということでもう少し手を入れてみる。

消費額を離散化して、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"))

f:id:WAFkw:20160221215645p:plain

清純な感じになりました。
親御さんもこれで一安心です。

結果については前回記事でもグダグダ言ってましたが、
納豆の地域差が東西ではっきり出ていて、豆腐はそんなに地域差は無さそう,,

早速描く~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)

f:id:WAFkw:20160221222900p:plain
f:id:WAFkw:20160221223004p:plain
f:id:WAFkw:20160221223029p:plain

描画のところのコードが汚くなってしまった。猛省。。
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)

f:id:WAFkw:20160221224605p:plain
f:id:WAFkw:20160221224624p:plain

※動きません


leafletは描くのも楽だし見た目も綺麗だし
少し透過させるとTileの地図情報も読み取れるのがいいですね。



ちからつきる

じつはcoloplethrパッケージで描くのもやりたかったんだけど、
ほんじつはここでちからつきた

ゆうしゃろと、ここにねむる。


P.S 書籍発売ラッシュ

R言語徹底解説

R言語徹底解説

かのHadley神の「AdvancedR」邦訳版がついに。


個人的には7章が好きです。
7章が好きです。

Rをはじめたての頃は、結果が返ればまあいいや、だったところ、
だんだん関数のソースを見てごにょごにょするようになってくると

ここのメソッドの呼び出しはこれ何をどう呼んでんだ...?
とか
結局クラスってなんだ(n回目)
とか

もやもやと悩んでいたので、
そういう「そもそもRって何がどうなっている言語なのか」について向き合える本だと思っています。

(目をそらしていたともいう)

のちのちゆっくりと書きます。



あとあと

岩波DSも2巻が発売されたようで、今回は自然言語処理周りとのこと。

岩波データサイエンス Vol.2

岩波データサイエンス Vol.2

欲しいな..
積読したくないので来月かな・・・

【TeachingDemos】誕生日なので大豆について考えてみた(2)

n+1歳したよポエムは後に回すとして、

世は節分なので大豆について考える。パート2。


前回
wafdata.hatenablog.com

テーマ

都道府県別の大豆データがあったので
サブプロットとして地図に埋め込んでみる。

実はちょっと夢だったサブプロット。わくわく。

使用データ

農林水産省/大豆関連データ集)より、


都道府県別生産状況(大豆収穫量(t))
・豆腐・納豆の都市別購入順位(豆腐都市別購入金額・納豆都市別購入金額)

このへんをお借りしてやってみる。

準備

シェープファイルはNipponパッケージのフォルダに含まれているやつを借りる。
CRSクラス作成では投影法を地球座標系、準拠楕円体をWGS84に指定した。

>map <- system.file("shapes/jpn.shp", package = "Nippon")[1] %>% 
   maptools::readShapePoly(., proj4string = CRS("+proj=longlat +datum=WGS84"))

> map@data
   SP_ID jiscode      name population           region
0      1      01  Hokkaido    5506419         Hokkaido
1      2      02    Aomori    1373339           Tohoku
2      3      03     Iwate    1330147           Tohoku
3      4      04    Miyagi    2348165           Tohoku
4      5      05     Akita    1085997           Tohoku
5      6      06  Yamagata    1168924           Tohoku
6      7      07 Fukushima    2029064           Tohoku
7      8      08   Ibaraki    2969770            Kanto
8      9      09   Tochigi    2007683            Kanto
9     10      10     Gunma    2008068            Kanto
10    11      11   Saitama    7194556            Kanto
11    12      12     Chiba    6216289            Kanto
12    13      13     Tokyo   13159388            Kanto
13    14      14  Kanagawa    9048331            Kanto
14    15      15   Niigata    2374450            Chubu
15    16      16    Toyama    1093247            Chubu
16    17      17  Ishikawa    1169788            Chubu
17    18      18     Fukui     806314            Chubu
18    19      19 Yamanashi     863075            Chubu
19    20      20    Nagano    2152449            Chubu
20    21      21      Gifu    2080773            Chubu
21    22      22  Shizuoka    3765007            Chubu
22    23      23     Aichi    7410719            Chubu
23    24      24       Mie    1854724            Chubu
24    25      25     Shiga    1410777            Kinki
25    26      26     Kyoto    2636092            Kinki
26    27      27     Osaka    8865245            Kinki
27    28      28     Hyogo    5588133            Kinki
28    29      29      Nara    1400728            Kinki
29    30      30  Wakayama    1002198            Kinki
30    31      31   Tottori     588667          Chugoku
31    32      32   Shimane     717397          Chugoku
32    33      33   Okayama    1945276          Chugoku
33    34      34 Hiroshima    2860750          Chugoku
34    35      35 Yamaguchi    1451338          Chugoku
35    36      36 Tokushima     785491          Shikoku
36    37      37    Kagawa     995842          Shikoku
37    38      38     Ehime    1431493          Shikoku
38    39      39     Kochi     764456          Shikoku
39    40      40   Fukuoka    5071968 Kyushu / Okinawa
40    41      41      Saga     849788 Kyushu / Okinawa
41    42      42  Nagasaki    1426779 Kyushu / Okinawa
42    43      43  Kumamoto    1817426 Kyushu / Okinawa
43    44      44      Oita    1196529 Kyushu / Okinawa
44    45      45  Miyazaki    1135233 Kyushu / Okinawa
45    46      46 Kagoshima    1706242 Kyushu / Okinawa
46    47      47   Okinawa    1392818 Kyushu / Okinawa

data部抜くとこんな感じ。
ちなみに人口も地域も入っているので相当便利。

で、
このデータ部の行番号の順を動かさずに、大豆のデータをくっつけたいので、
nameをキーにしてJOINすればいいかなーという方針を立てる。


各種大豆データを
PDF→CSVして、都道府県名をローマ字にするなどすると

>mame_sub<-read.csv("0203mame.csv",header = TRUE)
>knitr::kable(mame_sub)

|name      | natto| tofu| mame_t|
|:---------|-----:|----:|------:|
|Hokkaido  |  3847| 4362|  73600|
|Aomori    |  4807| 5140|   5370|
|Iwate     |  4654| 6945|   5470|
|Miyagi    |  5485| 6264|  19300|
|Akita     |  4477| 5123|   9640|
|Yamagata  |  5569| 5906|   7720|
|Fukushima |  5298| 5254|   2220|
|Ibaraki   |  5916| 4839|   5410|
|Tochigi   |  3962| 5111|   4250|
|Gunma     |  4808| 5226|    444|
|Saitama   |  4406| 6328|    786|
|Chiba     |  4064| 6151|   1080|
|Tokyo     |  3618| 5969|      4|
|Kanagawa  |  3418| 5820|     64|
|Niigata   |  4586| 5070|   8840|
|Toyama    |  4329| 4927|   7630|
|Ishikawa  |  3505| 5176|   2310|
|Fukui     |  3238| 4815|   2430|
|Yamanashi |  4501| 5274|    289|
|Nagano    |  4157| 4510|   3380|
|Gifu      |  3160| 5078|   3430|
|Shizuoka  |  3735| 6057|    317|
|Aichi     |  3348| 5350|   7180|
|Mie       |  2429| 5307|   3750|
|Shiga     |  2912| 5059|   9450|
|Kyoto     |  2680| 5982|    455|
|Osaka     |  1945| 5470|     20|
|Hyogo     |  2289| 5709|   3320|
|Nara      |  2421| 5408|    256|
|Wakayama  |  1994| 4864|     43|
|Tottori   |  2658| 6728|   1120|
|Shimane   |  2605| 6213|   1320|
|Okayama   |  2428| 4534|   2370|
|Hiroshima |  2592| 5258|    667|
|Yamaguchi |  2744| 4971|   1020|
|Tokushima |  2306| 6763|     45|
|Kagawa    |  2528| 5425|     97|
|Ehime     |  2628| 6170|    457|
|Kochi     |  2017| 6009|     61|
|Fukuoka   |  3125| 4563|  14300|
|Saga      |  3317| 6053|  15300|
|Nagasaki  |  2672| 4824|    496|
|Kumamoto  |  3861| 5417|   3710|
|Oita      |  3604| 6288|   1690|
|Miyazaki  |  2937| 5239|    317|
|Kagoshima |  3513| 5471|    301|
|Okinawa   |  2389| 7478|      1|

#natto<-納豆購入額
#tofu<-豆腐購入額
#mame_t<-大豆収穫量(t)


あとはJOINして準備するだけ。

#JOINしたものをデータ部に置き換え
map@data<- dplyr::left_join(map@data,mame_sub,"name") 

#各都道府県の代表点抽出
>map.c<-coordinates(map) %>% 
  as.data.frame(.)
>knitr::kable(map.c)


|   |       V1|       V2|
|:--|--------:|--------:|
|0  | 142.5637| 43.37322|
|1  | 140.8369| 40.78999|
|2  | 141.3672| 39.59491|
|3  | 140.9347| 38.45849|
|4  | 140.4040| 39.75225|
|5  | 140.1059| 38.44662|
|6  | 140.2207| 37.38554|
|7  | 140.3268| 36.33433|
|8  | 139.8238| 36.68965|
|9  | 138.9818| 36.51044|
|10 | 139.3504| 35.99388|
|11 | 140.2134| 35.53588|
|12 | 139.4481| 35.69395|
|13 | 139.3502| 35.40736|
|14 | 138.9373| 37.46545|
|15 | 137.2517| 36.63167|
|16 | 136.7571| 36.73964|
|17 | 136.2342| 35.83485|
|18 | 138.6130| 35.61755|
|19 | 138.0369| 36.12952|
|20 | 137.0573| 35.78676|
|21 | 138.3440| 35.00518|
|22 | 137.2163| 35.04840|
|23 | 136.3712| 34.54634|
|24 | 136.1268| 35.22382|
|25 | 135.4460| 35.24641|
|26 | 135.5121| 34.62303|
|27 | 134.8311| 35.08439|
|28 | 135.8637| 34.31070|
|29 | 135.5186| 33.90602|
|30 | 133.8624| 35.36655|
|31 | 132.4981| 35.00385|
|32 | 133.8122| 34.89499|
|33 | 132.7951| 34.63696|
|34 | 131.5512| 34.20905|
|35 | 134.2566| 33.91387|
|36 | 133.9805| 34.20724|
|37 | 132.8611| 33.60891|
|38 | 133.3812| 33.41964|
|39 | 130.6720| 33.53417|
|40 | 130.1211| 33.28651|
|41 | 129.9260| 32.95754|
|42 | 130.8424| 32.64685|
|43 | 131.4390| 33.18740|
|44 | 131.2934| 32.18176|
|45 | 130.6119| 31.62818|
|46 | 127.9566| 26.48630|

cols<-c("#E06A3B","#FBE481") #色スタンバイ

抽出した代表点に、各サブプロットを置いていくイメージです。

棒グラフで豆腐と納豆を比較する

分かりそうなことは、
豆腐と納豆の購入額を純粋に比較して、地域差がどう出るか、みたいなとこでしょうか。
ただの差の比較にしかならないけど。

#元となるplot
plot(map,col=gray(0.8))


#サブプロット用データと、Y軸の最大値決めるためのMax値
>map.d<-cbind(map@data$natto,map@data$tofu)
>map.m<-max(map.d)

#各ポリゴンの代表点の位置にsubplotをする
for(i in 1:nrow(map.c)){
  subplot(barplot(map.d[i,],ylim=c(0,map.m),col=cols,
                  names=c("",""),yaxt="n"),x=map.c[i,1],
          y=map.c[i,2],size=c(0.1,0.5),vadj=0)
}

locator(2)
legend(139,29,c("納豆","豆腐"),fill=cols,cex=0.8)

f:id:WAFkw:20160203001022p:plain

(forで回しちゃったけどいいやり方あるかしら...)


Y軸を共通にしているので一応単純比較出来るはず。

これが結構おもしろい。


圧倒的西日本の豆腐消費がまず目につく。
が、
圧倒的納豆消費の少なさも目に付く。

東日本は豆腐消費はほどほど、
納豆もほどほど、といった感じ。

これはそれぞれコロプレス図みたいにしてもパキっと分かれて綺麗そう。

特に納豆。納豆。

近畿・中国・四国は特に少ないけどたまたまかなあ、
経年でデータ持ってきてみれば良かった。失敗。

円グラフもしてみる

円グラフもかいてみた

plot(map,col=gray(0.8))

for(i in 1:nrow(map.c)){
  subplot(pie(cbind(map@data$natto[i],map@data$tofu[i]),col=cols,labels=c("","")),
          x=map.c[i,1],y=map.c[i,2],size=c(0.3,0.3),vadj=0)
}
legend(139,29,c("納豆","豆腐"),fill=cols,cex=0.8)

f:id:WAFkw:20160203001959p:plain


なんか凡例がうまく表示されない。はて。

円でも黄色比重が高まっているエリアを見つけられる。
このデータの作りだとあまり円にする意味は無かったかな...

円の大きさを大豆の収穫量で変化するようにしてみた

第3の要素として、ここに各都道府県の大豆の収穫量をちゃっかり入れて、
円グラフの半径として使う。

大豆収穫量の多いところと消費の関連とか掴めたらおもしろいなーと思ったんだけど...

#サイズ変える
r<-sqrt(map@data$mame_t/pi) #円の大きさは大豆収穫量/パイの二乗根
r<-r/max(r)*0.4

plot(map,col=gray(0.8))

for(i in 1:nrow(map.c)){
  subplot(pie(cbind(map@data$natto[i],map@data$tofu[i]),col=cols,labels=c("","")),
              x=map.c[i,1],y=map.c[i,2],size=c(r[i],r[i]),vadj=0)
}
legend(139,29,c("納豆","豆腐"),fill=cols,bty="n",cex=0.8)

f:id:WAFkw:20160203002446p:plain


...世の中そううまくはいかなかった。

作りなよ。みんなもっと大豆作りなよ。ってくらいの収穫量してた。
何も見えやしない。

豆の話(2)はおしまい

後味は悪いが今回はここまで。
まだ大豆データがあるから何か閃いたら(3)としてやろうかな



誕生日

あんまりめでたくないけどせっかく年が一つ増えて最初なので
なんか振り返りながらつらつらと書いてみる。

振り返る

・卒業:卒論を提出後完全にスイッチが切れて日本各地をぷらぷら
    大阪→京都 →秋田→(どこに居たか覚えてない10日)→
    気付いたら埼玉の某所で寿司食べてた

・仕事:こんなポンコツ根暗でも上司に恵まれてなんとか仕事が出来ている。
    頭が上がらない。

・趣味:卒業したらRを触る「理由」みたいなものが無くなってしまったので、
    趣味と開き直ってブログを始めたり、
    各所のRコミュニティーに勇気を振り絞って出たりした。
    塩を撒かれて追い返されるんだと思ってたのに、
    気が付いたらR繋がりで知り合いの方が増えてきて。
    
    Rについてはもっとわかるようになりたい。
    何をもってわかるとするのか、みたいのはあるけど。
    わかりたい、とは思う。し、
    そろそろ逃げない分析をやらなきゃなあ、とも思う。IRとかIRとか。

100万円貯まったら次の街へ行く


北の大地を出ることにした。


と、結構な決断をしたつもりだったのに、
どうやら自分から出ていかなくてもそうなりそうな予感。
多分近くも遠くもないうちに。

元々縁もゆかりも無かったけど、今じゃ永住出来ると思えるくらい一番好きな街なので、
札幌は心からおすすめ。

住んでるうちに変わるものも変わらないものもあって、
5年間少し穏やかに過ごし過ぎたので、次の街ではちょっと苦しみたいなというのが希望です。


まだもうしばらくは居る。予定。


ちなみにこれに感化されたわけではないけど良い映画。
とにもかくにも主題歌が名曲です。

百万円と苦虫女 [DVD]

百万円と苦虫女 [DVD]

おわり。