8: 立地分析

公開

2024年6月5日

更新日

2024年8月7日

はじめに

今回から3回にわたって、これまでの講義内容の復習を兼ねたRを使った立地分析の演習を行います。今日の演習では、いろいろな空間データをつかって、佐賀県の待機児童について考えてみることにします。

待機児童とは、保育園等の利用を申し込んだ児童数から、実際に保育園等を利用している児童の数を除いた人数のことをいいます。つまり、保育園等を利用したいのに利用できずにいて、保育園等に空きが出るのを待っている（待機している）人数のことです。

日本の待機児童は、2017年の26,081人をピークに減少を続けており、2023年の待機児童は、ピークのおよそ10分の1となる2,680人でした。

待機児童の定義では、①育児休業中の者、②特定の保育園等のみ希望している者、③地方単独事業を利用してる者、④求職活動を休止している者、は待機児童に含めないことになっています（「除外4類型」と呼ばれます）。しかし、これらの人数についても、潜在的には待機児童であると考えることもできます。ちなみに、全国の除外4類型の人数は、2020年の74,840年をピークに減少傾向にありますが、2023年は66,168人でした。

準備

はじめに、今日の演習で利用するパッケージをロードしておきましょう。

library(tidyverse)
library(readxl)
library(sf)

佐賀県内自治体の待機児童

待機児童に関する全国の待機児童の傾向はわかりました。それでは、佐賀県の状況はどのようになっているでしょうか？

演習①待機児童数の塗り分け地図

佐賀県内自治体の待機児童数および潜在的待機児童数を調べ、塗り分け地図を作成しましょう。図を書くためにはどのようなデータが必要でしょうか。自治体のポリゴンデータと、待機児童数のデータが必要ですね。

地図データは、佐賀県オープンデータのポリゴン（市町村別）を使いましょう。 GeoJSONファイルをダウンロードし、プロジェクトのdataフォルダに移動してください。

sf::st_read関数で読み込みます。 sf::st_transform関数を使って、「平面直角座標系II系」に投影変換しておきます。

saga <- st_read("data/410004saga.geojson") |> 
  st_transform(6670)

データの内容を確認しておきましょう。市郡名と超村名が異なるデータ列に格納されていることがわかります。

head(saga, 5)

Simple feature collection with 5 features and 4 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: -117268.6 ymin: -1173.744 xmax: -54075.69 ymax: 69215.3
Projected CRS: JGD2011 / Japan Plane Rectangular CS II
  KEN_NAME GST_NAME CSS_NAME KEY_CODE                       geometry
