1 目的

人口密度を加えたメッシュデータ（地図で見る統計）にコミュニティバス路線（国土数値情報）を重ねて，人口密度がバス事業成立の条件を満たしているかを可視化．

地方民間バス路線の廃止に伴い，コミュニティバスや乗合タクシーなどの形で自治体がバス事業を引き継ぐケースが見受けられる．
バス事業を成立（1便／1時間）するには，30人／ha程度の人口密度が求められる．
- 今回は500mメッシュ統計を利用．30人／haはメッシュ内750人程度に該当．
- 例として，富山県射水市のコミュニティバスを取り上げる．
  - 射水市コミュニティバス「きときとバス」（アクセス日：2022/05/06）
  - 参考文献：日立市交通交通協議（2015年）（アクセス日：2022/05/07）
  - （⇒日立市地域公共交通網形成計画（アクセス日：2024/03/25））

利用ライブラリ：sf，tidyverse, MetBrewer

library(sf)
library(tidyverse)

#Hokusai2カラーパレット利用のため（好み）
library(MetBrewer)

2 バス路線の可視化

富山県

国土数値情報ダウンロードサイト　>　行政区域（ポリゴン）　>　富山県　>　令和3年
- シェープファイルの読込（read_sf）．Toyama_mapと名付ける．

Toyama_map <-
  read_sf("N03-20210101_16_GML/N03-20210101_16_GML/N03-21_16_210101.shp")

射水市

パイプ演算子（%>%）を使いToyama_mapから「射水市」（Imizushi_shi_mapと名付ける）を選別（filter）．
- 射水市の場所を知らない人のために，その場所（ポリゴン）を富山県の地図上に示す．境界内（fill）は白．境界線は太く（linewidth）する．

#射水市の選別
Toyama_map %>%
  filter(N03_004=="射水市") -> 
  Imizu_shi_map 

#射水市の地図をImizuと命名
Imizu <-
  geom_sf(data=Imizu_shi_map, fill="white", 
          linewidth=0.8)
  
#射水市の可視化
ggplot()+
  geom_sf(data=Toyama_map)+
  Imizu

コミュニティバス路線

国土数値情報ダウンロードサイト　>　バスルート（ライン）　>　富山県　>　平成23年（のみ掲載，アクセス日：2022/05/06）¹
- このデータは日本語が文字化け．オプション（encoding="cp932"）で文字化けを防ぐ．

busline <-
  read_sf("N07-11_16_GML/N07-11_16.shp",
          options = c("encoding=CP932"),
          crs="EPSG:6668")

注意：ggplot()+geom_sf(data=busline)と指示してもエラーが出る．

バスルート（ライン）ファイルはGeometry type: MULTILINESTRINGでDimension:XYMのため地図に表せない模様（ここでは確認しないがhead(busline)でXYMであることが確認できる）．
- Mが3つ目の次元．2次元（Dimension:XY）にするためにst_xm()関数を使用して，Mを削除（head(busline)で再確認すると，Dimension:XYに変更される）．

#2次元に変更
busline <- st_zm(busline)

パイプ演算子（%>%）を使い，事業者（N07_002）の「射水市」を抜き出し，busline_Imizuと名付ける．
- バス路線（linetype）を破線（dashed）で表現．
- その他の形状：linetype参照．

#射水市コミュニティバスの選別
busline %>%
  filter(N07_002=="射水市") -> 
  busline_Imizu

#射水市コミュニティバス路線の地図をImizu_busと命名
Imizu_bus <-
  geom_sf(data=busline_Imizu, 
          linewidth=0.6, linetype="dashed")

#路線の可視化
ggplot()+
  Imizu+
  Imizu_bus+
  ggtitle("射水市コミュニティバス路線")

バス停留所

国土数値情報ダウンロードサイト　>　バス停留所（ポイント）　>　富山県　>　平成22年（のみ掲載，アクセス日：2022/05/06）²

busstop <-
  read_sf("P11-10_16_GML/P11-10_16-jgd-g_BusStop.shp", 
          crs="EPSG:6668")

射水市コミュニティバスのバス停のみ抜き出す．バス路線を運営する事業者名（コミュニティバスは自治体名）がP11_003_1に記されている．ただし，バス停が民間バスのバス停である場合もある．
- このため，P11_003_1==射水市と指定しても，不十分な数のバス停しか抜き出せない．
- ここでは，str_detect()を使い，P11_003_1列に「射水市」を含む文字列があればすべて抜き出すように指示．
  - 行の選択 - filter関数（アクセス日：2022/06/09）参考．

