保育園合格戦略×データサイエンス

on 2月 4, 2021

leaflet

保育園の合格戦略をデータサイエンスから考える

GISを使って攻略したいと思います。

今回は、東京都で待機児童が一番多い（後記　2020年には待機児童ゼロを達成されたそうです）と言われている世田谷区のデータを使って解析したいと思います。最終的に見たいデータは、競争率です。保育園の半径500ｍ円の範囲中の人口（０～４歳）と、保育園の収容人数を使って競争率を計算します。

GIS的にポイントとなるのは、半径500m円の中の人口を計算するところです。保育園の500m範囲の人口のデータそのものは存在しません。そのため、同じ範囲で人口データを有している行政境界データを用いて、交差した面積に応じて人口を計算します。この計算は、Areal weighted interpolationもしくは面積按分法と言われています。

使用データ

e-stat：人口統計データ　

１― 国勢調査　2015年　小地区　統計データ　年齢（５歳階級、４区分）別、男女別人口　13 東京都

２― 国勢調査　2015年　5次メッシュ（250mメッシュ）統計データ　年齢（５歳階級、４区分）別、男女別人口　M5339

e-stat：上記データのshape file

３― 国勢調査　2015年　小地区　境界データ　13112 世田谷区世界測地系緯度経度・Shapefile

４― 国勢調査　2015年　5次メッシュ（250mメッシュ）境界データ　13112 世田谷区世界測地系緯度経度・Shapefile

政府統計の総合窓口(e-Stat)（https://www.e-stat.go.jp/）を加工して作成

国土交通省：保育園の位置データ

５― 国土数値情報　福祉施設データ　世界測地系　平成27年

「国土数値情報（福祉施設データ）」（国土交通省）（https://nlftp.mlit.go.jp/ksj/gml/datalist/KsjTmplt-P14.html#prefecture13）を加工して作成

地図上に人口コログラムと保育園を描写

ライブラリtidyverseとsfを読み込みます。また、日本語フォントも読み込みます。

library(tidyverse)
library(sf)
windowsFonts(meiryo = windowsFont("Meiryo UI"))

まずは、データ３の世田谷区の行政境界データを読み込みます。

setagaya_gyosei <- st_read("gyoseikyokai/h27ka13112.shp",
                    options = "ENCODING=CP932")

つぎは、データ１の東京の人口データを読み込みます。そして、上記の世田谷行政境界データとKEY_CODEを使って結合させます。T000849002列に０～４歳人口が記載されています。

pop_gyosei <- read.delim("tokyo_jinko_gyosei.txt",
                         sep=",",
                         encoding = "CP932") 

# T000849002が0～4歳人口
pop_gyosei <- pop_gyosei %>% select(KEY_CODE,T000849002)

pop_gyosei$KEY_CODE <- as.character(pop_gyosei$KEY_CODE)

setagaya_gyosei <- left_join(setagaya_gyosei,pop_gyosei)

setagaya_gyosei$T000849002 <- as.numeric(setagaya_gyosei$T000849002)

さっそくggplotで地図を書いてみます。

ggplot()+
  geom_sf(data = setagaya_gyosei,
          aes(fill=T000849002))+
  scale_fill_gradientn(colors=c("deepskyblue","lightgreen","palevioletred","red"))+
  theme_void()+
  labs(fill="0~4歳人口")+
  theme(legend.text=element_text(family = "meiryo"))

これでは、どこかよくわからないのでベースマップを付けます。ベースマップはopenstreetmapを使います。ライブラリggspatialを読み込みます。

距離を示す棒（スケールバー）も付けました。

library(ggspatial)
ggplot()+
  annotation_map_tile(zoomin = 0) +
  geom_sf(data = setagaya_gyosei,
          aes(fill=T000849002),
          alpha=0.5)+
  scale_fill_gradientn(colors=c("deepskyblue","lightgreen","palevioletred","red"))+
  annotation_scale(location = "bl", width_hint = 0.5) +
  labs(caption = "\U00a9 OpenStreetMap contributors") +
  theme_void()+
  labs(fill="0~4歳人口")+
  theme(legend.text=element_text(family = "meiryo"))

データ５の福祉施設のポイントデータが入ったデータを読み込みます。この中に保育関係のデータが含まれています。列名P114_006の番号800～899が保育関係です。そして、世田谷区だけ抽出します。

fukushi <- st_read("fukushi/P14-15_13.shp",
                   options = "ENCODING=CP932")