1   佐賀県   鹿島市     <NA>   412074 MULTIPOLYGON (((-82226.11 3...
2   佐賀県   唐津市     <NA>   412023 MULTIPOLYGON (((-106350 451...
3   佐賀県   神埼市     <NA>   412104 MULTIPOLYGON (((-60455.97 2...
4   佐賀県   杵島郡   江北町   414247 MULTIPOLYGON (((-78005.98 2...
5   佐賀県   多久市     <NA>   412040 MULTIPOLYGON (((-76464.8 27...

後ほど、このポリゴンデータと待機児童のデータとを結合する必要があるのですが、結合する時のキーについて考えておかなくてはなりません。最も望ましいのは、結合するデータにもKEY_CODE列が存在して、それをキーにして結合することです。しかし、実は今回利用する待機児童データにはそのようなデータ列は存在しませんので、今回は市町名をキーにすることにします。

しかし先に見たように、このポリゴンデータは、市名と町名が異なる列に入っていますから、これを結合に使えるようするには、市名と町名をデータにもつ列を作成しなくてはいけません。ここでは、dplyr::if_else関数を使っています。新しいcity列を作成し、その内容は、CSS_NAME（町村名）がNA（すなわち、市）ならばGST_NAME（市郡名）を、そうでなければCSS_NAMEになるように指定しています。

saga <- saga |> 
  mutate(city = if_else(is.na(CSS_NAME), GST_NAME, CSS_NAME))
head(saga, 5)

Simple feature collection with 5 features and 5 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: -117268.6 ymin: -1173.744 xmax: -54075.69 ymax: 69215.3
Projected CRS: JGD2011 / Japan Plane Rectangular CS II
  KEN_NAME GST_NAME CSS_NAME KEY_CODE                       geometry   city
1   佐賀県   鹿島市     <NA>   412074 MULTIPOLYGON (((-82226.11 3... 鹿島市
2   佐賀県   唐津市     <NA>   412023 MULTIPOLYGON (((-106350 451... 唐津市
3   佐賀県   神埼市     <NA>   412104 MULTIPOLYGON (((-60455.97 2... 神埼市
4   佐賀県   杵島郡   江北町   414247 MULTIPOLYGON (((-78005.98 2... 江北町
5   佐賀県   多久市     <NA>   412040 MULTIPOLYGON (((-76464.8 27... 多久市

市町村別の待機児童数のデータは、こども家庭庁の保育所等関連状況取りまとめから入手できます。 Excelファイル「（参考）申込者の状況（令和5年4月1日）」をダウンロードし、プロジェクトのdataフォルダに移動してください。

readxl::read_excel関数で読み込みます。ここでは、Excelファイルの内容から、（最初に読み飛ばした2列を除いて）1列目に都道府県名、2列目に市町村名、3〜14列が数値データであることを指定しています。そして、dplyr::filter関数を使って、都道府県名（pref）が「佐賀県」であるデータ列だけを抽出しています。

data <- 
  read_excel('data/20230901_policies_hoiku_torimatome_r5_04.xlsx',
             skip = 7, 
             col_names = c('pref', 'city', str_c("X", 1:12)),
             col_types = c('skip', 'skip', 'text', 'text',
                           rep('numeric', 12))) |> 
  filter(pref == '佐賀県')
head(data, 5)

# A tibble: 5 × 14
  pref   city     X1    X2    X3    X4    X5    X6    X7    X8    X9   X10   X11
  <chr>  <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 佐賀県 佐賀…  6319  2986  2234   708   321     0     0     0     0    70     0
2 佐賀県 唐津…  3714  2496  1165     0    33     0     7     0     0    12     1
3 佐賀県 鳥栖…  1882  1524   211     0    51     0    46     0     0    50     0
4 佐賀県 多久…   494   308   183     1     2     0     0     0     0     0     0
5 佐賀県 伊万…  1680  1334   250     0    96     0     0     0     0     0     0
# ℹ 1 more variable: X12 <dbl>

以降の分析では、こども家庭庁のデータのうち、「申請者数」「潜在的待機児童数」「待機児童数」のデータのみ取り出しています。ここで「潜在的待機児童数」は、先に述べて「除外4要件」に当てはまる待機者数と「待機児童数」の合計にしています。

data <- data |> 
  mutate(potential = X8 + X9 + X10 + X11 + X12) |> 
  select(pref, city, applicant = X1, waiting = X12, potential)
head(data, 5)

# A tibble: 5 × 5
  pref   city     applicant waiting potential
  <chr>  <chr>        <dbl>   <dbl>     <dbl>
1 佐賀県 佐賀市        6319       0        70
2 佐賀県 唐津市        3714       0        13
3 佐賀県 鳥栖市        1882       0        50
4 佐賀県 多久市         494       0         0
5 佐賀県 伊万里市      1680       0         0

ポリゴンデータと統計データを結合する準備が整いました。 dplyr::left_join関数を使って結合します。

saga <- saga |> 
  left_join(data, by = join_by(city))

それでは、これをコロプレス図にしてみましょう。ここではggplot2::ggplotを使います。まずは、待機児童の塗り分け地図です。佐賀県東部（特に、上峰町・吉野ヶ里町・みやき町）の待機児童数が目立ちますね。

ggplot() + 
  geom_sf(aes(fill = waiting), data = saga) +
  scale_fill_distiller(direction = 1) +
  theme_minimal()

そして、潜在的待機児童の塗り分け地図です。先ほどの3町に加えて、佐賀市、鳥栖市、唐津市などで潜在的待機児童数が大きいことがわかります。

ggplot() + 
  geom_sf(aes(fill = potential), data = saga) +
  scale_fill_distiller(direction = 1) +
  theme_minimal()

さて、これで演習①の当初の目的は達成したのですが、せっかくデータを準備したので、もう少し掘り下げてみましょう。まず、塗り分け地図について、待機児童数そのままのデータではなく、申請者に占める待機児童の比率に直したら、結果はどのようになるでしょうか。

saga <- saga |> 
  mutate(rate = waiting / applicant * 100)
ggplot() + 
  geom_sf(aes(fill = rate), data = saga) +
  scale_fill_distiller(direction = 1) +
  theme_minimal()

そして申請者に占める潜在的待機児童の比率です。

saga <- saga |> 
  mutate(rate = potential / applicant * 100)
ggplot() + 
  geom_sf(aes(fill = rate), data = saga) +
  scale_fill_distiller(direction = 1) +
  theme_minimal()

やはり佐賀県の待機児童については、東高西低の傾向があるようです。佐賀県の平均年齢が西高東低であったことを考えると、平均年齢の低い県東部で待機児童が相対的に多い、というのはそれと整合的な結果だといえるかもしれません。

つまり、そもそも県西部では保育所に対する需要が少ないのではないか、という仮説です。そこで最後に、自治体別の0〜5最人口に占める保育所等申請者の比率を可視化してみましょう。

人口データは、佐賀県統計年鑑のデータを使ってみます。 Excelファイル「第4章　人口及び世帯」をダウンロードし、プロジェクトのdataフォルダに移動してください。

readxl::read_excel関数で読み込みましょう。保育所を利用するのは、6歳未満の乳幼児だと考えられます。本来であれば、0〜5歳の人口が分かればいいのですが、残念ながら5歳階級の人口データしか公表されていないようなので、ここでは、0〜4歳人口に、5〜9歳人口の5分の1を加えることで、0〜5歳人口を代替しています。

pop <- 
  read_excel('data/3_101414_309615_up_vpymzuh2.xlsx',
             sheet = '4-7(1)', skip = 13, 
             col_names = c('city', str_c('X', 1:17)),
             col_types = c('skip', 'text', rep('numeric', 17), 'skip')) |> 
  drop_na() |> 
  mutate(toddler = X4 + X6 / 5) |> 
  select(city, toddler)
head(pop, 5)

# A tibble: 5 × 2
  city     toddler
  <chr>      <dbl>
1 佐賀市    10617.
2 唐津市     5352.
3 鳥栖市     3832.
4 多久市      705 
5 伊万里市   2402.

dplyr::left_joinで結合しましょう。

saga <- saga |> 
  left_join(pop, by = join_by(city))

結合したら、0〜5歳人口あたりの保育所等申請者数を塗り分け地図にしてみます。

saga <- saga |> 
  mutate(rate = applicant / toddler * 100)
ggplot() + 
  geom_sf(aes(fill = rate), data = saga) +
  scale_fill_distiller(direction = 1) +
  theme_minimal()

意外にも、玄海町や太良町の申請率が高く、鳥栖市が最も低いという結果になりました。（なぜこのような結果になったのか、思いつく仮説の1つは幼稚園の存在です。佐賀市や鳥栖市などの都市部では、保育所等ではなく、幼稚園に通う子供が多いのかもしれません。他にどのような仮説があるか、皆さんも考えてみてください）。

さて、佐賀県内の潜在的待機児童の大半を占めるのは、「特定の保育園等のみ希望している」児童です。これは、他の（通える範囲と見なされる¹）保育園等には空きがあるのだけれど、希望する保育園等には空きがないために、通園できないでいる児童です。保育サービスの需要と供給のミスマッチが起こっていることが問題だといえます。そこで次に、佐賀市を事例として、空間的な需要と供給のミスマッチを可視化しましょう。具体的には、保育所等の立地と6歳未満人口の分布を可視化し、観察します。

演習②佐賀市の保育所立地と6歳未満人口分布

保育所のデータは、国土数値情報にある「福祉施設（ポイント）」を使いましょう。

佐賀県の令和3年のデータ（P14-21_41_GML.zip）をダウンロードしてください。 zipファイルを解凍し、できたP14-21_41_GMLフォルダをプロジェクトのdataフォルダに移動して下さい。

このファイルには、あらゆる福祉施設のデータが入っていますので、この中から保育所のデータだけを抜き出す必要があります。データには属性データとして、福祉施設大分類、福祉施設中分類、福祉施設小分類が入っていいますので、今回はこれを使って抜き出します。

福祉施設中分類コードを見ると、0504：保育所等、0505：地域型保育事業所、0506：認可外保育施設を抜き出せば良さそうです。属性データを使って一部のデータだけを抜き出すには、dplyr::filter関数を使います。福祉施設中分類（P14_006）が”0504”、“0505”、“0506”のいずれかであるデータだけを抽出しましょう。

nursery <- st_read("data/P14-21_41_GML/P14-21_41.geojson") |> 
  filter(P14_006 %in% c("0504", "0505", "0506")) |> 
  st_transform(6670)

Reading layer `P14-21_41' from data source 
  `/Users/kzktmr/Documents/LectureNotes/data/P14-21_41_GML/P14-21_41.geojson' 
  using driver `GeoJSON'
Simple feature collection with 2520 features and 10 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: 129.7636 ymin: 32.96552 xmax: 130.537 ymax: 33.5962
Geodetic CRS:  JGD2011

head(nursery, 5)

Simple feature collection with 5 features and 10 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: -94530.19 ymin: 11738.72 xmax: -53796.57 ymax: 35596.83
Projected CRS: JGD2011 / Japan Plane Rectangular CS II
  P14_001          P14_002 P14_003              P14_004 P14_005 P14_006 P14_007
1  佐賀県           小城市   41208        牛津町牛津664      05    0504  050402
2  佐賀県           嬉野市   41209 嬉野町大字下宿乙1279      05    0505  050501
3  佐賀県 神埼郡吉野ヶ里町   41327           字立野1077      05    0504  050402
4  佐賀県   三養基郡上峰町   41345          大字坊所710      05    0504  050402
5  佐賀県   三養基郡上峰町   41345          大字坊所699      05    0504  050402
                     P14_008 P14_009 P14_010                   geometry
1       牛津ルーテルこども園       5       1  POINT (-74531.6 27597.24)
2       うれしのつぼみ保育園       4       2 POINT (-94530.19 11738.72)
3       認定こども園　きらり       5       1 POINT (-55112.59 35344.28)
4 認定こども園かみみね幼稚園       5       1 POINT (-53796.57 35596.83)
5             ひかりこども園       5       1    POINT (-53804 35455.55)

佐賀県にある449の保育所データが登録されていることがわかります。

これを佐賀県の地図と合わせて表示てみましょう。佐賀県の地図データは、以前の講義で使った410004saga.geojsonを使いましょう（以前のプロジェクトフォルダからコピーしてください）。

自治体コードで色をつけて、地図に重ねて表示してみます。

ggplot() +
  geom_sf(aes(fill = KEY_CODE), data = saga, alpha = 0.5) +
  geom_sf(aes(color = P14_003), data = nursery) +
  theme_minimal()

市町によって保育所の数に差がありそうですね。

では、佐賀県には20市町ありますが、それぞれの市町にいくつの保育所があるのでしょうか。 dplyr::group_byとdplyr::summariseを使って数えてみましょう。

このような時には、dplyr::count関数が便利です。

nursery |>
  count(P14_002) |> 
  arrange(desc(n))

Simple feature collection with 20 features and 2 fields
Geometry type: MULTIPOINT
Dimension:     XY
Bounding box:  xmin: -114729.4 ymin: -3252.091 xmax: -43067.99 ymax: 66661.28
Projected CRS: JGD2011 / Japan Plane Rectangular CS II
First 10 features:
        P14_002   n                       geometry
1        佐賀市 138 MULTIPOINT ((-74988.15 4660...
2        唐津市  60 MULTIPOINT ((-114729.4 6428...
3        鳥栖市  43 MULTIPOINT ((-48606.03 4121...
4      伊万里市  38 MULTIPOINT ((-112373.3 3759...
5        武雄市  24 MULTIPOINT ((-99392.73 2193...
6        小城市  19 MULTIPOINT ((-76407.43 2639...
7        鹿島市  18 MULTIPOINT ((-87182.73 9215...
8        嬉野市  17 MULTIPOINT ((-95470.91 1170...
9        多久市  14 MULTIPOINT ((-87852.49 2976...
10 杵島郡白石町  14 MULTIPOINT ((-83397.25 2102...

count関数は、引数に指定したグループごとにデータの件数を数え、nという列に件数のデータを格納します²。ここでは、自治体名を使ってグルーピングしています。これで自治体ごとの保育所数がわかりますので、dplyr::arrange関数とdplyr::desc関数を使って、データを保育所数の降順に並び替えています。

これを見ると、最も多い佐賀市には138の保育所があることがわかります。最も少ないのは玄海町で、町内の保育所数は2です。

以降の分析では、佐賀市の保育所データだけを使いますので、dplyr::filter関数を使って、データを絞り込んでおきましょう。

nursery <- nursery |> filter(P14_002 == "佐賀市")
head(nursery, 5)

Simple feature collection with 5 features and 10 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: -66722.7 ymin: 22391.03 xmax: -61945.25 ymax: 35849.97
Projected CRS: JGD2011 / Japan Plane Rectangular CS II
  P14_001 P14_002 P14_003               P14_004 P14_005 P14_006 P14_007
1  佐賀県  佐賀市   41201            金立町2467      05    0506  050600
2  佐賀県  佐賀市   41201           若宮1-13-17      05    0504  050403
3  佐賀県  佐賀市   41201 川副町大字早津江560-1      05    0504  050402
4  佐賀県  佐賀市   41201       兵庫町大字渕892      05    0504  050402
5  佐賀県  佐賀市   41201            長瀬町2-18      05    0504  050402
                                     P14_008 P14_009 P14_010
1                                 金立保育園       5       1
2                   西九州大学附属三光保育園       5       1
3                           博愛の里こども園       5       1
4                           エミールこども園       5       1
5 日新こども園(日新こども園好生館分園きらら)       5       1
                    geometry
1 POINT (-65032.81 35849.97)
2 POINT (-65992.32 30151.54)
3 POINT (-61945.25 22391.03)
4 POINT (-63537.64 29852.58)
5  POINT (-66722.7 27928.32)

小地域の年齢別人口のデータでは、5歳階級のものが入手できる最も細かいデータですので、先ほどと同様に、0〜4歳の人口と、5〜9歳の人口の5分の1を足し合わせたものを、保育所の潜在的な需要とみなすことにしましょう。

使うデータは、2020年国勢調査の5歳階級人口データです、e-Statからダウンロードしてください。境界データは、第5回の講義で利用したものと同じですので、第5回のプロジェクトフォルダから、「A002005212020DDSWC41201-JGD2011」フォルダをコピーして持ってきても構いません。統計データは、前回の講義で使用したものと同じですので、前回のプロジェクトフォルダから、「tblT001082C41.txt」ファイルをコピーして使っても構いません。

境界データ
- [小地域]→[国勢調査]→[2020年]→[小地域（町丁・字等）（JGD2011）]→[世界測地系緯度経度・Shapefile]→[41 佐賀県]→[41201 佐賀市]からダウンロード
- 境界データ（A002005212020DDSWC41201-JGD2011.zip）を展開してできたフォルダを、プロジェクトのdataフォルダに移動
統計データ
- [国勢調査]→[2020年]→[小地域（町丁・字等）]→[年齢（５歳階級、４区分）別、男女別人口]→[41 佐賀県]からCSVファイルをダウンロード
- 統計データ（tblT001082C41.zip）を展開してできたファイル（tblT001082C41.txt）を、プロジェクトのdataフォルダに移動

まずは境界データをsf::st_read関数を使って読み込みます。

district <- 
  st_read("data/A002005212020DDSWC41201-JGD2011/r2ka41201.shp") |> 
  st_transform(6670)

Reading layer `r2ka41201' from data source 
  `/Users/kzktmr/Documents/LectureNotes/data/A002005212020DDSWC41201-JGD2011/r2ka41201.shp' 
  using driver `ESRI Shapefile'
Simple feature collection with 484 features and 29 fields
Geometry type: POLYGON
Dimension:     XY
Bounding box:  xmin: 130.1393 ymin: 33.14109 xmax: 130.3794 ymax: 33.48151
Geodetic CRS:  JGD2011

次に、統計データをreadr::read_csv関数を使って読み込みます。前回と同様、ファイルをヘッダ部分とデータ部分に分けて、read_csv関数を2回使って読み込みます。

まずは、先頭行だけ読み込んで、列名の文字列ベクトルを作ります。

col_names <- 
  read_csv("data/tblT001082C41.txt", col_types = "c", 
           col_names = FALSE, n_max = 1, 
           locale = locale(encoding = "cp932")) |> 
  as.character()
col_names

 [1] "KEY_CODE"   "HYOSYO"     "CITYNAME"   "NAME"       "HTKSYORI"  
 [6] "HTKSAKI"    "GASSAN"     "T001082001" "T001082002" "T001082003"
[11] "T001082004" "T001082005" "T001082006" "T001082007" "T001082008"
[16] "T001082009" "T001082010" "T001082011" "T001082012" "T001082013"
[21] "T001082014" "T001082015" "T001082016" "T001082017" "T001082018"
[26] "T001082019" "T001082020" "T001082021" "T001082022" "T001082023"
[31] "T001082024" "T001082025" "T001082026" "T001082027" "T001082028"
[36] "T001082029" "T001082030" "T001082031" "T001082032" "T001082033"
[41] "T001082034" "T001082035" "T001082036" "T001082037" "T001082038"
[46] "T001082039" "T001082040" "T001082041" "T001082042" "T001082043"
[51] "T001082044" "T001082045" "T001082046" "T001082047" "T001082048"
[56] "T001082049" "T001082050" "T001082051" "T001082052" "T001082053"
[61] "T001082054" "T001082055" "T001082056" "T001082057" "T001082058"
[66] "T001082059" "T001082060"

次に、列名の文字列ベクトルを利用して、CSVファイルを読み込みます。ここでは、空白「」と欠損値「-」および秘匿データ（X）をNAにするように設定しています。また、文字コードを指定して、日本語が正しく読み込まれるようにしておきます。

population <- 
  read_csv("data/tblT001082C41.txt", skip = 2,
           na = c("", "-", "X"), col_names = col_names,
           locale = locale(encoding = "cp932"))

Rows: 2567 Columns: 67
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr  (4): CITYNAME, NAME, HTKSAKI, GASSAN
dbl (63): KEY_CODE, HYOSYO, HTKSYORI, T001082001, T001082002, T001082003, T0...

ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

head(population, 5)

# A tibble: 5 × 67
    KEY_CODE HYOSYO CITYNAME NAME  HTKSYORI HTKSAKI GASSAN T001082001 T001082002
       <dbl>  <dbl> <chr>    <chr>    <dbl> <chr>   <chr>       <dbl>      <dbl>
1    4.12e 4      1 佐賀市   <NA>         0 <NA>    <NA>       233301       9001
2    4.12e 8      3 佐賀市   駅前…        0 <NA>    <NA>         2496        110
3    4.12e10      4 佐賀市   駅前…        0 <NA>    <NA>          583         12
4    4.12e10      4 佐賀市   駅前…        0 <NA>    <NA>         1081         53
5    4.12e10      4 佐賀市   駅前…        0 <NA>    <NA>          832         45
# ℹ 58 more variables: T001082003 <dbl>, T001082004 <dbl>, T001082005 <dbl>,
#   T001082006 <dbl>, T001082007 <dbl>, T001082008 <dbl>, T001082009 <dbl>,
#   T001082010 <dbl>, T001082011 <dbl>, T001082012 <dbl>, T001082013 <dbl>,
#   T001082014 <dbl>, T001082015 <dbl>, T001082016 <dbl>, T001082017 <dbl>,
#   T001082018 <dbl>, T001082019 <dbl>, T001082020 <dbl>, T001082021 <dbl>,
#   T001082022 <dbl>, T001082023 <dbl>, T001082024 <dbl>, T001082025 <dbl>,
#   T001082026 <dbl>, T001082027 <dbl>, T001082028 <dbl>, T001082029 <dbl>, …

たくさんのデータがありますが、このうち0〜4歳人口と5〜9歳人口のデータを使って、0〜5歳人口を推計しています。また同時に、列の名前をわかりやすく「toddler」にしておきます（同時に、あとの演習で利用するために、65歳以上人口の列も「elderly」という列名で取り出しておきます）。また、KEY_CODEを文字列型に変更しているのは、後でleft_joinのキーにする際に、saga_cityのKEY_CODEとデータ型が違うとエラーになるので、その対策です。

population <- population |> 
  filter(CITYNAME == "佐賀市") |> 
  mutate(KEY_CODE = as.character(KEY_CODE),
         toddler = T001082002 + T001082003 / 5,
         elderly = T001082019) |> 
  select(KEY_CODE, toddler, elderly)
head(population, 5)

# A tibble: 5 × 3
  KEY_CODE    toddler elderly
  <chr>         <dbl>   <dbl>
1 41201       11085.    64802
2 412010010     128.      349
3 41201001001    15.4      73
4 41201001002    61.2     132
5 41201001003    51.8     144

KEY_CODEをキーにして、境界データに統計データを結合します。 dplyr::left_join関数を使います。また、キーを指定するのに、dplyr::join_by関数を使っています。

district <- district |> 
  left_join(population, by = join_by(KEY_CODE))

0〜5歳人口と保育所の分布を重ねて表示してみましょう。

ggplot() + 
  geom_sf(aes(fill = toddler), data = district) +
  geom_sf(data = nursery, alpha = 0.5) +
  scale_fill_distiller(direction = 1) +
  theme_minimal()

さて、ここから、空間的な需要と供給のミスマッチを議論したいのですが、需要側のデータである保育園等利用申請者数の小地域データも、供給側である保育所等の定員のデータもありませんので、ここでは簡易的な分析になってしまいますが、保育所等と居住地との距離を用いた分析をしてみましょう。

まず、居住地から距離1,000 m以内に保育所が1つもない地区を可視化してみましょう。 sf::st_nearest_feature関数を使って、それぞれの小地域から最も近い保育所を特定しています。そして、sf::st_distance関数を使って、それぞれの小地域とそこから最も近い保育所等の距離を計算しています。最後に、nursery1というデータ列を作成しています。最寄りのの保育所までの距離が1,000 mを超えるかどうかという論理型データ（TRUEまたはFALSEのどちらかの値を持つ）になっています。

nearest_nursery <-
  st_nearest_feature(district, nursery)
nearest_nursery_dist <-
  st_distance(district, nursery[nearest_nursery,], by_element = TRUE) |> 
  units::drop_units()
district <- district |> 
  mutate(nursery1 = nearest_nursery_dist > 1000)

ggplot() +
  geom_sf(aes(fill = nursery1), data = district) +
  theme_minimal()

また、市全体の0〜6歳人口のうち、最寄り保育所までの距離が1,000 mを超える地域に住む人口の比率は、以下のようにして計算できます。

sum(district$toddler * district$nursery1, na.rm = TRUE) / sum(district$toddler, na.rm = TRUE) * 100

[1] 3.585765

さて、「潜在的待機児童」の大半は、通園可能と考えられる範囲に保育所等があるにも関わらず、何らかの理由によって待機しているでした。このような「潜在的待機児童」を減らすためには、何が必要となるのでしょうか。

最も必要となることの1つは、保育所等の選択肢を増やすことでしょう。通園可能な保育所が0であるよりも1つでもあったほうがいいですし、通園可能な保育所等が1つしかないよりも2つ以上あったほうがより良いといえます。

そこで次に、居住地から、距離1,000 m以内の保育所が2箇所未満である（1つもない、もしくは1つしかない）地区を可視化してみましょう。これを実現する方法は、いくつか考えられますが、ここでは、地区と保育所の距離行列を計算した後、地区ごとに、そこから距離1,000 mの範囲にある保育所等の数を数えることにします。

距離行列は、sf::st_distance関数を使って計算できます。これを、距離が1,000 m以下かどうかを表す論理行列に変換し、rowSums関数でその行方向の和を計算しています（論理型の値は、TRUE=1、FALSE=0として計算できますので、行和を計算することで、TRUEの個数をカウントすることができます）。最後に、この行和（つまり、距離1,000 m以下にある保育所等の数）が1以下かどうかという論理型データをnursery2として新しいデータ列にして追加しています。

distance_matrix <- st_distance(district, nursery) |> 
  units::drop_units()
n_nursery <- rowSums(distance_matrix <= 1000)
district <- district |> 
  mutate(nursery2 = n_nursery <= 1)

結果を塗り分け地図にしてみましょう。

ggplot() +
  geom_sf(aes(fill = nursery2), data = district) +
  theme_minimal()

先ほどと比べて、TRUEである地区、すなわち保育所までのアクセス性が悪いと評価される地区が増えていることがわかると思います。 0〜5最人口に占める比率も計算してみましょう。

sum(district$toddler * district$nursery2, na.rm = TRUE) / sum(district$toddler, na.rm = TRUE) * 100

[1] 7.064822

先ほどの最も近い保育所までの距離を考えた時よりも、比率が大きくなっていることが確認できました。

このように、\(n\)番目に近い施設までの距離を考えることは重要です。特に、近年多発する自然災害など、不確実性を考慮した際には、最寄りの施設が利用できなくなった時のバックアップをどうするかを考えておく必要があります。「潜在的待機児童」の、何らかの理由によって最寄りの保育所等に通園できないという問題も、同じ枠組みで考えることができると思います。

佐賀市のバスアクセス

佐賀市内には、路線バス事業者が4社バスを走らせていますし、市北部ではコミュニティバスも走っています。佐賀駅バスセンターはとてもきれいて使いやすいですし、バスセンターでの他社バス間の乗り継ぎ割引という制度は、全国的にみてもあまり例はないのではないかと思います。

私も毎週、佐賀駅バスセンターから佐賀大学までバスに乗ってきますが、バスはいつも満席です。学生と高齢者の利用が多いと感じます。

そこで、次のテーマとして、佐賀市内のバスアクセスを取り上げることにします。

演習③佐賀市内のバス停分布と高齢人口分布

バス停のデータは、国土数値情報からダウンロードしてください。バス停（ポイント）データと、バスルート（ライン）データが別々のファイルになっています。

バス停留所データから、佐賀県の2022年度のShapeFile形式・GeoJSON形式のデータ（P11-22_41_SHP.zip）をダウンロードします。
バスルートデータから、佐賀県の2022年度のShapeFile形式・GeoJSON形式のデータ（N07-22_41_SHP.zip）をダウンロードします。

ダウンロードしたら、ダウンロードフォルダで展開し、できたフォルダごと、プロジェクトのdataフォルダに移動してください。

sf::st_read関数で読み込みます。平面直角座標系II系に投影変換しておきます。

stops <- st_read('data/P11-22_41_SHP/P11-22_41.geojson') |> 
  st_transform(6670)

Reading layer `P11-22_41' from data source 
  `/Users/kzktmr/Documents/LectureNotes/data/P11-22_41_SHP/P11-22_41.geojson' 
  using driver `GeoJSON'
Simple feature collection with 3080 features and 73 fields
Geometry type: POINT
Dimension:     XY
Bounding box:  xmin: 129.7615 ymin: 32.95459 xmax: 130.5412 ymax: 33.55822
Geodetic CRS:  JGD2011

route <- st_read('data/N07-22_41_SHP/N07-22_41.geojson') |> 
  st_transform(6670)

Reading layer `N07-22_41' from data source 
  `/Users/kzktmr/Documents/LectureNotes/data/N07-22_41_SHP/N07-22_41.geojson' 
  using driver `GeoJSON'
Simple feature collection with 4021 features and 2 fields
Geometry type: LINESTRING
Dimension:     XY
Bounding box:  xmin: 129.7592 ymin: 32.95459 xmax: 130.5418 ymax: 33.55822
Geodetic CRS:  JGD2011

ggplotを使って、佐賀市内のバス路線図を描いてみましょう。 sf::st_filter関数を使って、佐賀県のバス停とバスルートデータから、佐賀市内にあるものだけを抽出します。 ggplotを使って、地図にします。「バス事業者名」で色を変えるように設定しています。

saga_city <- saga |> filter(GST_NAME == '佐賀市')
saga_stops <- stops |> st_filter(saga_city)
saga_route <- route |> st_filter(saga_city)
ggplot() +
  geom_sf(data = saga_city) +
  geom_sf(aes(color = P11_002), data = saga_stops, alpha = 0.5) +
  geom_sf(aes(color = N07_001), data = saga_route) +
  theme_minimal()

神埼市のコミュニティバスのバス停が、一部佐賀市内にあるようですね。神埼市のバスを除いた、それ以外のバス停を地図にしてみます。このデータを分析に使うことにしましょう。

saga_stops <- saga_stops |> 
  filter(P11_002 != '神埼市')
ggplot() +
  geom_sf(data = saga_city) +
  geom_sf(aes(color = P11_002), data = saga_stops, alpha = 0.5) +
  theme_minimal()

次は高齢人口分布です。これは、演習②で使った国勢調査の小地域人口データを使いましょう。先ほど作ったdistrictデータに、65歳以上人口のelderly列があります。 ggplotを使って、65歳以上人口の塗り分け地図を描いてみます。

ggplot() +
  geom_sf(aes(fill = elderly), data = district) +
  scale_fill_distiller(direction = 1) +
  theme_minimal()

これと、先ほどのバス停データを重ねてみましょう。

ggplot() +
  geom_sf(aes(fill = elderly), data = district) +
  geom_sf(data = saga_stops, color = 'blue', alpha = 0.5) +
  scale_fill_distiller(direction = 1) +
  theme_minimal()

さて、このバス停の位置と高齢者の分布を分析したいわけですが、先ほどの保育所の分析とは違うツールを使いましょう。ここでは、バッファーを使ってバス停からの距離と人口分布との関係を分析します。 sf::st_buffer関数を使って、バス停から500 mのバッファーを作成しています。各バス停から500 mの円形領域が作成されますので、これにsf::st_union関数を適用して、1つの領域に集計しています。

buffer <- saga_stops |> st_buffer(dist = 500) |> st_union()
ggplot() +
  geom_sf(aes(fill = elderly), data = district) +
  geom_sf(data = buffer, color = 'red', alpha = 0.5) +
  scale_fill_distiller(direction = 1) +
  theme_minimal()

このバッファー領域に含まれる高齢者人口を集計しましょう。色々な考え方があると思いますが、ここでは面積按分によって各小地域の人口を、バッファーに含まれる部分と含まれない部分とに分けてみましょう。

はじめに、小地域の面積を計算しています。国勢調査の境界データには、面積データ（AREA）が入っているのですが、ここでは面積按分する際に結果を整合させるために、同じ方法で計測した面積を使いたいために、あえて面積を計算しています。

sf::sf_intersection関数を使って、小地域領域とバス停のバッファー領域の共通部分を取り出しています。含まれる部分を取り出しています。

district <- district |> 
  mutate(area = st_area(district))
intersect <- district |> 
  st_intersection(buffer)

Warning: attribute variables are assumed to be spatially constant throughout
all geometries

ggplot() +
  geom_sf(data = saga_city) +
  geom_sf(aes(fill = elderly), data = intersect) +
  scale_fill_distiller(direction = 1) +
  theme_minimal()

小地域領域のうちバッファーに含まれる部分が抽出できるので、この面積を計算することで小地域のなかのバッファー部分の面積比率が計算できることになります。小地域全体の高齢者人口にこの面積比率を掛けることで、小地域のうちバス停バッファー部分に相当する人口を推計します。つまり、面積按分によって対象地域の人口を推計しています。

intersect <- intersect |> 
  mutate(new_area = st_area(intersect),
         ratio = new_area / area,
         new_elderly = elderly * ratio)
head(intersect, 5)

Simple feature collection with 5 features and 37 fields
Geometry type: POLYGON
Dimension:     XY
Bounding box:  xmin: -65669.76 ymin: 28615.66 xmax: -64125.32 ymax: 30085.38
Projected CRS: JGD2011 / Japan Plane Rectangular CS II
     KEY_CODE PREF CITY S_AREA PREF_NAME CITY_NAME         S_NAME KIGO_E HCODE
1 41201001001   41  201 001001    佐賀県    佐賀市 駅前中央一丁目   <NA>  8101
2 41201001002   41  201 001002    佐賀県    佐賀市 駅前中央二丁目   <NA>  8101
3 41201001003   41  201 001003    佐賀県    佐賀市 駅前中央三丁目   <NA>  8101
4 41201002001   41  201 002001    佐賀県    佐賀市     大財一丁目   <NA>  8101
5 41201002002   41  201 002002    佐賀県    佐賀市     大財二丁目   <NA>  8101
       AREA PERIMETER R2KAxx R2KAxx_ID KIHON1 DUMMY1 KIHON2  KEYCODE1  KEYCODE2
1 147281.68  1528.456      1         0   0010      -     01 201001001 201001001
2 108748.65  1371.930      2         1   0010      -     02 201001002 201001002
3 144423.54  1875.425      3         2   0010      -     03 201001003 201001003
4  74548.31  1150.161      4         3   0020      -     01 201002001 201002001
5 137859.16  1900.060      5         4   0020      -     02 201002002 201002002
  AREA_MAX_F KIGO_D N_KEN N_CITY KIGO_I KBSUM JINKO SETAI   X_CODE   Y_CODE
1          M   <NA>  <NA>   <NA>   <NA>    16   583   378 130.2978 33.26336
2          M   <NA>  <NA>   <NA>   <NA>    14  1081   576 130.2995 33.26637
3          M   <NA>  <NA>   <NA>   <NA>    19   832   432 130.3041 33.26694
4          M   <NA>  <NA>   <NA>   <NA>    11   337   202 130.3027 33.25728
5          M   <NA>  <NA>   <NA>   <NA>    10   797   413 130.3083 33.25828
   KCODE1 toddler elderly nursery1 nursery2            area
1 0010-01    15.4      73    FALSE    FALSE 147286.14 [m^2]
2 0010-02    61.2     132    FALSE    FALSE 108751.81 [m^2]
3 0010-03    51.8     144    FALSE    FALSE 144427.64 [m^2]
4 0020-01    12.8      93    FALSE    FALSE  74550.14 [m^2]
5 0020-02    35.6     200    FALSE    FALSE 137862.50 [m^2]
                        geometry        new_area         ratio  new_elderly
1 POLYGON ((-65198.19 29415.8... 147286.14 [m^2] 1.0000000 [1]  73.0000 [1]
2 POLYGON ((-65456.69 29590.8... 108751.81 [m^2] 1.0000000 [1] 132.0000 [1]
3 POLYGON ((-64843.62 29876.6... 144427.64 [m^2] 1.0000000 [1] 144.0000 [1]
4 POLYGON ((-64772.49 28755.2...  74550.14 [m^2] 1.0000000 [1]  93.0000 [1]
5 POLYGON ((-64354.59 28748.7... 135304.75 [m^2] 0.9814471 [1] 196.2894 [1]

面積按分によって求めた小地域の高齢人口を、全地域について足し上げることで集計しましょう。

city_elderly <- district |> pull(elderly) |> sum(na.rm = TRUE) 
bus_elderly <- intersect |> pull(new_elderly) |> sum(na.rm = TRUE)
city_elderly

[1] 68281

bus_elderly

52375.62 [1]

bus_elderly / city_elderly * 100

76.706 [1]

佐賀市の高齢人口は68,281ですが、そのうちバス停から500 mの範囲の人口は52,376人と推計されました。高齢人口の約76.7％（およそ4分の3）はバスが利用しやすい場所に住んでいると言えそうです。

この数字は高いと思いますか、それとも低いと思いますか？

脚注

「他に利用可能な特定教育・保育施設又は特定地域型保育事業等があるにも関わらず、特定の保育所等を希望し、保護者の私的な理由により待機している場合には待機児童数には含めない」とされています。ここで「他に利用可能な特定教育・保育施設又は特定地域型保育事業等」とは、「開所時間が保護者の需要に応えている（例えば、希望の保育所と開所時間に差異がない）、立地条件が登園するのに無理がない（例えば、通常の交通手段により、自宅から20〜30分未満で登園が可能）」などとされています。↩︎
つまり、nursery |> group_by(P14_002) |> summarise(n = n())とほぼ同じ意味です。↩︎