もしもGeoJSONをR上で扱ったら：国土交通データプラットフォームの事例
Keisuke Takano
4/25/2020
Introduction
2020年4月24日に，国土交通省による国土交通データプラットフォームの提供が開始されました．公式WEBサイトによると，このプラットフォームでは，国・地方自治体の保有する橋梁やトンネル、ダムや水門などの社会インフラ（施設）の諸元や点検結果に関するデータ約8万件と全国のボーリング結果等の地盤データ約14万件の計22万件を地図上に表示可能です．
このプラットフォームのありがたい点は，地図上に表示される各種位置情報を，各種GISソフトで扱うことができるGeoJSON形式で提供していることです（GeoJSONの詳細はこちらのWEBサイトを参照）．
このサイトでは，提供されているGeoJSONの使い心地を少しみてみます．
なお，このサイトで使用されている各種関数・パッケージの使い方に関しては，別途チュートリアルサイトを用意してあるので，そちらを参照ください．
Updates
- 橋梁ポイントの座標系と市区町村境界ポリゴンの座標系が異なる問題を修正しました（2020年4月26日）．
- 茨城県外にプロットされてしまっている橋梁ポイントを抽出するコードを追加しました（2020年4月26日）．
- パッケージのアップデートに対応しました（2021年1月2日）．
Load Packages
今回使うパッケージを全て読み込みます．
library(rvest)
library(dplyr)
library(geojsonsf)
library(stringr)
library(magrittr)
library(ggplot2)Download & Import GeoJSON Datasets
最初に，今回のメインデータセットであるインフラデータをダウンロード・読み込みします．今回は，一連のデータセットから，橋梁GeoJSONデータをダウンロード・読み込みします．
ダウンロードページでは，1県に対して1つのGeoJSONのリンクが置いてあります．しかし，リンクアドレスには都道府県コードなどに基づく規則性がないので，for文を使ってリンクアドレスを変えながら全県のGeoJSONを一括ダウンロードができません．
そこで，WEBスクレイピングのテクニックを使って，1県ずつリンクをクリック＆保存という煩雑な作業を避けつつ，一括ダウンロードを行う方法を用います（RによるWEBスクレイピングに関しては，このWEBサイトを参照）．
rvestパッケージの関数が適用されている部分が，WEBスクレイピングを行っている箇所になります．
各県のGeoJSONデータのダウンロードが済んだら，全県のレコードを縦結合し，後程使う「点検の結果_総合」変数がNAでないレコードのみを残します．
#「社会資本情報プラットフォーム（道路－橋梁）」ページのソースを読み込み
page <- xml2::read_html("https://www.geospatial.jp/ckan/dataset/ipf0101")
#ソースから，リンクを抽出
json_url <- page %>% rvest::html_nodes("a") %>%
#リンクのURLを抽出
rvest::html_attr("href") %>%
#URLが".geojson"で終わるものを抽出
stringr::str_subset("\\.geojson")
#上で作成したリスト内の全てのgeojsonをダウンロードし，縦結合
bridge <- do.call(rbind,lapply(json_url,function(x){geojsonsf::geojson_sf(x)})) %>%
magrittr::set_colnames(gsub(" ","_",colnames(.)))
#点検の結果がNAになっているレコードを削除した上で，ポイントデータに変換
bridge_noNA <- bridge %>% dplyr::filter(!is.na(点検の結果_総合)) %>%
sf::st_as_sf(sf_column_name="geometry") %>%
sf::st_transform(sf::st_crs(4612))
head(bridge_noNA)## Simple feature collection with 6 features and 24 fields
## geometry type: POINT
## dimension: XY
## bbox: xmin: 143.1011 ymin: 42.01722 xmax: 143.3161 ymax: 42.12472
## geographic CRS: JGD2000
## 分野名称_1 橋長 管理者名称_1 分野名称_4 幅員 分野名称_2
## 1 道路 24 国の機関-国土交通省北海道開発局 橋梁 8.5 -
## 2 道路 178.9 国の機関-国土交通省北海道開発局 橋梁 6 -
## 3 道路 8.4 国の機関-国土交通省北海道開発局 橋梁 8.5 -
## 4 道路 19.8 国の機関-国土交通省北海道開発局 橋梁 8.5 -
## 5 道路 35 国の機関-国土交通省北海道開発局 橋梁 9.5 -
## 6 道路 6.3 国の機関-国土交通省北海道開発局 橋梁 10.5 -
## 分野名称_3 管理者名称_2_code 分野名称_4_code 分野名称_1_code 管理者名称_2
## 1 国道336号 - D010 04 -
## 2 国道336号 - D010 04 -
## 3 国道336号 - D010 04 -
## 4 国道336号 - D010 04 -
## 5 国道336号 - D010 04 -
## 6 国道336号 - D010 04 -
## 管理者名称_3 都道府県_code 管理者名称_4 都道府県 管理者名称_1_code
## 1 - 01 - 北海道 022
## 2 - 01 - 北海道 022
## 3 - 01 - 北海道 022
## 4 - 01 - 北海道 022
## 5 - 01 - 北海道 022
## 6 - 01 - 北海道 022
## geometry 点検の結果_総合 位置_経度 施設所在地1 施設名称
## 1 POINT (143.2131 42.02583) Ⅱ 143.2131 えりも町 伏見橋
## 2 POINT (143.3161 42.12472) Ⅱ 143.3161 えりも町 猿留橋
## 3 POINT (143.1469 42.01722) Ⅰ 143.1469 えりも町 幌泉橋
## 4 POINT (143.2631 42.02306) Ⅰ 143.2631 えりも町 阿津橋
## 5 POINT (143.3069 42.075) Ⅰ 143.3069 えりも町 新咲梅橋
## 6 POINT (143.1011 42.04667) Ⅰ 143.1011 えりも町 笛舞橋
## 位置_緯度 時期_完成もしくは供用開始年 管理者名称_3_code 点検実施年度
## 1 42.02583 1972 - H30年度
## 2 42.12472 1960 - H29年度
## 3 42.01722 1965 - H30年度
## 4 42.02306 1973 - H30年度
## 5 42.07500 2004 - H28年度
## 6 42.04667 1972 - H27年度# #[NOTRUN]緯度経度変数を使ってもポイントデータに変換可能
# bridge_noNA <- bridge %>% dplyr::filter(!is.na(点検の結果_総合)) %>%
# sf::st_as_sf(coords=c("位置_経度","位置_緯度"),crs=sf::st_crs(4612))Download & Import Municipal Boundary Datasets
市区町村境界のshapefileをダウンロード・読み込みします．手続きに関しては，チュートリアルサイトに記したものと全く同じです．
読み込み後，市区町村境界を都道府県ごとに融合し，都道府県境界ポリゴンを作成します．その際，都道府県名を，後程都道府県境界ポリゴン内の変数として用いたいので，st_union関数ではなく，summarise関数に基づきポリゴンを融合しています．
#一時ディレクトリを作成
boundary <- tempfile()
#shapefileが入っているzipファイルを一時ディレクトリにダウンロード
curl::curl_download(url="https://www.esrij.com/cgi-bin/wp/wp-content/uploads/2017/01/japan_ver81.zip",destfile=boundary)
#zipファイルを解凍
unzip(zipfile=boundary,exdir=getwd())
#一時ディレクトリの削除
unlink(boundary);rm(boundary)
#市区町村境界shapefileの読み込み（一部不正なポリゴンがあったので，st_make_validで整形）
bound <- sf::st_read(dsn=getwd(),layer="japan_ver81",stringsAsFactors=F,options="ENCODING=CP932",quiet=TRUE) %>%
sf::st_make_valid(.)
head(bound)## Simple feature collection with 6 features and 9 fields
## geometry type: MULTIPOLYGON
## dimension: XY
## bbox: xmin: 140.9905 ymin: 42.78073 xmax: 141.4735 ymax: 43.18944
## geographic CRS: JGD2000
## JCODE KEN SICHO GUN SEIREI SIKUCHOSON CITY_ENG
## 1 01101 北海道 石狩振興局 <NA> 札幌市 中央区 Sapporo-shi, Chuo-ku
## 2 01102 北海道 石狩振興局 <NA> 札幌市 北区 Sapporo-shi, Kita-ku
## 3 01103 北海道 石狩振興局 <NA> 札幌市 東区 Sapporo-shi, Higashi-ku
## 4 01104 北海道 石狩振興局 <NA> 札幌市 白石区 Sapporo-shi, Shiroishi-ku
## 5 01105 北海道 石狩振興局 <NA> 札幌市 豊平区 Sapporo-shi, Toyohira-ku
## 6 01106 北海道 石狩振興局 <NA> 札幌市 南区 Sapporo-shi, Minami-ku
## P_NUM H_NUM geometry
## 1 230382 136592 MULTIPOLYGON (((141.335 43....
## 2 283392 147700 MULTIPOLYGON (((141.3982 43...
## 3 259826 138498 MULTIPOLYGON (((141.4472 43...
## 4 210479 117531 MULTIPOLYGON (((141.4735 43...
## 5 219071 121796 MULTIPOLYGON (((141.3663 43...
## 6 140809 72126 MULTIPOLYGON (((141.0964 43...#都道府県単位にポリゴンを融合
bound_pref <- bound %>% dplyr::group_by(KEN) %>% dplyr::summarize()Plot of Bridge Point Data in Ibaraki
各ポイントが，橋梁の位置を表しています．また，「点検の結果_総合」変数に基づいてポイントの色分けを行っています．
どういうわけか，茨城県外の場所にもポイントがプロットされているので，元データの座標値に誤りがあると思われます（GeoJSONに付与されているgeometry変数を用いても，緯度経度座標変数を用いても同じ結果でした）．
#市区町村境界ポリゴンから，茨城県に含まれるもののみ抽出
bound_ibaraki <- bound %>% dplyr::filter(KEN=="茨城県")
#橋梁ポイントから茨城県に該当するもののみ抽出
bridge_noNA_ibaraki <- bridge_noNA %>% dplyr::filter(都道府県=="茨城県")
#まずは市区町村境界ポリゴン（白地図）をプロットし，その上に橋梁ポイント（点検の結果カテゴリに基づき着色）を重ね描きする
ggplot2::ggplot() +
ggplot2::geom_sf(data=bound_ibaraki,fill="white") +
ggplot2::geom_sf(data=bridge_noNA_ibaraki,aes(color=点検の結果_総合)) +
ggplot2::theme_gray(base_family="HiraKakuPro-W3")
#なぜか茨城県外にプロットされてしまうポイントを特定する（もっと良い方法があるように思える）
bridge_noNA_ibaraki_ID <- bridge_noNA_ibaraki %>% dplyr::mutate(FID=dplyr::row_number())
bridge_noNA_ibaraki_ID_within_ibaraki <- sf::st_intersection(bridge_noNA_ibaraki_ID,bound_ibaraki)
bridge_noNA_ibaraki_ID %>% dplyr::filter(!FID%in%(bridge_noNA_ibaraki_ID_within_ibaraki$FID))## Simple feature collection with 2 features and 25 fields
## geometry type: POINT
## dimension: XY
## bbox: xmin: 140.0645 ymin: 36.90278 xmax: 141.3383 ymax: 37.28472
## geographic CRS: JGD2000
## 分野名称_1 橋長 管理者名称_1 分野名称_4 幅員 分野名称_2 分野名称_3
## 1 道路 11.6 国の機関-国土交通省 橋梁 1.5 - 国道6号
## 2 道路 13.5 国の機関-国土交通省 橋梁 1.8 - 国道6号
## 管理者名称_2_code 分野名称_4_code 分野名称_1_code 管理者名称_2 管理者名称_3
## 1 03 D010 04 関東地方整備局 -
## 2 03 D010 04 関東地方整備局 -
## 都道府県_code 管理者名称_4 都道府県 管理者名称_1_code
## 1 08 - 茨城県 021
## 2 08 - 茨城県 021
## geometry 点検の結果_総合 位置_経度 施設所在地1
## 1 POINT (141.3383 37.28472) Ⅱ 141.3383 北茨城市
## 2 POINT (140.0645 36.90278) Ⅱ 140.0645 取手市
## 施設名称 位置_緯度 時期_完成もしくは供用開始年
## 1 沢尻川側道橋 37.28472 2004
## 2 取手跨線橋側道橋（上） 36.90278 1984
## 管理者名称_3_code 点検実施年度 FID
## 1 - H26年度 24
## 2 - H30年度 26Plot Inspection Result by Prefecture
都道府県別に「点検の結果_総合」がⅢ（早期措置段階）かⅣ（緊急措置段階）である橋梁の割合を計算し，地図にプロットします．
#変数「点検の結果_総合」が"Ⅲ"もしくは"Ⅳ"のレコードで1をとり，さもなくば0をとる変数dangを作成
bridge_noNA_bin <- bridge_noNA %>%
dplyr::mutate(dang=ifelse(点検の結果_総合=="Ⅲ"|点検の結果_総合=="Ⅳ",1,0))
#変数dangが付与された橋梁ポイントからジオメトリを除去し，都道府県でグループ化した後に，変数dangの平均（＝dangが1になる割合）dang_ratioを計算
bridge_noNA_bin_sum <- bridge_noNA_bin %>% sf::st_set_geometry(value=NULL) %>% dplyr::group_by(都道府県) %>%
dplyr::summarise(dang_ratio=mean(dang)) %>% as.data.frame(stringsAsFactors=F)
#都道府県ポリゴンにdang_ratioを結合
bound_pref_dang_ratio <- bound_pref %>% dplyr::left_join(bridge_noNA_bin_sum,by=c("KEN"="都道府県"))
#dang_ratioを都道府県地図上にプロット
plot_pref <- ggplot2::ggplot() +
ggplot2::geom_sf(data=bound_pref_dang_ratio,aes(fill=dang_ratio),size=0.005)
plot_pref
# #[NOTRUN]プロットを画像に保存したければ実行
# ggplot2::ggsave(file="plot_pref_dang_ratio.png",plot=plot_pref,width=5,height=5,dpi=300)