ステップ0：パッケージの読み込み

library(sf)
library(tidyverse)
library(ggplot2)
library(openxlsx)
library(ggspatial)

ステップ1：e-Statの500mメッシュ統計・境界データのダウンロード

省略．

ステップ2：統計データを境界データに結合できる形式に加工

省略．

ステップ3：神奈川県三浦市の500mメッシュ境界データの作成

省略．

ステップ4：統計データを境界データに結合

まずは，全ての統計データ・境界データを読み込む．

#神奈川県三浦市周辺のバス停留所データ（bus_stops_miura.shp）
bus_stops_miura <- sf::st_read(dsn="bus_stops_miura.shp",
                      stringsAsFactors=FALSE,
                      quiet=TRUE,
                      options="ENCODING=UTF-8")
#神奈川県三浦市周辺のバスルートデータ（bus_routes_miura.shp）
bus_routes_miura <- sf::st_read(dsn="bus_routes_miura.shp",
                      stringsAsFactors=FALSE,
                      quiet=TRUE,
                      options="ENCODING=UTF-8") %>%
  #座標精度に関する属性は今回使わないので削除
  sf::st_zm()
#神奈川県三浦市の国勢調査小地域境界データ（census_sa_miura.shp）
#文字コードをUTF-8で，shapefile内部の変数をfactor化せずに読み込む
census_sa_miura <- sf::st_read(dsn="census_sa_miura.shp",
                               stringsAsFactors=FALSE,
                               quiet=TRUE,
                               options="ENCODING=UTF-8")
#国勢調査データ（elderly_pop.xlsx）
elderly_pop <- openxlsx::read.xlsx(xlsxFile="elderly_pop.xlsx",sheet="elderly_pop")
#神奈川県三浦市周辺の医療機関データ（med_miura.shp）
#文字コードをUTF-8で，shapefile内部の変数をfactor化せずに読み込む
med_miura <- sf::st_read(dsn="med_miura.shp",
                         stringsAsFactors=FALSE,
                         quiet=TRUE,
                         options="ENCODING=UTF-8")
#神奈川県三浦市に重なる4次メッシュ（mesh500m_miura.shp）
#shapefile内部の変数をfactor化せずに読み込む
mesh500m_miura <- sf::st_read(dsn="mesh500m_miura.shp",
                      stringsAsFactors=FALSE,
                      quiet=TRUE)

次に，結合基準の変数にKEY_CODEを用いて境界データへ統計データを結合．

mesh500m_miura <- mesh500m_miura %>%
  #境界データへ統計データを結合
  dplyr::left_join(y=elderly_pop,by="KEY_CODE") %>%
  #統計データの変数「６５歳以上人口総数」結合されたレコードのみ残す
  dplyr::filter(!is.na(６５歳以上人口総数))

ステップ5：高齢者人口分布図の作成

結合された統計データ（６５歳以上人口総数）の値を地図上にプロットする．

#６５歳以上人口総数を分級する階級を設定
#６５歳以上人口総数の最小値・最大値を取得
min(mesh500m_miura$６５歳以上人口総数)
## [1] 1
max(mesh500m_miura$６５歳以上人口総数)
## [1] 970
#下限値を６５歳以上人口総数の最小値，上限値を６５歳以上人口総数の最大値とし，100刻みで分級
range <- c(1,100,200,300,400,500,600,700,800,900,970)

ggplot2::ggplot() +
  #最初に小地域境界を白塗りでプロット
  ggplot2::geom_sf(data=census_sa_miura,fill="white") +
  #６５歳以上人口総数が結合された500mメッシュをプロット
  ggplot2::geom_sf(data=mesh500m_miura,
                   #上で作成した分級
                   aes(fill=cut(６５歳以上人口総数,breaks=range)),
                   #不透明度80%で塗り分け
                   alpha=0.8) +
  #階級の色分けには「Reds」というパレットを使う
  ggplot2::scale_fill_brewer(palette="Reds") +
  #図の凡例のタイトルを手動で指定する
  ggplot2::guides(fill=guide_legend(title="６５歳以上人口総数")) +
  #日本語が文字化けしない為の設定（Windowsで実行の際は不要）
  ggplot2::theme_gray(base_family="HiraKakuPro-W3")

ステップ6：国土数値情報の医療機関データのダウンロードと加工

ggplotでは，条件式を組み込んだ可視化が可能．医療機関データに対して適用してプロットしてみる．