hoiku <- fukushi %>% 
  filter(P14_006>=800) %>% 　　#800番台が保育関係
  filter(P14_006<900) %>% 
  filter(P14_002=="世田谷区")　　#世田谷区のみ抽出

さきほど作った地図の上にポイントデータを載せて描写します。

ggplot()+
  annotation_map_tile(zoomin = 0) +
  geom_sf(data = setagaya_gyosei,
          aes(fill=T000849002),
          alpha=0.5)+
  scale_fill_gradientn(colors=c("deepskyblue","lightgreen","palevioletred","red"))+
  geom_sf(data = hoiku)+
  annotation_scale(location = "bl", width_hint = 0.5) +
  labs(caption = "\U00a9 OpenStreetMap contributors") +
  theme_void()+
  labs(fill="0~4歳人口")+
  theme(legend.text=element_text(family = "meiryo"))

とりあえず、人口コログラムと保育園をプロットできました。

保育園を中心とした半径500ｍの円を描写

保育園から半径500ｍ以内の人口を計算するために、円を作ります。距離を扱う場合には、ライブラリunitsを読み込むのと、crsを直交系に変更する必要があります。東京なのでCRSはEPSG：2451にしました。（JGD2000。。）

st_bufferでポイントからの円を作れます。

library(units)
radius <- set_units(0.5,km)

# 直交系に変換
hoiku <- st_transform(hoiku, crs = 2451)
# 円を作る。
hoiku_500m <- st_buffer(hoiku, dist = radius)
# 直交系から緯度経度に戻す（WGS84）
hoiku_500m <- st_transform(hoiku_500m, crs = 4326)

とりあえず円を書いてみます。

ggplot()+
  annotation_map_tile(zoomin = 0) +
  geom_sf(data = setagaya_gyosei,
          aes(fill=T000849002),
          alpha=0.5)+
  scale_fill_gradientn(colors=c("deepskyblue","lightgreen","palevioletred","red"))+
  geom_sf(data = hoiku_500m,alpha=0.2)+
  annotation_scale(location = "bl", width_hint = 0.5) +
  labs(caption = "\U00a9 OpenStreetMap contributors") +
  theme_void()+
  labs(fill="0~4歳人口")+
  theme(legend.text=element_text(family = "meiryo"))

半径500m円内の人口を計算（Areal weighted interpolation面積按分法）

最初の４行は、必要なデータのみを選抜しています。

次の行は、st_areaを使ってポリゴンの面積を計算しています。新たな列をつくってそこに面積を代入します。

hoiku_500m_selected <- hoiku_500m %>%
  select(P14_007,P14_009,geometry) #名称、収容人数、geometory

setagaya_selected <- setagaya_gyosei %>%
  select(S_NAME,T000849002,geometry) #地域名、０～４歳人口、geometry

setagaya_selected$area <- st_area(setagaya_selected) # 面積計算

hoiku_500m_selected <- st_transform(hoiku_500m_selected,　
                     crs = st_crs(setagaya_selected))　# crsを合わせる

intersect_500m <- st_intersection(hoiku_500m_selected,  # 重なりを計算
                            setagaya_selected)

intersect_500m$now_area <- st_area(intersect_500m)　#　面積を計算

次は、重なりを使って新たなシェープファイルを作ります。

st_intersectionを使って重なりを検出します。しかしCRSが揃ってないといけないので、その前にCRSの変更をしています。

st_intersection後は新たなシェープファイルが出来ましたが、ここでもまた面積を計算します。

hoiku_500m_selected <- st_transform(hoiku_500m_selected,　
                     crs = st_crs(setagaya_selected))　# crsを合わせる

intersect_500m <- st_intersection(hoiku_500m_selected,  # 重なりを計算
                            setagaya_selected)

intersect_500m$now_area <- st_area(intersect_500m)　#　面積を計算

intersect_500m$areal_weight <- intersect_500m$now_area/intersect_500m$area

intersect_500m$population_0to4 <- intersect_500m$areal_weight*intersect_500m$T000849002 #行政エリアの0~4歳人口をかける

intersect_500m <- intersect_500m %>%
  group_by(P14_007,P14_009) %>%　#名前と収容人数を残す
  summarise(population_0to4 = sum(population_0to4,
                      na.rm = T))

４コード目がareal weightを計算しているところです。行政境界のときの面積と、今の面積の割合を計算します。

５コード目が、面積按分を行っているところです。行政境界に住んでいる人数を、面積割合に応じて分配しています。

