library(sf)
library(tidyverse)

#インタラクティブな地図作成
library(mapview)

1 目的

ある地物（地上にある物）や地点周辺の統計地図の作成，および統計量（合計値，平均値など）の計算．ここでは，多摩ニュータウン（ポイントデータ）周辺の人口（メッシュデータ）を可視化．また多摩ニュータウンから500mの円を描き，その円と重なるメッシュデータの人口数を計算．

知りたい統計量の例

駅から1km（半径500m）内のメッシュデータの住宅数，ポイントデータの平均地価．
危険箇所と重なるメッシュデータの人口数．

2 多摩ニュータウンと多摩市周辺の人口分布

2.1 ダウンロード

行政区域

国土数値情報ダウンロードサービス　>　行政区域（ポリゴン）　>　東京都　>　世界測地系　令和4年
- ダウンロードしたZIPファイルは同じフォルダ（ディレクトリ）に納め，ZIPファイルを選択 >　右クリック　>　すべて展開　>　展開先の選択とファイルの選択　>　展開を実行する．展開後は最初にダウンロードしたZIPファイルを削除．
- シェープファイルの読込（read_sf()）：行政区域をTokyo_mapと名付ける．
  - 多摩市周辺の行政区域の抜き出し（filter()）：Tama_shi_mapと名付ける．

#東京の行政区域（令和4年）
Tokyo_map <-
  read_sf("N03-20220101_13_GML/N03-22_13_220101.shp")

#多摩市周辺（多摩，町田，日野，稲城，八王子）
Tokyo_map %>%
  filter(N03_007=="13224" | N03_007=="13209" |
           N03_007=="13212" | N03_007=="13225" |
           N03_007=="13201") ->
  Tama_shi_map

#多摩市の行政区域の可視化
ggplot()+
  geom_sf(data=Tama_shi_map)

ニュータウン

国土数値情報ダウンロードサービス　>　ニュータウン（ポイント）　>　東京都　>　世界測地系　平成25年
- シェープファイルの読込（read_sf()）：行政区域をTokyo_newtownと名付ける．
- missing crsのエラーが出るため，crs（座標参照系）を設定．
  - WGS84：1984年に定められた世界測地系．
- 多摩ニュータウンを抜き出し，Tama_newtownと名付ける．

#東京ニュータウン
Tokyo_newtown <-
  read_sf("P26-13_13/P26-13_13/P26-13_13.shp", 
          crs="WGS84")

#多摩ニュータウンの選抜
Tokyo_newtown %>%
  filter(P26_005=="多摩ニュータウン") ->
  Tama_newtown

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

市区町村別メッシュ・コード一覧（残念ながら3次メッシュの情報だけの模様）から多摩市のコードが5339からはじまることがわかる．

人口

統計地理情報システム　>　統計データダウンロード　>　国勢調査　>　2020年　>　4次メッシュ（500mメッシュ）　>　人口及び世帯　>　都道府県で絞込みはコチラ　>　13 東京都　>　M5339 ¹

read.table()を用いてデータを読込．Tokyo_popと名付ける．

#人口
Tokyo_pop<-read.table("tblT001101H5339/tblT001101H5339.txt",
                      header=TRUE, sep=",",
                      fileEncoding="CP932" )

#便宜上1行目削除（データ名が2行にわたるため）
Tokyo_pop %>% 
  slice(-1) ->
  Tokyo_pop

#メッシュデータと結合するため，KEY_CODEのデータ型変換
Tokyo_pop %>%
  mutate(KEY_CODE=as.character(KEY_CODE)) ->
  Tokyo_pop

境界

メッシュデータの人口総数を地図上に可視するには，シェープファイルと結合する必要がある．そこで，メッシュデータの境界データをダウンロードし，結合する．

統計地理情報システム　>　境界データダウンロード　>　4次メッシュ（500mメッシュ）　>　世界測地系緯度経度・Shapefile　>　都道府県で絞込みはコチラ　>　13 東京都　>　M5339
- read.sf()を用いてデータを読込．meshと名付ける
- meshとTokyo_popをKEY_CODEで結合．