#med_miura内で，医療機関分類コードP04_001の値に対応したラベル変数を作成
med_miura <- med_miura %>%
  dplyr::mutate(class=dplyr::case_when(
    P04_001==1 ~ "病院",
    P04_001==2 ~ "診療所",
    P04_001==3 ~ "歯科診療所"
  ))
  
ggplot2::ggplot() +
  #色分けとポイントの大きさは下の行で設定したいので，aes内で色・サイズ分けの基準になる変数名（変数class）のみ指定
  ggplot2::geom_sf(data=med_miura,aes(color=class,size=class)) +
  #変数classに基づき，病院なら青，診療所なら緑，歯科診療所なら赤でプロット
  scale_color_manual(name="class",
                     breaks=c("病院","診療所","歯科診療所"),
                     values=c("blue","green","red")) +
  #変数classに基づき，病院ならサイズ3，診療所なら2，歯科診療所なら1でプロット
  scale_size_manual(name="class",
                    breaks=c("病院","診療所","歯科診療所"),
                    values=c(3,2,1)) +
  #日本語が文字化けしない為の設定（Windowsで実行の際は不要）
  ggplot2::theme_gray(base_family="HiraKakuPro-W3")

先ほどの６５歳以上人口総数のプロットと重ねてみる．

ggplot2::ggplot() +
  #最初に小地域境界を白塗りでプロット
  ggplot2::geom_sf(data=census_sa_miura,fill="white") +
  #６５歳以上人口総数が結合された500mメッシュをプロット
  ggplot2::geom_sf(data=mesh500m_miura,
                   #上で作成した分級
                   aes(fill=cut(６５歳以上人口総数,breaks=range)),
                   #不透明度80%で塗り分け
                   alpha=0.8) +
  #階級の色分けには「Reds」というパレットを使う
  ggplot2::scale_fill_brewer(palette="Reds") +
  #図の凡例のタイトルを手動で指定する
  ggplot2::guides(fill=guide_legend(title="６５歳以上人口総数")) +
  #色分けとポイントの大きさは下の行で設定したいので，aes内で色・サイズ分けの基準になる変数名（変数class）のみ指定
  ggplot2::geom_sf(data=med_miura,aes(color=class,size=class)) +
  #変数classに基づき，病院なら青，診療所なら緑，歯科診療所なら赤でプロット
  scale_color_manual(name="class",
                     breaks=c("病院","診療所","歯科診療所"),
                     values=c("blue","green","red")) +
  #変数classに基づき，病院ならサイズ3，診療所なら2，歯科診療所なら1でプロット
  scale_size_manual(name="class",
                    breaks=c("病院","診療所","歯科診療所"),
                    values=c(3,2,1)) +
  #日本語が文字化けしない為の設定（Windowsで実行の際は不要）
  ggplot2::theme_gray(base_family="HiraKakuPro-W3")

ステップ7：背景地図とバス停・ルートの追加

まずはOSMを背景地図として読み込んでみる．

ggplot2::ggplot() +
  #背景にOSMのタイルを読み込み
  ggspatial::annotation_map_tile(zoom=13) +
  #６５歳以上人口総数が結合された500mメッシュをプロット
  ggplot2::geom_sf(data=mesh500m_miura,
                   #上で作成した分級
                   aes(fill=cut(６５歳以上人口総数,breaks=range)),
                   #不透明度80%で塗り分け
                   alpha=0.8) +
  #階級の色分けには「Reds」というパレットを使う
  ggplot2::scale_fill_brewer(palette="Reds") +
  #図の凡例のタイトルを手動で指定する
  ggplot2::guides(fill=guide_legend(title="６５歳以上人口総数")) +
  #色分けとポイントの大きさは下の行で設定したいので，aes内で色・サイズ分けの基準になる変数名（変数class）のみ指定
  ggplot2::geom_sf(data=med_miura,aes(color=class,size=class)) +
  #変数classに基づき，病院なら青，診療所なら緑，歯科診療所なら赤でプロット
  scale_color_manual(name="class",
                     breaks=c("病院","診療所","歯科診療所"),
                     values=c("blue","green","red")) +
  #変数classに基づき，病院ならサイズ3，診療所なら2，歯科診療所なら1でプロット
  scale_size_manual(name="class",
                    breaks=c("病院","診療所","歯科診療所"),
                    values=c(3,2,1)) +
  #日本語が文字化けしない為の設定（Windowsで実行の際は不要）
  ggplot2::theme_gray(base_family="HiraKakuPro-W3")