６コード目は、円が複数の行政境界にまたがっているので、その境界ごとに計算された人数を合計しているところです。group_byのP14_007が保育園の名前です。P14_009は収容人数で、グループ化に重要ではなく、ただ変数を残したいからgroup_byに入っています。

円内の人口をわかるようにプロットします。その前に、population_0to4には、面積計算からの流れで単位が入っています。その単位があるとグラフが書けないのでただの数字にします。

#単位が入った形からただの数字にする
intersect_500m$population_0to4 <- as.numeric(intersect_500m$population_0to4)

ggplot()+
  annotation_map_tile(zoomin = 0) +
  geom_sf(data = intersect_500m,
          aes(fill=population_0to4),
          alpha=0.5)+
  scale_fill_gradientn(colors=c("deepskyblue","lightgreen","palevioletred","red"))+
  annotation_scale(location = "bl", width_hint = 0.5) +
  labs(caption = "\U00a9 OpenStreetMap contributors") +
  theme_void()+
  labs(fill="保育園の半径500m圏内の0~4歳人口")+
  theme(legend.text=element_text(family = "meiryo"))

競争率を計算

０～４歳人口÷収容人数で競争率を計算します。

なぜか収容人数がマイナスなのがあるので最初に除外しています。


kyosou <- intersect_500m %>% 
  filter(P14_009>0)

kyosou$kyosouritu <- kyosou$population_0to4/kyosou$P14_009

ggplot()+
  annotation_map_tile(zoomin = 0) +
  geom_sf(data = kyosou,
          aes(fill=kyosouritu),
          alpha=0.5)+
  scale_fill_gradientn(colors=c("deepskyblue","lightgreen","palevioletred","red"))+
  annotation_scale(location = "bl", width_hint = 0.5) +
  labs(caption = "\U00a9 OpenStreetMap contributors") +
  theme_void()+
  labs(fill="保育園の半径500m圏内の0~4歳人口")+
  theme(legend.text=element_text(family = "meiryo"))

円がでかすぎてわかりにくいので、保育園のポイントデータを使って競争率を表す。

hoiku <- left_join(hoiku,
                   as_tibble(kyosou),
                   by="P14_007")

hoiku$kyosouritu <- as.numeric(hoiku$kyosouritu)

ggplot()+
  annotation_map_tile(zoomin = 0) +
  geom_sf(data = hoiku,
          aes(color=kyosouritu),
          size = 3,
          shape=15)+
  scale_color_gradientn(colors=c("deepskyblue","lightgreen","palevioletred","red"))+
  annotation_scale(location = "bl", width_hint = 0.5) +
  labs(caption = "\U00a9 OpenStreetMap contributors") +
  theme_void()+
  labs(fill="保育園の半径500m圏内の0~4歳人口")+
  theme(legend.text=element_text(family = "meiryo"))

今度は、データ２と４メッシュデータを使って解析してみます。メッシュデータは、250mの正方形ごとに人口が計算されているのでもっと正確です。しかし、０～４歳人口のデータがないので６歳未満世帯員のいる一般世帯数をつかいます。

データ４のメッシュデータを読み込みます。しかし、このメッシュは東京全体をすべて網羅するものなので世田谷に限定します。行政境界と同じ範囲のメッシュに限定したいと思います。メッシュには、土地の名称等のデータがないので地図上で重なりを判定しなければなりません。

st_intersectsを使います。さきほど使ったのは、st_intersectionsです。st_intersectsとst_intersectionsの違いは、st_intersectsは重なっているかTRUEかFALSEで判定されたデータが帰ってきます。一方、st_intersectionsは重なったシェープファイルが帰ってきます。今回は、メッシュデータに傷をつけたくないので、st_intersectsでロジカル判定列を使ってスライスします。

st_intersectsの返り値は特殊な行列なので、lengthsを使った変換を経て、スライスに使います。

mesh <- st_read("tokyo_mesh\\MESH05339.shp",
                options = "ENCODING=CP932")

mesh <- st_transform(mesh,
                     crs = st_crs(setagaya_gyosei))

mesh_logical <- st_intersects(mesh,setagaya_gyosei)
mesh_logical = lengths(mesh_logical) > 0
mesh_setagaya <- mesh[mesh_logical,]

ggplot()+
  annotation_map_tile(zoomin = 0) +
  geom_sf(data = mesh_setagaya)+
  scale_color_gradientn(colors=c("deepskyblue","lightgreen","palevioletred","red"))+
  annotation_scale(location = "bl", width_hint = 0.5) +
  labs(caption = "\U00a9 OpenStreetMap contributors") +
  theme_void()

