環境ガバナンスのデータサイエンス/空間モデリング特論 R演習03
演習の概要
- 主題図の作成
注意事項
- 慶應義塾大学SFCで開講している「環境ガバナンスのデータサイエンス」「空間モデリング特論」の授業履修者向けの演習用ページです。
- 必ずしも全てのバージョンのRやOSで動作確認を行っていません。この演習用ページを作成している段階では、R3.5.3を使っています。
- Rの更新などにより、Rコードなどを予告無しに変更する場合があります。
データ分析の準備
パッケージのインストール
install.packages()でパッケージをインストールし、library()でパッケージを呼び出します。ここでは以下のパッケージを用います
install.packages(c("maptools","classInt","RColorBrewer", "readr",
"magrittr","dplyr","kokudosuuchi", "estatapi",
"sf","purrr", "needs", "DT"))
devtools::install_github("tidyverse/ggplot2")library(maptools)
library(classInt)
library(RColorBrewer)
library(readr)
library(magrittr)
library(dplyr)
library(kokudosuuchi)
library(estatapi)
library(sf)
library(purrr)
library(needs)
library(ggplot2)
library(DT)needsパッケージのmagrittr() はパイプ %>% の優先順位を高めるために用いられます
needs::prioritize(magrittr) GISメッシュデータのダウンロードと可視化
国土数値情報の活用
まず最初に、国土交通省の国土数値情報のウェブサイトから1kmメッシュ人口データをダウンロードします。さらに、神奈川県の2020年の将来人口推計データを可視化します。
国土数値情報ダウンロードサービスは以下のURLからアクセスできます。 http://nlftp.mlit.go.jp/ksj/index.html
kokudosuuchiパッケージのgetKSJSummary()を用いると、APIを使って国土数値情報の要約情報を得ることができます。
kokudosuuchi::getKSJSummary() %>%
filter(title=="1kmメッシュ別将来推計人口（H29国政局推計）（shape形式版）")## # A tibble: 1 x 5
## identifier title field1 field2 areaType
## <chr> <chr> <chr> <chr> <chr>
## 1 m1000 1kmメッシュ別将来推計人口（H29国政局推計）（shape形式版）… 各種統計… " " 3getKSJURL()関数で国土数値情報のURLを取得したあと、curlパッケージのcurl_download()関数で国土数値情報をダウンロードできます。
url <- kokudosuuchi::getKSJURL("m1000", prefCode = 14) %>%
dplyr::select(zipFileUrl) %>%
dplyr::pull(1)
tf <- tempfile(fileext = ".zip")
curl::curl_download(url, tf)
unzip(tf, exdir = getwd(), junkpaths = TRUE)
unlink(tf)
rm(tf)sfパッケージのst_read()関数でArcGISのシェープファイルを読み込みます。
d <- sf::st_read('Mesh3_POP_14.shp')ggplot2を使って、神奈川県の2020年将来人口推計データを以下のようにして可視化します。
ggplot2::ggplot(data = d, aes(fill=POP2020)) +
geom_sf() +
scale_fill_distiller(palette="RdYlGn") +
theme_bw() +
labs(title = "Population of Kanagawa in 2020 (est.)")
classIntパッケージのclassIntervals()関数のquantilestyleに従って定義されています。
他のclass styleを使いたい場合は、help(classIntervals)または以下のURLを参照してください。
http://web.sfc.keio.ac.jp/~maunz/wiki/index.php?asakura_sp_chap04
2020年推計人口の主題図を作成するためにの階級区分をn=5とした。またカラーパレットをfindClours()を用いて定義した。ここでは、グレースケールのカラーパレットを用いている。
pal0 <- c("white","grey","grey2")
ci <- classInt::classIntervals(d$POP2020, n=5, style = "quantile")
ci_Col <- classInt::findColours(ci,pal0)
境界線の色はborder=""により指定できる。
plot(sf::st_geometry(d), border="grey", col="white")
タイトルと凡例を、それぞれtitle()およびlegend()を使って指定し、主題図を以下のように描画することができる。
plot(sf::st_geometry(d), border="grey" ,col=ci_Col)
title("Population of Kanagawa in 2020 (Estimated)", cex=1.4)
legend("topleft",fill=attr(ci_Col,"palette"), cex=0.8,
legend=names(attr(ci_Col,"table")),bty="n")
以下の方法は、パイプ%>%を用いて主題図を作成する方法です。 この方法では、カラーパレットを作成するのにRColorBrewerを用いています。(以下のURLを参照 https://moderndata.plot.ly/create-colorful-graphs-in-r-with-rcolorbrewer-and-plotly/)
d %>%
sf::st_geometry(d) %>%
plot(col=classInt::findColours(ci, rev(RColorBrewer::brewer.pal(7, "RdYlGn"))))
グレースケールの主題図は、以下のように作成できる。
d %>%
sf::st_geometry(d) %>%
plot(col=classInt::findColours(ci, pal0), border="grey")
dplyrパッケージのfilter()関数を用いて、地域を抽出することができる。ここでは、市区町村コードを意味するCITY_CODEにより横浜市を抽出する。
「市区町村コード」は以下のURLを参照のこと。 http://www.soumu.go.jp/denshijiti/code.html
yoko <- dplyr::filter(d, CITY_CODE>=14101&CITY_CODE<=14143)
yoko %>%
sf::st_geometry(yoko) %>%
plot(col=classInt::findColours(ci, rev(RColorBrewer::brewer.pal(7, "RdYlGn"))))
500メートルメッシュデータを国土数値情報からダウンロードし、将来推計人口の主題図を作成する。
kokudosuuchi::getKSJSummary() %>%
filter(title=="500mメッシュ別将来推計人口（H29国政局推計）（shape形式版）")
url <- kokudosuuchi::getKSJURL("m500", prefCode = 14) %>%
select(zipFileUrl) %>%
pull(1)
tf <- tempfile(fileext = ".zip")
curl::curl_download(url, tf)
unzip(tf, exdir = getwd(), junkpaths = TRUE)
unlink(tf)
rm(tf)
d2 <- sf::st_read('Mesh4_POP_14.shp')d2 %>%
sf::st_geometry(d2) %>%
plot(col=classInt::findColours(
classInt::classIntervals(d2$POP2020, n=5, style = "sd"),
rev(RColorBrewer::brewer.pal(7, "RdYlGn"))),
border="grey")
視覚的に分かりやすい主題図を作成するために、横浜市の人口データに加え、鉄道や行政境界の地図をオーバーレイする。鉄道データや行政境界GISデータも国土数値情報から取得できる。
# kokudosuuchi::getKSJSummary() %>% filter(title=="鉄道")
url <- kokudosuuchi::getKSJURL("N02") %>%
select(zipFileUrl) %>%
pull(1)
tf <- tempfile(fileext = ".zip")
curl::curl_download(url, tf)
unzip(tf, exdir = getwd(), junkpaths = TRUE)
unlink(tf)
rm(tf)
rw_ln <- sf::st_read('N02-05-g_RailroadSection.shp')# kokudosuuchi::getKSJSummary() %>% filter(title=="行政区域")
kokudosuuchi::getKSJURL("N03", prefCode = 14) %>%
filter(year==2015) %>% select(zipFileUrl) %>%
pull(1) -> url
tf <- tempfile(fileext = ".zip")
curl::curl_download(url, tf)
unzip(tf, exdir = getwd(), junkpaths = TRUE)
unlink(tf)
rm(tf)
# Read a boundary of Kanagawa Pref.
ad <- sf::st_read('N03-15_14_150101.shp')
# Select Yokohama City boundary and plot it
ad %>%
filter(N03_007 %in% c(14101:14118)) -> y_ad# Select Yokohama city: City codes is from 14101 to 14118
y2 <- dplyr::filter(d2, CITY_CODE>=14101&CITY_CODE<=14118)
y2 %>%
sf::st_geometry(y2) %>%
plot(col=classInt::findColours(
classInt::classIntervals(y2$POP2020, n=5, style = "sd"),
rev(RColorBrewer::brewer.pal(7, "RdYlGn"))),
border="grey")
plot(rw_ln, add=TRUE)## Warning in plot.sf(rw_ln, add = TRUE): ignoring all but the first attributeplot(y_ad, col="transparent", border="red", add=TRUE)## Warning in plot.sf(y_ad, col = "transparent", border = "red", add = TRUE):
## ignoring all but the first attributetitle("Population of Yokohama CIty in 2020 (Estimated)")
legend("bottomright",
fill=attr(
classInt::findColours(
classInt::classIntervals(d$POP2020, n=5, style = "quantile"),
rev(RColorBrewer::brewer.pal(7, "RdYlGn"))),"palette"),
cex=1.0,
legend=names(attr(
classInt::findColours(
classInt::classIntervals(d$POP2020, n=5, style = "quantile"),
rev(RColorBrewer::brewer.pal(7, "RdYlGn"))),"table")),bty="n")
国土数値情報から複数のGISデータをダウンロードする
国土数値情報から複数の地域（ここでは、埼玉・千葉・東京・神奈川）の行政境界データをダウンロードし、関東地域の地図を作成してみよう。
preflist <- 11:14
urls <- kokudosuuchi::getKSJURL("N03", prefCode = preflist) %>%
filter(year==2015) %$% zipFileUrl
d <- purrr::map(urls, getKSJData,cache_dir = "cached_zip")
zip_exdir <- tempfile()
dir.create(zip_exdir)
setwd("cached_zip")
map(list.files(), unzip, exdir = zip_exdir)
mdata <- file.path(zip_exdir,
list.files(zip_exdir),
paste("N03-15_",11:14,"_150101.shp", sep="")) %>%
map(st_read) %>%
do.call(rbind, .)
setwd("../")
st_write(mdata, "N03-15_Kanto_150101.shp")シェープファイルと属性テーブルの結合
次に、東京２３区の行政境界GISデータと、住宅地地価（公示地価）データを結合して主題図を作成する方法を示します。
2018年の東京都の公示地価データは、東京都のウェブサイトからダウンロードできる。ここで、地価データは地点データの属性情報として提供されている。
CSVファイル（カンマ区切りファイル）はreadrのread_csv()関数を用いて読み込むことがでる。
select()関数とfilter()関数を用いて、住宅地の地価に関する地点データを選択することができる。
lp_url <- sprintf("http://www.zaimu.metro.tokyo.jp/kijunchi/30kouji/30kouji_chiten_data.csv")
curl::curl_download(lp_url, "lp_tky.csv")
lph_tky <- readr::read_csv("lp_tky.csv", skip=1,
col_types = cols(都道府県市区町村コード = "c"),
locale = locale(encoding = "CP932")) %>%
select(CODE = 2,
LP2018 = 10,
Chiseki = 13,
Yoto = 34) %>%
filter(Yoto %in% c("１住居", "１中専", "１低専", "２住居", "２中専", "２低専", "準住居"))さらに、区単位で平均地価を計算するために、dplyパッケージのsummarize()関数を用いて、CODE毎の統計量を集約する。エクセルのピボットテーブルのようなものです。
lph_tky %>%
dplyr::group_by(CODE) %>%
dplyr::summarise(LP2018.Mean = mean(LP2018))東京都の行政境界も国土数値情報から入手可能です。
kokudosuuchi::getKSJURL("N03", prefCode = 13) %>%
filter(year==2015) %>% select(zipFileUrl) %>%
pull(1) -> url
tf <- tempfile(fileext = ".zip")
curl::curl_download(url, tf)
unzip(tf, exdir = getwd(), junkpaths = TRUE)
unlink(tf)
rm(tf)東京都の行政境界を読み込み、２３区の行政境界を抽出する。
tky_bd <- sf::st_read('N03-15_13_150101.shp') %>%
dplyr::group_by(N03_007) %>%
dplyr::summarise(geometry = st_union(geometry)) %>%
dplyr::ungroup()
tky_bd %>%
filter(N03_007 %in% c(13101:13123)) -> tky23ku_bddplyrパッケージのinner_join()関数を用いてシェープファイルとテーブルファイルを結合します。
# Join a shape file and a table
tky23ku_lph <- dplyr::inner_join(tky23ku_bd, lph_tky, by = c("N03_007" = "CODE"))## Warning: Column `N03_007`/`CODE` joining factor and character vector,
## coercing into character vector最後に、2018年の東京23区の住宅地地価データ（平均値）の主題図を以下のように作成します。
plot(tky23ku_lph["LP2018"])
総務省統計局のe-statから統計データをダウンロードする
e-Stat URL: https://www.e-stat.go.jp/
データをダウンロードする前に、e-StatのAPIを以下のWebサイト上で登録する。
https://www.e-stat.go.jp/mypage/user/preregister
ここで登録した、あなたのe-Stat AIP IDをmyappIDに入力する必要があります。
e-Statの統計表IDから統計表を検索することができます。 estatapiパッケージのestat_getStatsList()関数を用いて、統計表IDのリストを返すことができます。
# Search stats list by stats data ID of e-stat if known
StatsList <- estatapi::estat_getStatsList(
appId = myappID,
lang = c("J"),
searchWord = "",
statsDataID = "0000010109"
)
str(StatsList)## Classes 'tbl_df', 'tbl' and 'data.frame': 10000 obs. of 13 variables:
## $ @id : chr "0003110904" "0003110914" "0003110917" "0003110918" ...
## $ STAT_NAME : chr "国民経済計算" "国民経済計算" "国民経済計算" "国民経済計算" ...
## $ GOV_ORG : chr "内閣府" "内閣府" "内閣府" "内閣府" ...
## $ STATISTICS_NAME : chr "国民経済計算年次推計 平成2年基準 平成10年度確報 フロー編" "国民経済計算年次推計 平成2年基準 平成10年度確報 フロー編" "国民経済計算年次推計 平成2年基準 平成10年度確報 フロー編" "国民経済計算年次推計 平成2年基準 平成10年度確報 フロー編" ...
## $ TITLE : chr "1.統合勘定 2. 国民可処分所得と処分勘定 年度" "1.統合勘定 1. 国内総生産と総支出勘定 年度" "1.統合勘定 2. 国民可処分所得と処分勘定 暦年" "1.統合勘定 3. 資本調達勘定 (1)実物取引 暦年" ...
## $ CYCLE : chr "年次" "年次" "年次" "年次" ...
## $ SURVEY_DATE : chr "199801-199812" "199801-199812" "199801-199812" "199801-199812" ...
## $ OPEN_DATE : chr "2015-09-30" "2015-09-30" "2015-09-30" "2015-09-30" ...
## $ SMALL_AREA : chr "0" "0" "0" "0" ...
## $ MAIN_CATEGORY : chr "企業・家計・経済" "企業・家計・経済" "企業・家計・経済" "企業・家計・経済" ...
## $ SUB_CATEGORY : chr "国民経済計算" "国民経済計算" "国民経済計算" "国民経済計算" ...
## $ OVERALL_TOTAL_NUMBER: chr "616" "748" "616" "396" ...
## $ UPDATED_DATE : chr "2016-08-15" "2016-08-15" "2016-08-15" "2016-08-15" ...estat_getMetaInfo()関数は、統計情報のメタデータを返します。
# Get metadata information
meta_info <- estatapi::estat_getMetaInfo(
appId = myappID,
statsDataId = "0000010109"
)
meta_info## $tab
## # A tibble: 1 x 3
## `@code` `@name` `@level`
## <chr> <chr> <chr>
## 1 00001 観測値 1
##
## $cat01
## # A tibble: 384 x 4
## `@code` `@name` `@level` `@unit`
## <chr> <chr> <chr> <chr>
## 1 I1101 I1101_平均余命（０歳）（男） 1 年
## 2 I1102 I1102_平均余命（０歳）（女） 1 年
## 3 I1201 I1201_平均余命（２０歳）（男） 1 年
## 4 I1202 I1202_平均余命（２０歳）（女） 1 年
## 5 I1301 I1301_平均余命（４０歳）（男） 1 年
## 6 I1302 I1302_平均余命（４０歳）（女） 1 年
## 7 I1401 I1401_平均余命（６０歳）（男） 1 年
## 8 I1402 I1402_平均余命（６０歳）（女） 1 年
## 9 I1501 I1501_平均余命（６５歳）（男） 1 年
## 10 I1502 I1502_平均余命（６５歳）（女） 1 年
## # ... with 374 more rows
##
## $area
## # A tibble: 48 x 3
## `@code` `@name` `@level`
## <chr> <chr> <chr>
## 1 00000 全国 1
## 2 01000 北海道 2
## 3 02000 青森県 2
## 4 03000 岩手県 2
## 5 04000 宮城県 2
## 6 05000 秋田県 2
## 7 06000 山形県 2
## 8 07000 福島県 2
## 9 08000 茨城県 2
## 10 09000 栃木県 2
## # ... with 38 more rows
##
## $time
## # A tibble: 42 x 3
## `@code` `@name` `@level`
## <chr> <chr> <chr>
## 1 1975100000 1975年度 1
## 2 1976100000 1976年度 1
## 3 1977100000 1977年度 1
## 4 1978100000 1978年度 1
## 5 1979100000 1979年度 1
## 6 1980100000 1980年度 1
## 7 1981100000 1981年度 1
## 8 1982100000 1982年度 1
## 9 1983100000 1983年度 1
## 10 1984100000 1984年度 1
## # ... with 32 more rows
##
## $.names
## # A tibble: 4 x 2
## id name
## <chr> <chr>
## 1 tab 観測値
## 2 cat01 Ｉ 健康・医療
## 3 area 地域
## 4 time 調査年estat_getStatsData()関数を用いて統計データを取得できる。 statsDataId、cdCat01あるいはlvAreaで条件を指定すれば、検索条件を狭めることができます。
# Get stats from e-Stat
SDS_df <- estatapi::estat_getStatsData(
appId = myappID,
statsDataId = "0000010109", # 社会・人口統計体系 都道府県データ 健康・医療
cdCat01 = "I1601", # 健康寿命
# cdCat01 = c("I1601", "I1602"), # 健康寿命
lvArea = 2 # Level of Area
# lvTime = 1 # Level of Time
# cdtime = "2010100000" # 2010年度
# surveyYears = 2010
) ## Fetching record 1-94... (total: 94 records)h2013_df <- subset(SDS_df, SDS_df$調査年 %in% "2013年度") %>%
as.data.frame()
h2013_df$KENCODE <- as.numeric(h2013_df$area_code)/1000
head(h2013_df)## tab_code 観測値 cat01_code Ｉ 健康・医療 area_code 地域
## 1 00001 観測値 I1601 I1601_健康寿命（男） 01000 北海道
## 2 00001 観測値 I1601 I1601_健康寿命（男） 02000 青森県
## 3 00001 観測値 I1601 I1601_健康寿命（男） 03000 岩手県
## 4 00001 観測値 I1601 I1601_健康寿命（男） 04000 宮城県
## 5 00001 観測値 I1601 I1601_健康寿命（男） 05000 秋田県
## 6 00001 観測値 I1601 I1601_健康寿命（男） 06000 山形県
## time_code 調査年 unit value KENCODE
## 1 2013100000 2013年度 年 71.11 1
## 2 2013100000 2013年度 年 70.29 2
## 3 2013100000 2013年度 年 70.68 3
## 4 2013100000 2013年度 年 71.99 4
## 5 2013100000 2013年度 年 70.71 5
## 6 2013100000 2013年度 年 71.34 6str(h2013_df)## 'data.frame': 47 obs. of 11 variables:
## $ tab_code : chr "00001" "00001" "00001" "00001" ...
## $ 観測値 : chr "観測値" "観測値" "観測値" "観測値" ...
## $ cat01_code : chr "I1601" "I1601" "I1601" "I1601" ...
## $ Ｉ 健康・医療: chr "I1601_健康寿命（男）" "I1601_健康寿命（男）" "I1601_健康寿命（男）" "I1601_健康寿命（男）" ...
## $ area_code : chr "01000" "02000" "03000" "04000" ...
## $ 地域 : chr "北海道" "青森県" "岩手県" "宮城県" ...
## $ time_code : chr "2013100000" "2013100000" "2013100000" "2013100000" ...
## $ 調査年 : chr "2013年度" "2013年度" "2013年度" "2013年度" ...
## $ unit : chr "年" "年" "年" "年" ...
## $ value : num 71.1 70.3 70.7 72 70.7 ...
## $ KENCODE : num 1 2 3 4 5 6 7 8 9 10 ...全国の都道府県境界GISデータは、ここからダウンロードできます。
# Download GIS data set from my web page
url.gis <- sprintf("http://web.sfc.keio.ac.jp/~maunz/asakura_sp/asakura_sp_data_1101.zip")
tf <- tempfile(fileext = ".zip")
curl::curl_download(url.gis, tf)
unzip(tf, exdir = getwd(), junkpaths = TRUE)
# Read jpn_pref.shp = prefecture boundary data
jpn_pref <- sf::st_read('jpn_pref.shp')最後に、都道府県境界GISデータと健康寿命データを結合させることにより、主題図を作成することができる。ここでは、2013年の都道府県レベルでの健康寿命に関する主題図を作成する。
# Join SDS health care data and jpn_pref.shp
jpn_pref_h <- dplyr::inner_join(jpn_pref, h2013_df, by = c("PREF_CODE" = "KENCODE"))
# Plot a map
ggplot(data = jpn_pref_h, aes(fill= value)) +
geom_sf() +
scale_fill_distiller(palette="RdYlGn") +
theme_bw() +
labs(title = "Health life expectancy of male in Japan (2013)")