市区町村役所間の距離行列を求める~pdftoolsでのデータ抽出とsfによる算出編~
はじめに
Qiitaのなかで@3tkyさんが、都道府県庁間の距離行列をRで求める記事を書かれていました。
そして、挑戦状のような一文が残されています。
市区町村の役所の緯度経度は 国土地理院 がまとめているので、PDFを整形し、60進法表記を @uri さんのkuniezuパッケージ や pazerパッケージを活用して10進法に変換後、上のコードから市区町村役所間の距離行列が作成できそう。
これはやらねばなるまい（謎の使命感）。 本記事では、[@3tky](https://qiita.com/3tky) さんの書かれていた都道府県庁間の距離行列を算出する方法の別解と同様に市区町村役所の距離行列を求める方法を書きます。 大元の市区町村の役所の緯度経度データは国土地理院が公開するPDFデータです。そのためPDFからデータを抽出、加工する方法についても言及します。
記事の中で以下のパッケージを利用します。予め読み込んでおきましょう。
library(sf) library(ggplot2) library(rvest) # library(units) # library(tidyverse)
都道府県庁間の距離行列
記事を見たあと、CRANにリリース準備だった kuniezu パッケージに都道府県の県庁位置のデータを収納しました。 このデータの取得・整形するコードは
https://github.com/uribo/kuniezu/blob/master/data-raw/office_locaiton.R
に置いています。
data("jp47prefectural_offices", package = "kuniezu") # 47都道府県の県庁位置データ jp47prefectural_offices
## Simple feature collection with 47 features and 1 field
## geometry type: POINT
## dimension: XY
## bbox: xmin: 127.681114197 ymin: 26.2124996185 xmax: 141.346939087 ymax: 43.0641670227
## CRS: EPSG:4326
## # A tibble: 47 x 2
## office geometry
## <chr> <POINT [°]>
## 1 北海道庁 (141.346939087 43.0641670227)
## 2 青森県庁 (140.740005493 40.8244438171)
## 3 岩手県庁 (141.152496338 39.7036094666)
## 4 宮城県庁 (140.871948242 38.2688903809)
## 5 秋田県庁 (140.102493286 39.7186126709)
## 6 山形県庁 (140.363327026 38.2405548096)
## 7 福島県庁 (140.467773438 37.75)
## 8 茨城県庁 (140.446670532 36.3413887024)
## 9 栃木県庁 (139.883605957 36.5658340454)
## 10 群馬県庁 (139.060836792 36.3911094666)
## # … with 37 more rows
まずは対角要素が0となる距離行列を作成します。sfオブジェクトでの地物間の距離は st_distance() を使って求められます。
dist_pref <- jp47prefectural_offices %>% st_distance()
colnames(dist_pref) <- jp47prefectural_offices$office rownames(dist_pref) <- jp47prefectural_offices$office dist_pref[seq_len(6), seq_len(6)]
## Units: [m]
## 北海道庁 青森県庁 岩手県庁 宮城県庁
## 北海道庁 0.000000 253808.763222 373582.0205556 534013.529177
## 青森県庁 253808.763222 0.000000 129308.2699883 283959.259182
## 岩手県庁 373582.020556 129308.269988 0.0000000 161119.538069
## 宮城県庁 534013.529177 283959.259182 161119.5380689 0.000000
## 秋田県庁 385850.690256 134229.605213 90055.4114243 174198.478652
## 山形県庁 542058.359064 288699.213957 176229.8003421 44629.402174
## 秋田県庁 山形県庁
## 北海道庁 385850.6902557 542058.359064
## 青森県庁 134229.6052134 288699.213957
## 岩手県庁 90055.4114243 176229.800342
## 宮城県庁 174198.4786523 44629.402174
## 秋田県庁 0.0000000 165635.769370
## 山形県庁 165635.7693696 0.000000
また、距離の単位はmですが units パッケージの関数を使って任意の単位に変更可能です。
dist_pref <- dist_pref %>% units::set_units(km) dist_pref[seq_len(6), seq_len(6)]
## Units: [km]
## [,1] [,2] [,3] [,4] [,5]
## [1,] 0.000000000 253.808763222 373.5820205556 534.013529177 385.8506902557
## [2,] 253.808763222 0.000000000 129.3082699883 283.959259182 134.2296052134
## [3,] 373.582020556 129.308269988 0.0000000000 161.119538069 90.0554114243
## [4,] 534.013529177 283.959259182 161.1195380689 0.000000000 174.1984786523
## [5,] 385.850690256 134.229605213 90.0554114243 174.198478652 0.0000000000
## [6,] 542.058359064 288.699213957 176.2298003421 44.629402174 165.6357693696
## [,6]
## [1,] 542.058359064
## [2,] 288.699213957
## [3,] 176.229800342
## [4,] 44.629402174
## [5,] 165.635769370
## [6,] 0.000000000
この値が @3tky さんの結果や国土地理院の公表データと一致していることを確認します。
市区町村役所の距離行列
続いて市区町村役所の距離行列を求めます。 こちらのデータも国土地理院が位置情報を整理しています。 ですがPDFなので、Rで扱う際にはデータを抽出する必要が生じます。
PDFファイルのダウンロードは以下のコードで行います。
df_link <- read_html("https://www.gsi.go.jp/KOKUJYOHO/center.htm") %>% html_nodes(css = "div.base_txt > div:nth-child(5) > table > tbody > tr > td > a") %>% { tibble::tibble( name = stringr::str_remove(html_text(.), "\\[.+\\]"), link = html_attr(., "href")) } df_link$link %>% purrr::walk( ~ download.file(url = .x, destfile = basename(.x))
RでPDFのデータ抽出を行うパッケージはいくつかありますが、今回は
pdftools を使いました。
pdftools::pdf_text()によりテキストを抽出、若干血生臭い文字列処理を行いデータフレーム化します。
これらの処理を一括で実行する関数を書きました。
リンク先のコード（関数部分）をコピーペーストで読み込んでください。 最終的には各都道府県のPDFを引数に与える
gsi_office_extract() を実行することでデータ抽出が完了します。
# 茨城県の市区町村役所の位置情報 sf_pref08office <- gsi_office_extract("ibaraki_heso.pdf")
sf_pref08office
## Simple feature collection with 45 features and 1 field
## geometry type: POINT
## dimension: XY
## bbox: xmin: 139.745285034 ymin: 35.8577766418 xmax: 140.751113892 ymax: 36.8019447327
## CRS: EPSG:4326
## # A tibble: 45 x 2
## office geometry
## <chr> <POINT [°]>
## 1 茨城県庁 (140.446670532 36.3413887024)
## 2 水戸市役所 (140.471664429 36.3658332825)
## 3 日立市役所 (140.650558472 36.5988883972)
## 4 土浦市役所 (140.204162598 36.0783348083)
## 5 古河市役所 (139.755004883 36.1783332825)
## 6 石岡市役所 (140.286941528 36.1905555725)
## 7 結城市役所 (139.876663208 36.305557251)
## 8 龍ケ崎市役所 (140.182220459 35.9116668701)
## 9 下妻市役所 (139.967498779 36.1844444275)
## 10 常総市役所 (139.993896484 36.0236129761)
## # … with 35 more rows
この位置情報データに対して、都道府県庁間の距離行列を求めた時と同じコードを実行することで市区町村役所の距離行列も求められます。 これらの処理も関数化しておきました。可視化と合わせてどうぞ。
# 距離行列を作成 st_distmatrix <- function(data, var) { res <- data %>% st_distance() %>% units::set_units(km) vars <- data %>% purrr::pluck(var) colnames(res) <- vars rownames(res) <- vars res } # ggplot2で描画するためのデータ整形。三角行列に変換して縦長のデータにします matrix_to_longer <- function(data) { res_mt <- as.matrix(units::drop_units(data)) res_mt[lower.tri(res_mt)] <- NA res_mt %>% as.data.frame() %>% tibble::rownames_to_column(var = "from") %>% tibble::as_tibble() %>% tidyr::pivot_longer(cols = seq.int(2, ncol(.)), names_to = "to", values_to = "dist") %>% dplyr::filter(!is.na(dist)) %>% dplyr::mutate_at(dplyr::vars(from), list(~ forcats::fct_inorder(.))) %>% dplyr::mutate(to = forcats::fct_rev(forcats::fct_inorder(to))) } # ggplot2でのプロット plot_distmatrix <- function(data) { data %>% ggplot(aes(from, to)) + geom_tile(aes(fill = dist), color = "white") + scale_fill_viridis_c() + guides(fill = guide_colorbar(title = "distance (km)")) + theme_bw(base_family = "IPAexGothic") + theme( axis.text.x = element_text(angle = 90), axis.ticks = element_blank(), axis.line = element_blank(), panel.border = element_blank(), panel.grid.major = element_line(color = '#eeeeee') ) }
st_distmatrix(sf_pref08office, "office") %>% matrix_to_longer() %>% plot_distmatrix()
Enjoy!