mesh <-
  read_sf("HDDSWH5339/MESH05339.shp")

#データ結合
Tokyo_mesh <-
  left_join(mesh, Tokyo_pop, by=c("KEY_CODE"))

#人口総数（T001101001）を数値（integer：整数）に変更
Tokyo_mesh %>%
  mutate(T001101001=as.integer(T001101001)) ->
  Tokyo_mesh

2.2 多摩NTと周辺人口分布の可視化

Tokyo_meshは東京都全体をカバー．ここでは，多摩市周辺だけを示したい．そこで，多摩市周辺の行政区域（Tama_shi_map）とTokyo_meshの重複する部分を抽出する．

st_intersection()を用いて，2つのシェープファイル（Tokyo_meshとTama_shi_map）の重複部分を抽出．
- そのためには以下の準備で示されている指示が必要．²

#準備（投影変換→WGS84）
Tokyo_mesh_84<-
  Tokyo_mesh %>%
  st_transform("WGS84") 

Tama_shi_84<-
  Tama_shi_map %>% 
  st_transform("WGS84") 

st_agr(Tokyo_mesh_84)="constant"
st_agr(Tama_shi_84)="constant"

#多摩周辺の地図に人口が加わったメッシュデータの完成
Tama_mesh <-
  st_intersection(x=Tokyo_mesh_84, y=Tama_shi_84)

ggplot()