#射水市コミュニティバスのバス停の選別
busstop %>% 
  filter(str_detect(P11_003_1, "射水市")) ->
  busstop

ポイントデータはデフォルトが黒丸．枠と枠内の色を変えるため，shapeの21（円）を選び，枠を黒，枠内を白（fill）に指定．
- その他の形状：shape参照．

#射水市コミュニティバスのバス停の地図をImizu_busstopと命名
Imizu_busstop <-
  geom_sf(data=busstop, fill="white", 
          size=0.5, shape=21)

#バス停の可視化
ggplot()+
  Imizu+
  Imizu_bus+
  Imizu_busstop+
  ggtitle("射水市コミュニティバス路線")

鉄道路線

地図がわかりやすくなるように，鉄道路線を追加．

国土数値情報ダウンロードサイト　>　鉄道（ライン）　>　全国　>　令和元年（2019年）

railline <-
  read_sf("N02-19_GML/N02-19_RailroadSection.shp", 
                    options = c("encoding=CP932"))

あいの風とやま鉄道

あいの風とやま鉄道は富山県内をはしる第三セクター鉄道．北陸新幹線の延伸開業に伴い，JR西日本から経営分離された．

パイプ演算子（%>%）を使い，路線名（鉄道路線の名称，N02_003）の「あいの風とやま鉄道」を抜き出しAinokazeと名付ける．
- 線路（linetype）の形状をf8）にする．

#あいの風とやま鉄道の選別
railline %>%
  filter(N02_003=="あいの風とやま鉄道線") -> 
  Ainokaze

#あいの風とやま鉄道の地図をAinokaze_railwayと命名
Ainokaze_railway <-
  geom_sf(data=Ainokaze, linewidth=0.8,
          linetype="f8")

駅

あいの風とやま鉄道の名称と駅（呉羽，小杉，越中大門）の位置を表記．

annotate()：地図上の位置を指示．
- xは経度，yは緯度に対応．
- shapeの22（四角）を選び，枠を黒，枠内を白（fill）に指定．
- 緯度経度は地図検索サイトMapFanなどから特定．
geom_text()：地図上に文字の挿入を指示．

#あいの風とやま鉄道の可視化
ggplot()+
  Imizu+
  Ainokaze_railway+
  Imizu_bus+
  Imizu_busstop+
  annotate("point", x=137.0923306, y=36.7208275, 
           fill="white", size=2.5, shape=22)+
  annotate("point", x=137.0547234, y=36.734352, 
           fill="white", size=2.5, shape=22)+
  annotate("point", x=137.1648832, y=36.7188936, 
           fill="white", size=2.5, shape=22)+
  geom_text(aes(x=137.175, y=36.71), size=3, 
              label="あいの風とやま鉄道")+
  coord_sf(xlim=c(137.01, 137.24), 
           ylim=c(36.65, 36.79))+
  ggtitle("射水市コミュニティバス路線")

3 地域メッシュ統計による人口密度の可視化

境界データ：500メッシュ

統計地理情報システム　>　境界データダウンロード　>　4次メッシュ（500mメッシュ）　>　世界測地系緯度経度・Shapefile　>　都道府県で絞込みはコチラ　>　富山県
- 1県あたりシェープ・ファイルが4つほどあるが，今回は下記の2つで射水市をカバーできる．ここではmesh_2とmesh_4と名付ける．
- CRSがJGD2000のため，st_transform()で（JGD2011”EPSG:6668”）に揃える．
  - これをしないと，のちほどの作業できない．

mesh_2 <-
  read_sf("HDDSWH5437/MESH05437.shp") %>% 
  st_transform("EPSG:6668")

mesh_4 <-
  read_sf("HDDSWH5537/MESH05537.shp") %>% 
  st_transform("EPSG:6668")

メッシュの範囲が射水市をカバーしているか確認．
- 射水市（ポリゴン）の境界枠内を塗りつぶさず（fill，NAで透明に），境界を青（color，色は#1e4294），太線（linewidth）で示す．

ggplot()+
  geom_sf(data=mesh_2)+
  geom_sf(data=mesh_4)+
  geom_sf(data=Toyama_map, fill="NA")+
  geom_sf(data=Imizu_shi_map, fill="NA", 
          color="#1e4294", linewidth=0.8)

統計データ：500mメッシュ

e-stat　>　統計地理情報システム　>　統計データダウンロード　>　国勢調査　>　2015年　>　4次メッシュ（500mメッシュ）　>　その１　人口等基本集計に関する事項　>　都道府県で絞込みはコチラ　>　富山県　>　上記の番号と合致するファイル
人口データはテキスト形式．read.tableで読込．
- ここではPop_2とPop_4と名付ける．
- header=TRUEは列名が含まれていることを意味する．仮に列名がなければFALSE．
- sep=","はデータがカンマで区切られていることを意味する．skip=1は1番最初の行を削除することを意味する．