続いて，バス路線と一緒にプロット．
※環境によってはエラーが出ることがあるが，もう一度実行すればエラーは表示されず正しくプロットされる．

#バスルートデータをメッシュが存在する領域（バウンディングボックス）のみに絞り込む
bus_routes_miura <- sf::st_crop(x=bus_routes_miura,
                                y=sf::st_bbox(obj=mesh500m_miura))
#バス停留所データをメッシュが存在する領域（バウンディングボックス）のみに絞り込む
bus_stops_miura <- sf::st_crop(x=bus_stops_miura,
                               y=sf::st_bbox(obj=mesh500m_miura))

ggplot2::ggplot() +
  #背景にOSMのタイルを読み込み
  ggspatial::annotation_map_tile(zoom=13) +
  #６５歳以上人口総数が結合された500mメッシュをプロット
  ggplot2::geom_sf(data=mesh500m_miura,
                   #上で作成した分級
                   aes(fill=cut(６５歳以上人口総数,breaks=range)),
                   #不透明度80%で塗り分け
                   alpha=0.8) +
  #階級の色分けには「Reds」というパレットを使う
  ggplot2::scale_fill_brewer(palette="Reds") +
  #図の凡例のタイトルを手動で指定する
  ggplot2::guides(fill=guide_legend(title="６５歳以上人口総数")) +
  #バスルートのプロット
  ggplot2::geom_sf(data=bus_routes_miura) +
  #バス停留所のプロット
  ggplot2::geom_sf(data=bus_stops_miura) +
  #色分けとポイントの大きさは下の行で設定したいので，aes内で色・サイズ分けの基準になる変数名（変数class）のみ指定
  ggplot2::geom_sf(data=med_miura,aes(color=class,size=class)) +
  #変数classに基づき，病院なら青，診療所なら緑，歯科診療所なら赤でプロット
  scale_color_manual(name="class",
                     breaks=c("病院","診療所","歯科診療所"),
                     values=c("blue","green","red")) +
  #変数classに基づき，病院ならサイズ3，診療所なら2，歯科診療所なら1でプロット
  scale_size_manual(name="class",
                    breaks=c("病院","診療所","歯科診療所"),
                    values=c(3,2,1)) +
  #日本語が文字化けしない為の設定（Windowsで実行の際は不要）
  ggplot2::theme_gray(base_family="HiraKakuPro-W3")

ステップ8：高齢者率の計算と地図の作成

高齢化率を計算し，変数として追加する．

mesh500m_miura <- mesh500m_miura %>%
  dplyr::mutate(高齢化率=６５歳以上人口総数/人口総数)

高齢化率をプロット．

ggplot2::ggplot() +
  #背景にOSMのタイルを読み込み
  ggspatial::annotation_map_tile(zoom=13) +
  #高齢化率が結合された500mメッシュをプロット
  ggplot2::geom_sf(data=mesh500m_miura,
                   #高齢化率で塗り分け
                   aes(fill=高齢化率),
                   #不透明度80%で塗り分け
                   alpha=0.8) +
  #最小値に近いほど白，近いほど濃い赤に近くなるよう塗り分け
  ggplot2::scale_fill_gradient(low="white",high="red4") +
  #バスルートのプロット
  ggplot2::geom_sf(data=bus_routes_miura) +
  #バス停留所のプロット
  ggplot2::geom_sf(data=bus_stops_miura) +
  #色分けとポイントの大きさは下の行で設定したいので，aes内で色・サイズ分けの基準になる変数名（変数class）のみ指定
  ggplot2::geom_sf(data=med_miura,aes(color=class,size=class)) +
  #変数classに基づき，病院なら青，診療所なら緑，歯科診療所なら赤でプロット
  scale_color_manual(name="class",
                     breaks=c("病院","診療所","歯科診療所"),
                     values=c("blue","green","red")) +
  #変数classに基づき，病院ならサイズ3，診療所なら2，歯科診療所なら1でプロット
  scale_size_manual(name="class",
                    breaks=c("病院","診療所","歯科診療所"),
                    values=c(3,2,1)) +
  #日本語が文字化けしない為の設定（Windowsで実行の際は不要）
  ggplot2::theme_gray(base_family="HiraKakuPro-W3")