1 目的

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

  • 地方民間バス路線の廃止に伴い,コミュニティバスや乗合タクシーなどの形で自治体がバス事業を引き継ぐケースが見受けられる.

  • バス事業を成立(1便/1時間)するには,30人/ha程度の人口密度が求められる.

利用ライブラリ:sftidyverse, MetBrewer

library(sf)
library(tidyverse)

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

2 バス路線の可視化

富山県

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)1

    • このデータは日本語が文字化け.オプション(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: MULTILINESTRINGDimension: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("射水市コミュニティバス路線")

バス停留所

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列に「射水市」を含む文字列があればすべて抜き出すように指示.
#射水市コミュニティバスのバス停の選別
busstop %>% 
  filter(str_detect(P11_003_1, "射水市")) ->
  busstop
  • ポイントデータはデフォルトが黒丸.枠と枠内の色を変えるため,shape21(円)を選び,枠を黒,枠内を白(fill)に指定.
    • その他の形状:shape参照.
#射水市コミュニティバスのバス停の地図をImizu_busstopと命名
Imizu_busstop <-
  geom_sf(data=busstop, fill="white", 
          size=0.5, shape=21)

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

鉄道路線

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

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は緯度に対応.
    • shape22(四角)を選び,枠を黒,枠内を白(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_2mesh_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")
  • メッシュの範囲が射水市をカバーしているか確認.
    • 射水市(ポリゴン)の境界枠内を塗りつぶさず(fillNAで透明に),境界を青(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メッシュ) > その1 人口等基本集計に関する事項 > 都道府県で絞込みはコチラ > 富山県 > 上記の番号と合致するファイル

  • 人口データはテキスト形式.read.tableで読込.

    • ここではPop_2Pop_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_CODEcharacter(文字データ),人口データはinteger(数データ)のため結合できない.
    • そこで,Pop_#KEY_CODEcharacterに変換(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_2mesh_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=""を加えないと,横軸と縦軸にそれぞれxyと表記される.

#人口の間隔
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による地理空間データの可視化


  1. 2024/03/25(改定日現在)2011年(平成23年)のほか2022年(令和4年)も用意されている.↩︎

  2. 2024/03/25(改定日現在)2010年度(平成22年度)のほか2022年度(令和4年度)も用意されている.↩︎