Pop_2 <- 
  read.table("tblT000847H5437/tblT000847H5437.txt",
                  header=TRUE, sep=",", skip=1,
                  fileEncoding="CP932")
Pop_4 <- 
  read.table("tblT000847H5537/tblT000847H5537.txt",
                  header=TRUE, sep=",", skip=1,
                   fileEncoding="CP932")

3.1 列名の変更

境界データと人口データに結合するためにrename()で1列目の列名をKEY_CODEに変更．
5列目の列名を「人口総数」として可視化する際に利用．

Pop_2 %>% 
  rename(KEY_CODE=1,
         人口総数=5) ->
  Pop_2

Pop_4 %>% 
  rename(KEY_CODE=1,
         人口総数=5) ->
  Pop_4

注意：境界データのKEY_CODEがcharacter（文字データ），人口データはinteger（数データ）のため結合できない．
- そこで，Pop_#のKEY_CODEをcharacterに変換（mutate）する．再びPop_#と名付ける．

Pop_2 %>%
  mutate(KEY_CODE=as.character(KEY_CODE)) ->
  Pop_2
Pop_4 %>%
  mutate(KEY_CODE=as.character(KEY_CODE)) ->
  Pop_4

3.2 データ結合1

境界データと人口データを結合する．

境界データ（mesh_#）と人口データ（Pop_#_1）をKEY_CODEで結合（left_join）し，改めてmeshPop_#と名付ける．

meshPop_2 <-
  left_join(mesh_2, Pop_2, by=c("KEY_CODE"))
meshPop_4 <-
  left_join(mesh_4, Pop_4, by=c("KEY_CODE"))

3.3 データ結合2

meshPop_#とImizu_shi_mapを結合．st_intersection()を用いて共通部分を抽出．mesh_intersection_#と名付ける．
- 属性がジオメトリ全体で空間的に一定と仮定しないと警告が出る？？？．一定（constant）と指定．

#Warning（警告）を出さないため
st_agr(meshPop_2)="constant"
st_agr(meshPop_4)="constant"
st_agr(Imizu_shi_map)="constant"

#射水市の地図に人口総数が加わったメッシュデータの完成
mesh_intersection_2 <-
  st_intersection(x=meshPop_2, y=Imizu_shi_map)
mesh_intersection_4 <-
  st_intersection(x=meshPop_4, y=Imizu_shi_map)

データ結合3

2つの結合データを行（縦）方向に結合．

人口総数を階級区分（グループ）などで表現したい場合，データが分かれている（現在，mesh_intersection_2とmesh_intersection_4）と不都合が生じる場合がある．そこで便宜上，一つにまとめる（rbind）．まとめたデータをImizu_mesh_popと名付ける．

Imizu_mesh_pop <-
  rbind(mesh_intersection_2, mesh_intersection_4)

人口密度の可視化

行方向に結合したデータ（Imizu_mesh_pop）を用いて，500Mメッシュ内の人口総数を地図に表現．

scale_fill_met_c()を利用して，連続変数である人口総数を色付け．
- 欠損値（NA）は色が付かない．

ggplot()+ 
  geom_sf(data=Imizu_mesh_pop, 
          aes(fill=人口総数))+
  scale_fill_met_c("Hokusai2")+
  ggtitle("射水市の人口分布")

人口総数NAの扱い

欠損値（NA）は人口が0と仮定．

Imizu_mesh_pop %>%
  mutate(人口総数=
           if_else(is.na(人口総数), 0, 人口総数)) ->  
  Imizu_mesh_pop

4 コミュニティバス路線と人口密度の可視化

これまで作成してきた地図を重ねて，目的の地図を完成する．

バスの採算性は30／ha．500メッシュのため25倍の大きさ．1マスのメッシュ内に750人必要．人口を250人で区切る（my_breaks）．
labsは判例の表記の変更のために利用．X=""，y=""を加えないと，横軸と縦軸にそれぞれx，yと表記される．

#人口の間隔
my_breaks<-c(0, 250, 500, 750, 1000, 1250)

