読者です 読者をやめる 読者になる 読者になる

次元の海で溺れる

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

【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日です。ぜひ。

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