【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が何をどうやって求めてるのか分かればすっきり出来るのかしら。
おわり
ねむいのでおわります。
足りないものが足りないままなのすごくいらいらするなあ、というまいにち。