#可視化 
ggplot()+
  geom_sf(data=Toyama_map, fill="white")+
  Imizu+
  geom_sf(data=Imizu_mesh_pop, 
          aes(fill=人口総数))+
  scale_fill_met_c("Hokusai2",
                   breaks=my_breaks)+
  Ainokaze_railway+
  Imizu_bus+
  Imizu_busstop+
  annotate("point", x=137.0923306, y=36.7208275, 
           fill="white", size=2.5, shape=22)+
  annotate("point", x=137.0547234, y=36.734352, 
           fill="white", size=2.5, shape=22)+
  annotate("point", x=137.1648832, y=36.7188936, 
           fill="white", size=2.5, shape=22)+
  geom_text(aes(x=137.175, y=36.71), size=3, 
              label="あいの風とやま鉄道")+
  coord_sf(xlim=c(137.01, 137.19), ylim=c(36.65, 36.79))+
  labs(fill="人", x="", y="")+
  ggtitle("射水市コミュニティバス路線と人口",
  subtitle="（2015年500mメッシュ人口）")

凡例

凡例がわかりづらいため，修正.

パイプ演算子（%>%）を使いImizu_mesh_popに新しい変数「人口総数2」を作成（mutate）し，改めてImizushi_mesh_pop2と名付ける．
- case_when()を利用して，条件に合致した値を言葉に置き換えられる（条件式~“変更名称”）．それ以外はTRUE．
- factor以降は順番を指定するための指令．

Imizu_mesh_pop %>%
  mutate(人口総数2=case_when(人口総数<250~"250人（10人／ha）未満",
                             人口総数<500~"500人（20人／ha）未満",
                             人口総数<750~"750人（30人／ha）未満",
                             人口総数<1000~"1000人（40人／ha）未満",
                             TRUE~"1000人（40人／ha）以上"),
         人口総数2=factor(人口総数2,
                          levels=c("250人（10人／ha）未満", 
                                   "500人（20人／ha）未満", 
                                   "750人（30人／ha）未満", 
                                   "1000人（40人／ha）未満", 
                                   "1000人（40人／ha）以上"))) -> 
  Imizu_mesh_pop2

完成図

先ほどまでは人口総数が連続変数であったが，人口総数2は階層化された離散変数のためscale_fill_met_d()を利用して，色付け．

ggplot()+
  geom_sf(data=Toyama_map, fill="white")+
  Imizu+
  geom_sf(data=Imizu_mesh_pop2, aes(fill=人口総数2))+
  scale_fill_met_d("Hokusai2")+
  Ainokaze_railway+
  Imizu_bus+
  Imizu_busstop+
  annotate("point", x=137.0923306, y=36.7208275, 
           fill="white", size=2.5, shape=22)+
  annotate("point", x=137.0547234, y=36.734352, 
           fill="white", size=2.5, shape=22)+
  annotate("point", x=137.1648832, y=36.7188936, 
           fill="white", size=2.5, shape=22)+
  geom_text(aes(x=137.175, y=36.71), size=3, 
              label="あいの風とやま鉄道")+
  coord_sf(xlim=c(137.01, 137.19), ylim=c(36.65, 36.79), 
           datum=NA)+
  labs(fill="人", x="", y="")+
  ggtitle("射水市コミュニティバス路線と人口",
  subtitle="（2015年500mメッシュ人口）")

5 応用例（250mメッシュ）

統計地理情報システムには5次メッシュ（250mメッシュ）も用意されている．

500mメッシュの代わりに250mメッシュによって表現．
- 色付けはHokusai3．
- 注意：PCによっては情報量が増えるためデータ処理に時間・負荷がかかる．
コードは繰り返しのため，非表示．

Rによる地理空間データの可視化

チュートリアルホーム

2024/03/25（改定日現在）2011年（平成23年）のほか2022年（令和4年）も用意されている．↩︎
2024/03/25（改定日現在）2010年度（平成22年度）のほか2022年度（令和4年度）も用意されている．↩︎

コミュニティバス路線上の人口密度の可視化

射水市の例

岩田真一郎（　神奈川大学経済学部）

2022-06-11作成，2024-03-25改訂

1 目的

2 バス路線の可視化

3 地域メッシュ統計による人口密度の可視化

3.1 列名の変更

3.2 データ結合1

3.3 データ結合2

4 コミュニティバス路線と人口密度の可視化

5 応用例（250mメッシュ）

コミュニティバス路線上の人口密度の可視化

射水市の例

岩田真一郎（ 神奈川大学経済学部）

2022-06-11作成，2024-03-25改訂

1 目的

2 バス路線の可視化

3 地域メッシュ統計による人口密度の可視化

3.1 列名の変更

3.2 データ結合1

3.3 データ結合2

4 コミュニティバス路線と人口密度の可視化

5 応用例 （250mメッシュ）

　岩田真一郎（　神奈川大学経済学部）

5 応用例（250mメッシュ）