#完成図
ggplot()+
  geom_sf(data=Tama_mesh, aes(fill=T001101001))+
  scale_fill_viridis_c(option="G", direction=-1)+
  geom_sf(data=Tama_newtown, fill="yellow", color="yellow")+
  geom_sf(data=Tama_shi_map, fill="NA",
          linewidth=0.8)+
  labs(fill="人口総数",
       caption="出典：国土数値情報（国土交通省）,
       国勢調査4次メッシュ（総務省統計局）")+
  ggtitle("多摩NT位置と多摩市周辺人口（2020年）")+
  theme_void()

mapview()

mapview()を用いる利点：マウスカーソルを自ら動かすことで知りたいメッシュ内の人口総数を調べられるようになる．³

layer.nameを用いて凡例を変更．
- 凡例にあるNA（欠損値）は情報（ここでは人口総数）が含まれないことを意味する．

#viridisカラーの使用（好み）
library(scales)

#mapview表記
#col.regionsは好み（削除しても問題ない）
mapview(Tama_mesh, zcol="T001101001", 
        alpha.regions=0.6,
        col.region=viridis_pal(option="G",
                               direction=-1),
        layer.name="人口総数")+
  mapview(Tama_newtown, cex=2.5,
          color="yellow", col.regions="yellow")

3 周辺500mの可視化

ニュータウン（ポイント）の周辺500mを可視化．
- st_buffer()：円（ポリゴン）を作成．ここでは，Tama_newtownのポイント（geometry列に納められている）を中心として500mの円を作成．bufferと名付ける．

#単位を付与（ここでは周辺500m）
library(units)

Tama_newtown %>%
  st_buffer(geometry, dist=set_units(500, m)) ->
  buffer

#可視化1（500mの円）
ggplot()+
  geom_sf(data=Tama_mesh, aes(fill=T001101001))+
  scale_fill_viridis_c(option="G", direction=-1)+
  geom_sf(data=buffer, color="deeppink", 
          fill="NA", linewidth=0.7)+
  geom_sf(data=Tama_newtown, fill="yellow", color="yellow")+
  labs(fill="人口総数",
       caption="出典：国土数値情報（国土交通省）,
       国勢調査4次メッシュ（総務省統計局）")+
  ggtitle("多摩ニュータウン周辺500mの範囲")+
  theme_void()

4 人口総数（の合計値）の計算

st_intersects()：Tokyo_mesh_84（メッシュデータ）から円（buffer）ポリゴンと重なるメッシュを検索．
lengths()>0: 重なるメッシュを抽出しTRUEに．
- NT_Popと名付ける．

#エラーが出ないためのおまじない
sf_use_s2(FALSE)

NT_Pop <-
  lengths(st_intersects(Tokyo_mesh_84, buffer))>0

Tokyo_mesh_84からNT_Pop==TRUEを満たすメッシュを抽出（filter()）．
- 第1引数.ですべての変数を残す．
- NT_meshと名付ける．

Tokyo_mesh_84 %>%
  filter(., NT_Pop==TRUE) ->
  NT_mesh

#可視化2（円と重複するメッシュの可視化）
ggplot()+
  geom_sf(data=Tama_mesh, aes(fill=T001101001))+
  scale_fill_viridis_c(option="G", direction=-1)+
  geom_sf(data=buffer, color="deeppink", 
          fill="NA")+
  geom_sf(data=NT_mesh, color="pink", 
          fill="NA", linewidth=0.7)+
  geom_sf(data=Tama_newtown, fill="yellow", color="yellow")+
  labs(fill="人口総数",
       caption="出典：国土数値情報（国土交通省）,
       国勢調査4次メッシュ（総務省統計局）")+
  ggtitle("多摩ニュータウン周辺メッシュ")+
  theme_void()

周辺500mのメッシュに含まれる人口総数の計算

上の図の通り，重複するデータ（ここではメッシュ）が存在すると，統計量を計算できない．重複するデータはdistinct()で削除．
- 欠損値（NA）が含まれる場合も統計値を計算できない．欠損値を無視して計算する（na.rm=TRUE）．

#重複するメッシュ（KEY_CODE）を取り除く
#（.keep_all=TRUE）その他の変数残す
#欠損値無視オプション（na.rm）

NT_mesh %>% 
  distinct(KEY_CODE,.keep_all=TRUE) %>%
  summarize(sum_pop=sum(T001101001, na.rm=TRUE))

sum_pop	geometry
147078	MULTIPOLYGON (((139.3562 35…

mapview()で確認

mapview()を用いて，周辺500mと重なるメッシュのみを可視化．

lwd：枠の太さを指示．
凡例（legend）とホームボタン（homebutton）を削除（FALSE）．

#mapview
mapview(NT_mesh, zcol="T001101001", 
          alpha.regions=0.6, 
          col.region=viridis_pal(option="G",
                                 direction=-1),
          layer.name="人口総数")+
mapview(buffer, color="deeppink", 
        alpha.regions=0, lwd=1.2,
        legend=FALSE, homebutton=FALSE)+
  mapview(Tama_newtown, cex=2.5,
          color="yellow", col.regions="yellow",
          legend=FALSE, homebutton=FALSE)

謝辞

ポリゴンと重なるメッシュの検索，抽出方法について山本雅資氏にRコードを教示いただいた．

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

チュートリアルホーム

統計地理情報システムの国勢調査からは1995年から2020年（5年ごと）までの人口総数（男女），世帯総数などがわかる（アクセス日：2022-10-29）．2015年からは年齢別人口，外国人人口など詳しく記録されている．↩︎
準備に示されているコードはコミュニティバス路線上の人口密度の可視化に説明あり．↩︎
mapview()を用いた地理空間データの可視化についてはインタラクティブな地図によるデータの可視化参照．↩︎

地物周辺の地理空間データの可視化

多摩ニュータウン周辺500mの総人口

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

2022-10-29作成，2024-04-13改訂

1 目的

2 多摩ニュータウンと多摩市周辺の人口分布

2.1 ダウンロード

2.2 多摩NTと周辺人口分布の可視化

3 周辺500mの可視化

4 人口総数（の合計値）の計算

地物周辺の地理空間データの可視化

多摩ニュータウン周辺500mの総人口

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

2022-10-29作成，2024-04-13改訂

1 目的

2 多摩ニュータウンと多摩市周辺の人口分布

2.1 ダウンロード

2.2 多摩NTと周辺人口分布の可視化

3 周辺500mの可視化

4 人口総数（の合計値）の計算

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