データ２の人口メッシュデータを読み込みます。４コード目で、上のメッシュと結合させます。

結合するにあたっては、keyの型を合わせる必要があるのでcharacterに変更しています。T000876037は、６歳未満世帯員のいる一般世帯数の列です。

T000876037は、なぜかcharacter型で読み込まれていたので数にします。GISのデータ読み込みはこの手のエラーが多いので注意しています。グラフが変だなと思ったら、文字だったということがよくあります。

pop_mesh <- read.delim("tokyo_jinko_mesh.txt",
                       sep = ",",
                       encoding = "CP932") 

pop_mesh <- pop_mesh %>% select("KEY_CODE","T000876037")

pop_mesh$KEY_CODE <- as.character(pop_mesh$KEY_CODE)

mesh_setagaya <- left_join(mesh_setagaya,pop_mesh,by="KEY_CODE")

mesh_setagaya$T000876037 <- as.numeric(mesh_setagaya$T000876037) # ６歳未満世帯員のいる一般世帯数

ggplot()+
  annotation_map_tile(zoomin = 0) +
  geom_sf(data = mesh_setagaya,
          aes(fill=T000876037),
          alpha=0.5) +
  scale_fill_gradientn(colors=c("deepskyblue","lightgreen","palevioletred","red"))+
  geom_sf(data = hoiku,
          aes(color=kyosouritu),
          size = 3,
          shape=15)+
  scale_color_gradientn(colors=c("deepskyblue","lightgreen","palevioletred","red"))+
  annotation_scale(location = "bl", width_hint = 0.5) +
  labs(caption = "\U00a9 OpenStreetMap contributors") +
  theme_void()+
  labs(fill="６歳未満世帯員のいる一般世帯数",
       color = "保育園競争率")+
  theme(legend.text=element_text(family = "meiryo"))

ggsave("mesh6.png",
       plot = last_plot(),
       device = "png",
       units="in",
       dpi = 500)

最後に、インタラクティブなleafletを使った地図を作ります。色付けるのに独自のカラーパレットfunctionを作る必要があるのと、ＣＲＳを全て4326にするところがポイントです。

T000849002が0～4歳人口です。

library(leaflet)
setagaya_gyosei <- st_transform(setagaya_gyosei, crs = 4326)
hoiku <- st_transform(hoiku, crs = 4326)

setagaya_gyosei_leaf <- setagaya_gyosei %>% 
  select(T000849002,geometry)

setagaya_gyosei_leaf$T000849002 <- as.numeric(setagaya_gyosei_leaf$T000849002)

hoiku_leaf <- hoiku %>% 
  select(P14_007,kyosouritu)

hoiku_leaf$kyosouritu <- as.numeric(hoiku_leaf$kyosouritu)

pal1 <- colorNumeric(
  palette = "RdBu",
  reverse = T,
  domain = setagaya_gyosei_leaf$T000849002)
pal2 <- colorNumeric(
  palette = "RdBu",
  reverse = T,
  domain = hoiku_leaf$kyosouritu)


leaf <-leaflet() %>% 
  addTiles() %>% 
  addPolygons(data = setagaya_gyosei_leaf,
              stroke = T,
              color = "grey",
              weight = 2,
               fillColor = ~pal1(T000849002)) %>% 
  addLegend(data=setagaya_gyosei_leaf,
            "bottomright",
            pal = pal1,
            values = ~T000849002,
            title = "人口",
            opacity = 0.8) %>% 
  addCircles(data = hoiku_leaf,
             stroke = T,
             radius = 50,
             popup = ~P14_007,
             color = ~pal2(kyosouritu)) %>% 
  addLegend(data=hoiku_leaf,
            "topright",
            pal = pal2,
            values = ~kyosouritu,
            title = "競争率",
            opacity = 0.8)

htmlwidgets::saveWidget(leaf,"leaf_setagaya.html")

Categories:

GIS R その他

保育園の合格戦略をデータサイエンスから考える

使用データ

e-stat：人口統計データ

e-stat：上記データのshape file

国土交通省：保育園の位置データ

５― 国土数値情報 福祉施設データ 世界測地系 平成27年

地図上に人口コログラムと保育園を描写

保育園を中心とした半径500ｍの円を描写

半径500m円内の人口を計算（Areal weighted interpolation面積按分法）

競争率を計算

category

e-stat：人口統計データ　

５― 国土数値情報　福祉施設データ　世界測地系　平成27年