このRPubsでは，Rによる市場アクセス分析で用いられるデータセットの整備を行います．Rを用いて実践できる以下のテクニックを網羅しています．
RPubsを作成する際には，手動でファイルをダウンロード・解凍したり，それらをローカル環境から読み込んだりする作業を極力削減するという点を特に重視しました．
このWEBサイトに設置したリンクは，新しいタブ・ウィンドウで開くことを推奨します．バグの影響で，新しいタブ・ウィンドウでないと開けないリンクが多いためです．予めご容赦ください．
Rによる空間分析の基礎事項を知りたい場合は，以下の文献が参考になります．このページで用いているテクニックは，「『事例で学ぶ経済・政策分析のためのGIS入門』のためのR入門」でおおむねカバーされています．
このRPubsではAPIやWebスクレイピングの詳細な解説は行いませんが，基礎知識を身に付けたい場合は以下の文献が参考になります．このページで用いているテクニックは，「Rでスクレイピングするならrvest」でおおむねカバーされています．
以下で用いられるデータの原典は以下の通りです．
今回の演習では，以下のパッケージを用います．未インストールのものについては，適宜インストールを行ってください．インストールの過程で「Do you want to install from sources the packages which need compilation? (Yes/no/cancel)」というメッセージが表示された場合には，とりあえず「no」を選択することをお勧めします．
library(tidyverse)
library(sf)
library(openxlsx)
library(estatapi)
library(magrittr)
library(xml2)
library(rvest)
library(geojsonio)
library(nngeo)動作環境によっては，sfパッケージの読み込み時にエラーが出る場合があります．その原因は，sfパッケージが依存しているRcppのバージョンが古いことが多いようです．もしエラーが出る場合には，以下のコードを実行して見てください（筆者のPCではエラーが出なかったので，コメント化してあります）．
# install.packages("Rcpp")
# library(Rcpp)estatapiパッケージを用いて，e-Stat（政府統計の総合窓口）から政府統計データをダウンロードし，それらを分析に用いることができる形に加工する方法を示します．e-Stat上のデータは，e-Stat APIというサービスを通じ，RやPython等のソフトウェア上でのコマンド操作を用いて，直接ダウンロードすることができます．一方，APIの扱いには若干の専門知識が必要となるので，初心者が用いるにはハードルが高いです．estatapiは，そうした専門知識がなくても，APIを用いたデータのダウンロードを可能にするパッケージです．
事前にこちらのpdfに示した手順に従い，APIキーの発行を済ませておいてください． APIキーは以後何度も繰り返して使うので，いちいち入力しないで済むよう，オブジェクトappIdに格納します．格納するキーは，ご自身で取得したものに必ず置き換えて使うようにしてください．
appId <- "XXXXXXXXX"estatapiのestat_getStatsList関数を用いて，全ての社会・人口統計体系の統計表を検索します．統計表の検索方法はいくつかありますが，今回はキーワード検索を用います．具体的には，引数searchWordに「社会・人口統計体系」を指定します．引数appIdには先程入力したappIdを指定します．
ちなみに，別の検索方法として，政府統計コードを用いることもできます．その際には，こちらのページにある政府統計コードの一覧から，自分が取得したい統計のコードを調べ，引数statsCodeでそちらを指定してください．
tbl_info <- estatapi::estat_getStatsList(appId=appId,
searchWord="社会・人口統計体系",
searchKind=1)
head(tbl_info)## # A tibble: 6 × 19
## `@id` STAT_…¹ GOV_ORG STATI…² TITLE CYCLE SURVE…³ OPEN_…⁴ SMALL…⁵ COLLE…⁶
## <chr> <chr> <chr> <chr> <chr> <chr> <chr> <chr> <chr> <chr>
## 1 0000010101 社会・… 総務省 都道府… Ａ … 年度… 0 2023-0… 0 全国
## 2 0000010101 社会・… 総務省 都道府… Ａ … 年度… 0 2023-0… 0 都道府…
## 3 0000010102 社会・… 総務省 都道府… Ｂ … 年度… 0 2023-0… 0 全国
## 4 0000010102 社会・… 総務省 都道府… Ｂ … 年度… 0 2023-0… 0 都道府…
## 5 0000010103 社会・… 総務省 都道府… Ｃ … 年度… 0 2023-0… 0 全国
## 6 0000010103 社会・… 総務省 都道府… Ｃ … 年度… 0 2023-0… 0 都道府…
## # … with 9 more variables: MAIN_CATEGORY <chr>, SUB_CATEGORY <chr>,
## # OVERALL_TOTAL_NUMBER <chr>, UPDATED_DATE <chr>, TABULATION_CATEGORY <chr>,
## # TABULATION_SUB_CATEGORY1 <chr>, TABULATION_CATEGORY_EXPLANATION <chr>,
## # TABLE_NAME <chr>, TABULATION_SUB_CATEGORY_EXPLANATION1 <chr>, and
## # abbreviated variable names ¹STAT_NAME, ²STATISTICS_NAME, ³SURVEY_DATE,
## # ⁴OPEN_DATE, ⁵SMALL_AREA, ⁶COLLECT_AREA市区町村集計の統計表のリストを作成します．その際，合併処理後の統計表のみに候補を絞ります．
tbl_info_muni <- tbl_info %>%
#集計単位が市区町村のものに絞る
dplyr::filter(COLLECT_AREA=="市区町村") %>%
#合併処理後の統計表のみ残す
dplyr::filter(grepl("廃置分合処理済",STATISTICS_NAME))社会・人口統計体系は，統計表の生の値を収録した「基礎データ」と，基礎データ内の複数の変数を用いて計算される各種指標（例：可住地面積と総人口を用いて計算される「単位可住地面積あたり人口（人口密度）」）を収録した「社会生活統計指標」に分かれます．今回用いるのは，「基礎データ」の中の「Ａ 人口・世帯」「Ｂ 自然環境」「Ｃ 経済基盤」のみですので，それら統計表に関するレコードのみにリストを絞り込みます．
#基礎データの統計表リスト
tbl_info_muni_kiso <- tbl_info_muni %>%
#STATISTICS_NAMEに「基礎」を含むレコードを抽出
dplyr::filter(grepl("基礎",STATISTICS_NAME)) %>%
#TITLEが「Ａ 人口・世帯」「Ｂ 自然環境」「Ｃ 経済基盤」のどれかに一致するレコードのみ残す
dplyr::filter(TITLE%in%c("Ａ 人口・世帯","Ｂ 自然環境","Ｃ 経済基盤"))
head(tbl_info_muni_kiso)## # A tibble: 3 × 19
## `@id` STAT_…¹ GOV_ORG STATI…² TITLE CYCLE SURVE…³ OPEN_…⁴ SMALL…⁵ COLLE…⁶
## <chr> <chr> <chr> <chr> <chr> <chr> <chr> <chr> <chr> <chr>
## 1 0000020201 社会・… 総務省 市区町… Ａ … 年度… 0 2022-0… 0 市区町…
## 2 0000020202 社会・… 総務省 市区町… Ｂ … 年度… 0 2022-0… 0 市区町…
## 3 0000020203 社会・… 総務省 市区町… Ｃ … 年度… 0 2022-0… 0 市区町…
## # … with 9 more variables: MAIN_CATEGORY <chr>, SUB_CATEGORY <chr>,
## # OVERALL_TOTAL_NUMBER <chr>, UPDATED_DATE <chr>, TABULATION_CATEGORY <chr>,
## # TABULATION_SUB_CATEGORY1 <chr>, TABULATION_CATEGORY_EXPLANATION <chr>,
## # TABLE_NAME <chr>, TABULATION_SUB_CATEGORY_EXPLANATION1 <chr>, and
## # abbreviated variable names ¹STAT_NAME, ²STATISTICS_NAME, ³SURVEY_DATE,
## # ⁴OPEN_DATE, ⁵SMALL_AREA, ⁶COLLECT_AREAestatapiのestat_getStatsData関数でリスト内にある統計表をダウンロードし，行方向に結合します．estat_getStatsData関数の引数statsDataIdには，リスト内の変数@idを指定します．まず，lapply関数を用いて，並列処理で各@idに対応する統計表をダウンロードした上でlistの各要素として格納します．その結果をdo.call関数に渡し，行結合の関数rbindを適用することで，要素を行方向に結合します．
#tbl_info_muni_kisoの変数@idをベクトルとして取り出し，ベクトルの各要素について以下の処理を実行
table_kiso <- lapply(X=tbl_info_muni_kiso$`@id`,FUN=function(x){
#ベクトルの[x]番目の要素（すなわち@id）に対応する統計表を取得
dat <- estatapi::estat_getStatsData(appId=appId,statsDataId=x) %>%
#不要な変数を落とす
dplyr::select(-c(tab_code,観測値,cat01_code,annotation)) %>%
#変数名を強制的に変更
magrittr::set_colnames(c("table_id","area_code","area_name","time_code","year","unit","value"))
#取得結果を返す
return(dat)
}) %>%
#listの要素を行方向に結合
dplyr::bind_rows()
head(table_kiso)## # A tibble: 6 × 7
## table_id area_code area_name time_code year unit value
## <chr> <chr> <chr> <chr> <chr> <chr> <dbl>
## 1 A1101_総人口 01100 北海道 札幌市 1980100000 1980年度 人 1401757
## 2 A1101_総人口 01100 北海道 札幌市 1985100000 1985年度 人 1542979
## 3 A1101_総人口 01100 北海道 札幌市 1990100000 1990年度 人 1671742
## 4 A1101_総人口 01100 北海道 札幌市 1995100000 1995年度 人 1757025
## 5 A1101_総人口 01100 北海道 札幌市 2000100000 2000年度 人 1822368
## 6 A1101_総人口 01100 北海道 札幌市 2005100000 2005年度 人 1880863パネルデータ構築に用いる変数を，2000年度，2005年度，2010年度の3時点で抽出します．
kiso_vars <- table_kiso %>%
#必要な変数に該当するレコードのみ抽出
dplyr::filter(table_id=="A1101_総人口"|
table_id=="B1101_総面積（北方地域及び竹島を除く）"|
table_id=="C120110_課税対象所得"|
table_id=="C120120_納税義務者数（所得割）") %>%
#yearが「2000年度」「2005年度」「2010年度」のどれかのレコードのみ抽出
dplyr::filter(year%in%c("2000年度","2005年度","2010年度"))
head(kiso_vars)## # A tibble: 6 × 7
## table_id area_code area_name time_code year unit value
## <chr> <chr> <chr> <chr> <chr> <chr> <dbl>
## 1 A1101_総人口 01100 北海道 札幌市 2000100000 2000年度 人 1822368
## 2 A1101_総人口 01100 北海道 札幌市 2005100000 2005年度 人 1880863
## 3 A1101_総人口 01100 北海道 札幌市 2010100000 2010年度 人 1913545
## 4 A1101_総人口 01101 北海道 札幌市 中央区 2000100000 2000年度 人 181383
## 5 A1101_総人口 01101 北海道 札幌市 中央区 2005100000 2005年度 人 202801
## 6 A1101_総人口 01101 北海道 札幌市 中央区 2010100000 2010年度 人 220189現状取得されたデータはlong型で格納されています．すなわち，全変数の観測値が行方向に並んでいるような形式です．これを，tidyrのpivot_wider関数を用いて，市区町村毎に各変数の観測値が横並びになるwide型に変形します．どの変数を基準にデータを横並びにするかを指定するのは引数names_fromで，今回は各変数の名前table_idで指定します．values_fromには，横並びにする値が入った列を指定します．今回のケースではvalueです．
kiso_vars_wide <- kiso_vars %>%
#不要な変数を落とす
dplyr::select(-c(time_code,unit)) %>%
#データをwide型に変換
tidyr::pivot_wider(names_from=table_id,
values_from=value)
head(kiso_vars_wide)## # A tibble: 6 × 7
## area_code area_name year A1101_総人口 B1101_…¹ C12011…² C1201…³
## <chr> <chr> <chr> <dbl> <dbl> <dbl> <dbl>
## 1 01100 北海道 札幌市 2000年度 1822368 112112 2.59e9 736662
## 2 01100 北海道 札幌市 2005年度 1880863 112112 2.45e9 749998
## 3 01100 北海道 札幌市 2010年度 1913545 112112 2.41e9 797193
## 4 01101 北海道 札幌市 中央区 2000年度 181383 4642 NA NA
## 5 01101 北海道 札幌市 中央区 2005年度 202801 4642 NA NA
## 6 01101 北海道 札幌市 中央区 2010年度 220189 4642 NA NA
## # … with abbreviated variable names ¹`B1101_総面積（北方地域及び竹島を除く）`,
## # ²C120110_課税対象所得, ³`C120120_納税義務者数（所得割）`パネルデータ内の変数を作成します．まず，総面積あたり人口密度densと，納税義務者数あたり課税対象所得incomeを新たな変数として追加します．その後，総面積と総人口について，それぞれareaとpopに変数名を変更します．欠損値を含まないレコードのみを残した上で，以後用いない変数を落とします．
kiso_vars_wide <- kiso_vars_wide %>%
#変数を計算
dplyr::mutate(dens=A1101_総人口/`B1101_総面積（北方地域及び竹島を除く）`,
income=C120110_課税対象所得/`C120120_納税義務者数（所得割）`) %>%
#変数の名前を変更
dplyr::rename(area=`B1101_総面積（北方地域及び竹島を除く）`,
pop=A1101_総人口) %>%
#incomeとdens両方について欠損がないレコードのみ残す
dplyr::filter(!is.na(income)&!is.na(dens)) %>%
#不要な変数を削除
dplyr::select(-c("C120110_課税対象所得","C120120_納税義務者数（所得割）"))
head(kiso_vars_wide)## # A tibble: 6 × 7
## area_code area_name year pop area dens income
## <chr> <chr> <chr> <dbl> <dbl> <dbl> <dbl>
## 1 01100 北海道 札幌市 2000年度 1822368 112112 16.3 3510.
## 2 01100 北海道 札幌市 2005年度 1880863 112112 16.8 3272.
## 3 01100 北海道 札幌市 2010年度 1913545 112112 17.1 3017.
## 4 01202 北海道 函館市 2000年度 305311 67748 4.51 3197.
## 5 01202 北海道 函館市 2005年度 294264 67782 4.34 2981.
## 6 01202 北海道 函館市 2010年度 279127 67793 4.12 2715.データをバランス化するため，観測値が3時点で変数が完全には揃わない市区町村を調べます．その上で，該当する市区町村コード（変数area_code）をベクトルとして抽出します．
area_code_na <- kiso_vars_wide %>%
#市区町村コードでグループ化
dplyr::group_by(area_code) %>%
#市区町村毎のレコード数を集計
dplyr::summarise(count=n()) %>%
dplyr::ungroup() %>%
#レコード数が3に満たない（3時点で変数が完全には揃わない）市区町村コードを取得
dplyr::filter(count<3) %>%
#area_codeをベクトルとして取得
dplyr::pull(area_code)先ほどwide型に変換したデータについて，変数area_codeがベクトルarea_code_na（欠測値がある市区町村のコード）の要素に一致しないレコードを除きます．これにより，データをバランスしたパネルデータにできます．加えて，area_codeの先頭を0埋めした変数JCODEや，その自治体が市区部か町村部かを判別するための変数municipality_classを作成します．最後に，変数yearを数値として扱えるように変換します．
kiso_vars_wide_balance <- kiso_vars_wide %>%
#上で取得した欠損値有りのarea_codeに該当する市区町村のレコードを除く
dplyr::filter(!area_code%in%area_code_na) %>%
#area_codeの先頭を0埋めしたものを変数JCODEとして追加
dplyr::mutate(JCODE=formatC(x=area_code,width=5,flag="0")) %>%
#area_codeを削除
dplyr::select(-area_code) %>%
#JCODEを先頭に移動
dplyr::relocate(JCODE) %>%
#市区・町村区分変数を作成
dplyr::mutate(municipality_class=dplyr::case_when(
substr(area_name,nchar(area_name),nchar(area_name))%in%c("市","区") ~ "市区",
substr(area_name,nchar(area_name),nchar(area_name))%in%c("町","村") ~ "町村"
)) %>%
#時点変数を作成
dplyr::mutate(year=dplyr::case_when(
year=="2000年度" ~ 2000,
year=="2005年度" ~ 2005,
year=="2010年度" ~ 2010
))データには，東京都特別区部全体に関するレコードも入ってしまっているので，それを除きます．
kiso_vars_wide_balance <- kiso_vars_wide_balance %>%
dplyr::filter(JCODE!="13100")以後の作業で使わないオブジェクトを削除します．
rm(list=ls()[ls()!="kiso_vars_wide_balance"])
gc();gc()## used (Mb) gc trigger (Mb) limit (Mb) max used (Mb)
## Ncells 1977867 105.7 3625485 193.7 NA 3625485 193.7
## Vcells 3386680 25.9 65428045 499.2 16384 80549412 614.6## used (Mb) gc trigger (Mb) limit (Mb) max used (Mb)
## Ncells 1978002 105.7 3625485 193.7 NA 3625485 193.7
## Vcells 3387038 25.9 52342436 399.4 16384 80549412 614.6市区町村ポリゴンをTopoJSON形式で取得します．歴史的行政区域データセットの場合，ファイルはURLを指定することで直接取得が可能です．TopoJSONをsf形式で読み込む場合，geojsonioのtopojson_read関数を用います．市区町村ポリゴンには投影法が定義されていないので，sfのst_set_crs関数で投影法を定義します．その後，st_transform関数を用いて投影座標系に変換を行います．変数JCODEが欠損しているレコードを除いた上で，st_make_valid関数で不整地物の修正を行います．不整が含まれていると，空間データの各種処理の際にエラーが出るので，可能な限り修正を行うことをお勧めします．
##市区町村ポリゴンの読み込み
poly <- geojsonio::topojson_read(x="https://geoshape.ex.nii.ac.jp/city/topojson/20210101/jp_city_dc.c.topojson") %>%
#投影法をWGS84で定義
sf::st_set_crs(value=sf::st_crs(4326)) %>%
#WGS84/UTM54Nに投影変換
sf::st_transform(crs=sf::st_crs(32654)) %>%
#JCODEがNAのレコードを削除
dplyr::filter(!is.na(N03_007)) %>%
#不整地物の修正
sf::st_make_valid()演習で市区町村ポリゴンを用いるので，write_sf関数を用いてGeoJSON形式で書き出します．
sf::write_sf(obj=poly,
dsn="ab_jp_dc.geojson",
#既に同名のファイルがあれば削除
delete_dsn=TRUE)国土数値情報からのshapefileのダウンロードと，ダウンロードされたデータのRへの読み込み作業を，Webスクレイピングを用いて行う方法を示します．手続きが少々煩雑なので，工程をふたつの関数に分けた上で，各データについて適用します．
以下で示す手順を用いれば，手動でzipファイルをダウンロード・解凍したり，それらをローカル環境から読み込んだりする手間がかかりません．
まず，get_zip_links関数を定義します．
get_zip_links <- function(url,root="https://nlftp.mlit.go.jp/ksj/gml/"){
#urlのHTMLファイルを取得
all_attriv <- xml2::read_html(url)
#zipファイルへのリンクがある場所を探す
onclick <- all_attriv %>%
#aタグに絞る
rvest::html_nodes("a") %>%
#onclick属性に絞る
rvest::html_attr("onclick") %>%
#スペースを削除
gsub(" ","",.)
#zipファイルへのリンクを含む箇所を抽出
zip_links <- sapply(onclick[grepl(".zip",onclick)] %>% stringr::str_split("','"),function(x){x[3]}) %>%
#いらない文字列を削除
gsub("',this);","",.) %>%
gsub("\\../","",.) %>%
#RootのURLと文字列結合
paste0(root,.)
#作成されたURLを返す
return(zip_links)
}次に，sf_from_zip関数を定義します．
sf_from_zip <- function(links){
#linksの各要素について以下の処理
dat_out <- lapply(X=links,FUN=function(x){
#一時ファイルを作成
temp <- tempfile()
#一時ファイル内にzipをダウンロード
download.file(x,temp)
#zip展開先の一時ファイルを作成
temp_zip <- tempfile()
#一時ファイルにzipを解凍
unzipped <- unzip(temp,exdir=temp_zip)
#shpの名前を特定した上で，shpを読み込み
dat <- sf::read_sf(dsn=unzipped[grepl(".shp",unzipped)])
#使い終わったファイルの削除
unlink(temp);rm(temp)
unlink(temp_zip);rm(temp_zip)
#読み込んだデータを返す
return(dat)
}) %>%
#listの要素を行方向に結合
dplyr::bind_rows()
#結合されたデータを返す
return(dat_out)
}上で定義した関数を用いて，標高・傾斜データ（ポイント形式）の読み込みを行います．
slope <- sf_from_zip(links=get_zip_links(url="https://nlftp.mlit.go.jp/ksj/gml/datalist/KsjTmplt-G04-d.html"))
head(slope)## Simple feature collection with 6 features and 10 fields
## Geometry type: POLYGON
## Dimension: XY
## Bounding box: xmin: 136.0625 ymin: 20.425 xmax: 136.075 ymax: 20.42917
## CRS: NA
## # A tibble: 6 × 11
## G04d_001 G04d_002 G04d_003 G04d_004 G04d_005 G04d_…¹ G04d_…² G04d_…³ G04d_…⁴
## <chr> <chr> <chr> <chr> <chr> <chr> <chr> <chr> <chr>
## 1 3036501511 unknown unknown unknown unknown unknown unknown unknown unknown
## 2 3036501512 unknown unknown unknown unknown unknown unknown unknown unknown
## 3 3036501521 1.0 1.0 1.0 0 0.0 0 0.0 0
## 4 3036501522 unknown unknown unknown unknown unknown unknown unknown unknown
## 5 3036501513 unknown unknown unknown unknown unknown unknown unknown unknown
## 6 3036501514 unknown unknown unknown unknown unknown unknown unknown unknown
## # … with 2 more variables: G04d_010 <chr>, geometry <POLYGON>, and abbreviated
## # variable names ¹G04d_006, ²G04d_007, ³G04d_008, ⁴G04d_009同様にして，海岸線データ（ライン形式）の読み込みを行います．
coast <- sf_from_zip(links=get_zip_links(url="https://nlftp.mlit.go.jp/ksj/gml/datalist/KsjTmplt-C23.html"))
head(coast)## Simple feature collection with 6 features and 7 fields
## Geometry type: LINESTRING
## Dimension: XY
## Bounding box: xmin: 140.6925 ymin: 41.71013 xmax: 141.1874 ymax: 42.00935
## CRS: NA
## # A tibble: 6 × 8
## C23_001 C23_002 C23_003 C23_004 C23_005 C23_006 C23_007
## <chr> <chr> <chr> <chr> <chr> <chr> <chr>
## 1 01202 6 <NA> <NA> <NA> <NA> false
## 2 01202 4 9999 臼尻漁港 臼尻西 9 不明 false
## 3 01202 6 <NA> <NA> <NA> <NA> false
## 4 01202 6 <NA> <NA> <NA> <NA> false
## 5 01202 6 <NA> <NA> <NA> <NA> false
## 6 01202 6 <NA> <NA> <NA> <NA> false
## # … with 1 more variable: geometry <LINESTRING>手動でzipファイルのダウンロードと解凍を行った場合には，以下のコードを実行することで，上で示したものと同じファイルが作成できます．なお，それらデータはサブフォルダから取り出され，このRmdと同じ階層のディレクトリ（もしくは自分で指定した作業ディレクトリ）にあるものとします．なお，Rmd実行の都合上，コードはコメント化しています．
最初に標高・傾斜ポイントのshapefileに該当するファイル名を取得し，lapply関数を用いて並列処理でファイルを読み込むというふたつのステップで処理を実行します．
# #ディレクトリ内にある全ファイルのリスト
# slope_shp_list <- list.files() %>%
# #そのうち名前に「ElevationAndSlopeAngleFifthMesh.shp」が含まれるもの
# .[grepl("ElevationAndSlopeAngleFifthMesh.shp",.)]
# #各データを読み込んで行方向に結合
# slope <- lapply(slope_shp_list,function(x){
# dat <- sf::read_sf(dsn=x)
# return(dat)
# }) %>%
# #listの要素を行方向に結合
# dplyr::bind_rows()標高・傾斜データと同じ手順で処理を実行します．
# #ディレクトリ内にある全ファイルのリスト
# coast_shp_list <- list.files() %>%
# #そのうち名前に「Coastline.shp」が含まれるもの
# .[grepl("Coastline.shp",.)]
# #各データを読み込んで行方向に結合
# coast <- lapply(coast_shp_list,function(x){
# dat <- sf::read_sf(dsn=x)
# return(dat)
# }) %>%
# #listの要素を行方向に結合
# dplyr::bind_rows()各種空間分析が可能なように，標高・傾斜ポイントの前処理を行います．詳細はコード内のコメントを参照してください．
slope <- slope %>%
#座標系がJGD2000であると指定
sf::st_set_crs(sf::st_crs(4612)) %>%
#WGS84/UTM54Nに投影変換
sf::st_transform(sf::st_crs(32654)) %>%
#重心座標を取得
sf::st_centroid() %>%
#必要な変数だけ残す
dplyr::select(G04d_002,G04d_010) %>%
#変数名をつけ直す
dplyr::rename(altitude=G04d_002,slope=G04d_010) %>%
#平均標高が欠損のレコードを除去
dplyr::filter(altitude!="unknown") %>%
#平均傾斜度が欠損のレコードを除去
dplyr::filter(slope!="unknown") %>%
#文字列型になってしまっている変数を数値型に変換
dplyr::mutate(across(is.character,as.numeric))市区町村毎に，標高・傾斜に関する変数を作成します．標高変数は市区町村内の平均標高，傾斜変数は傾斜度の合計で定義されます．まずはsfのst_intersection関数を用いて市区町村ポリゴンと標高・傾斜ポイントの共通部分を取り，各ポイントに市区町村コード（N03_007）を空間結合します．結合結果にst_drop_geometry関数を適用してgeometry属性を削除した上で，データを市区町村コードでグループ化し，平均標高・傾斜度の合計を各々計算します．最後に市区町村コードの変数名をN03_007からJCODEに変更します．
poly_slope <- poly %>%
#市区町村ポリゴンと標高・傾斜ポイントの共通部分を取る
sf::st_intersection(y=slope) %>%
#geometry属性を削除
sf::st_drop_geometry() %>%
#JCODEでグループ化
dplyr::group_by(N03_007) %>%
#標高の平均値を計算
dplyr::summarise(altitude=mean(altitude),
#傾斜度の合計を計算
slope_sum=sum(slope)) %>%
#グループ化の解除
dplyr::ungroup() %>%
#変数名を変更
dplyr::rename(JCODE=N03_007)
#不要なオブジェクトの削除
rm(list=c("slope"))
gc();gc()## used (Mb) gc trigger (Mb) limit (Mb) max used (Mb)
## Ncells 2466864 131.8 89835929 4797.8 NA 140368638 7496.5
## Vcells 17381449 132.7 526525084 4017.1 16384 658156290 5021.4## used (Mb) gc trigger (Mb) limit (Mb) max used (Mb)
## Ncells 2466870 131.8 71868744 3838.3 NA 140368638 7496.5
## Vcells 17381537 132.7 421220068 3213.7 16384 658156290 5021.4各種空間分析が可能なように，海岸線ラインの前処理を行います．詳細はコード内のコメントを参照してください．
coast <- coast %>%
#座標系がJGD2000であると指定
sf::st_set_crs(sf::st_crs(4612)) %>%
#WGS84/UTM54Nに投影変換
sf::st_transform(sf::st_crs(32654))市区町村ポリゴンの重心ポイントから最寄りの海岸線ラインまでの距離を計算します．ある地物から最寄りの別の地物を探索する際には，nngeoのst_nn関数を用います．
poly_cent_coast <- sf::st_centroid(poly) %>%
#重心ポイントから最寄り海岸線ラインまでの距離を計算し，変数として追加
dplyr::mutate(d_coast=unlist(nngeo::st_nn(.,coast,returnDist=TRUE)[[2]])) %>%
#geometry属性の削除
sf::st_drop_geometry() %>%
#JCODEと海岸線距離のみを残す
dplyr::select(N03_007,d_coast) %>%
#変数名を変更
dplyr::rename(JCODE=N03_007)##
|
| | 0%
|
| | 1%
|
|= | 1%
|
|= | 2%
|
|== | 2%
|
|== | 3%
|
|== | 4%
|
|=== | 4%
|
|=== | 5%
|
|==== | 5%
|
|==== | 6%
|
|===== | 6%
|
|===== | 7%
|
|===== | 8%
|
|====== | 8%
|
|====== | 9%
|
|======= | 9%
|
|======= | 10%
|
|======= | 11%
|
|======== | 11%
|
|======== | 12%
|
|========= | 12%
|
|========= | 13%
|
|========= | 14%
|
|========== | 14%
|
|========== | 15%
|
|=========== | 15%
|
|=========== | 16%
|
|============ | 16%
|
|============ | 17%
|
|============ | 18%
|
|============= | 18%
|
|============= | 19%
|
|============== | 19%
|
|============== | 20%
|
|============== | 21%
|
|=============== | 21%
|
|=============== | 22%
|
|================ | 22%
|
|================ | 23%
|
|================ | 24%
|
|================= | 24%
|
|================= | 25%
|
|================== | 25%
|
|================== | 26%
|
|=================== | 26%
|
|=================== | 27%
|
|=================== | 28%
|
|==================== | 28%
|
|==================== | 29%
|
|===================== | 29%
|
|===================== | 30%
|
|===================== | 31%
|
|====================== | 31%
|
|====================== | 32%
|
|======================= | 32%
|
|======================= | 33%
|
|======================= | 34%
|
|======================== | 34%
|
|======================== | 35%
|
|========================= | 35%
|
|========================= | 36%
|
|========================== | 36%
|
|========================== | 37%
|
|========================== | 38%
|
|=========================== | 38%
|
|=========================== | 39%
|
|============================ | 39%
|
|============================ | 40%
|
|============================ | 41%
|
|============================= | 41%
|
|============================= | 42%
|
|============================== | 42%
|
|============================== | 43%
|
|============================== | 44%
|
|=============================== | 44%
|
|=============================== | 45%
|
|================================ | 45%
|
|================================ | 46%
|
|================================= | 46%
|
|================================= | 47%
|
|================================= | 48%
|
|================================== | 48%
|
|================================== | 49%
|
|=================================== | 49%
|
|=================================== | 50%
|
|=================================== | 51%
|
|==================================== | 51%
|
|==================================== | 52%
|
|===================================== | 52%
|
|===================================== | 53%
|
|===================================== | 54%
|
|====================================== | 54%
|
|====================================== | 55%
|
|======================================= | 55%
|
|======================================= | 56%
|
|======================================== | 56%
|
|======================================== | 57%
|
|======================================== | 58%
|
|========================================= | 58%
|
|========================================= | 59%
|
|========================================== | 59%
|
|========================================== | 60%
|
|========================================== | 61%
|
|=========================================== | 61%
|
|=========================================== | 62%
|
|============================================ | 62%
|
|============================================ | 63%
|
|============================================ | 64%
|
|============================================= | 64%
|
|============================================= | 65%
|
|============================================== | 65%
|
|============================================== | 66%
|
|=============================================== | 66%
|
|=============================================== | 67%
|
|=============================================== | 68%
|
|================================================ | 68%
|
|================================================ | 69%
|
|================================================= | 69%
|
|================================================= | 70%
|
|================================================= | 71%
|
|================================================== | 71%
|
|================================================== | 72%
|
|=================================================== | 72%
|
|=================================================== | 73%
|
|=================================================== | 74%
|
|==================================================== | 74%
|
|==================================================== | 75%
|
|===================================================== | 75%
|
|===================================================== | 76%
|
|====================================================== | 76%
|
|====================================================== | 77%
|
|====================================================== | 78%
|
|======================================================= | 78%
|
|======================================================= | 79%
|
|======================================================== | 79%
|
|======================================================== | 80%
|
|======================================================== | 81%
|
|========================================================= | 81%
|
|========================================================= | 82%
|
|========================================================== | 82%
|
|========================================================== | 83%
|
|========================================================== | 84%
|
|=========================================================== | 84%
|
|=========================================================== | 85%
|
|============================================================ | 85%
|
|============================================================ | 86%
|
|============================================================= | 86%
|
|============================================================= | 87%
|
|============================================================= | 88%
|
|============================================================== | 88%
|
|============================================================== | 89%
|
|=============================================================== | 89%
|
|=============================================================== | 90%
|
|=============================================================== | 91%
|
|================================================================ | 91%
|
|================================================================ | 92%
|
|================================================================= | 92%
|
|================================================================= | 93%
|
|================================================================= | 94%
|
|================================================================== | 94%
|
|================================================================== | 95%
|
|=================================================================== | 95%
|
|=================================================================== | 96%
|
|==================================================================== | 96%
|
|==================================================================== | 97%
|
|==================================================================== | 98%
|
|===================================================================== | 98%
|
|===================================================================== | 99%
|
|======================================================================| 99%
|
|======================================================================| 100%#不要なオブジェクトの削除
rm(list=c("coast"))
gc();gc()## used (Mb) gc trigger (Mb) limit (Mb) max used (Mb)
## Ncells 2100423 112.2 57494996 3070.6 NA 140368638 7496.5
## Vcells 8097507 61.8 336976055 2571.0 16384 658156290 5021.4## used (Mb) gc trigger (Mb) limit (Mb) max used (Mb)
## Ncells 2096217 112.0 45995997 2456.5 NA 140368638 7496.5
## Vcells 8090575 61.8 269580844 2056.8 16384 658156290 5021.4ここまで作成してきた各種データをひとつにまとめます．それに先立ち，所得・人口パネルデータ内で各種変数を追加します．なお，三宅島噴火によって，当該期間に人口が0となっていた東京都三宅村のレコードをこの段階で除外します．
income_dens_panel <- kiso_vars_wide_balance %>%
#所得・人口密度の対数値を変数として追加
dplyr::mutate(ln_income=log(income),
ln_dens=log(dens)) %>%
#三宅村のレコードを除外
dplyr::filter(!JCODE%in%c("13381")) %>%
#総面積の対数値を計算
dplyr::mutate(ln_area=log(area))パネルデータにコントロール変数のデータを結合します．
income_dens_panel <- income_dens_panel %>%
dplyr::left_join(y=poly_slope,by="JCODE") %>%
dplyr::left_join(y=poly_cent_coast,by="JCODE") %>%
#海岸線までの距離の対数を取る
dplyr::mutate(ln_d_coast=log(d_coast))統合されたデータを，openxlsxのwrite.xlsx関数を用いて書き出します．
openxlsx::write.xlsx(x=income_dens_panel,
file="income_dens_panel.xlsx")