cucumber flesh

Rを中心としたデータ分析・統計解析らへんの話題をしていくだけ

Tokyo.R#87でdata.tableパッケージの発表をしました

先日開催された第87回Tokyo.Rにおいてdata.tableパッケージに関する発表を行いました。 ここでは発表資料へのリンクと補足を書きます。

今回のスライド作成に関して、GitHubスポンサーの皆さんからの援助を受けています。 どうもありがとうございました。この場を借りて感謝いたします。

発表内容

当日の資料をspeakerdeckにアップロードしています。

speakerdeck.com

今こそ、data.tableを学ぼう! / datatable1130 - Speaker Deck

共有した画面上でカラーピッカーがデカデカと表示されたまま進行してしまったようで、どうもすみませんでした。 ご指摘ありがとうございました。

dplyrとの使い分け

data.tableパッケージはしばしば、データ操作という観点でdplyrパッケージとともに引き合いに出されます。 Rでのデータ操作といえばdplyr!そういう流れを感じますし、自分自身もdplyrファンで普段のコードもdplyr含めたtidyverseのものを使っています。 一方で時々、人様のコードを読む機会があるとdata.tableパッケージを使うものに遭遇します。 その時、data.tableで処理されたコードが読めないと困ります。

data.tableとdplyrを使ったデータ操作のコードの大きな違いは、

  • data.table: データを宣言し、[の内部でデータへの処理を記述する
  • dplyr: verbと呼ばれるデータに対する振る舞いを関数化、必要に応じてパイプ演算子%>%)によって処理を繋げていく

です。dplyrに馴染んでいる私は、パイプ演算子の出てこないコードに立ちすくんでしまうことがありました。 (逆にdata.tableをメインに使うのであったらパイプ演算子が出てくるコードを見ても謎に感じるのではないでしょうか)

data.tableでの処理は(覚えてしまうと)シンプルです。dplyrのように複数の関数の挙動を覚える必要がありません。 対象のデータ(data.tableオブジェクト)に対して [演算子を使い、その中で次の画像に示す規則で処理を記述することになります。

f:id:u_ribo:20200806085002j:plain

data.tableの文法を効率的に学習するにはdplyrでの処理と比較したら良いのではと思い、今回のスライドではdplyrのコードも掲載しています。

現実的な使い分けについてはTwitterで頂いたコメントへの返信の通り

普段使いのコードは馴染みのあるものを、ただ必要に迫られた際のためにどちらも覚えておくと良い、 と身も蓋もない話になってしまいますがそう思います。

スライド作成

フォント

発表後のチャット欄で「フォントは何を使っているか」という質問がありました。今回のスライドで用いたフォントは以下の通りです。

  • Canela Text
  • Font Awesome 5 Brands
  • FOT-ラグランパンチ
  • FOT-ロックンロール
  • Hack
  • Rockwell

このうち、「ラグランパンチ」、「ロックンロール」はフォントワークスmojimoサービスが提供するフォントパックに含まれるものです。「ラグランパンチ」はアニメ「キルラキル」やゲームソフトの「ゼルダの伝説 Breath of the world」で使われているものです。インパクトがあって見出しに適していますね。お気に入りです。

f:id:u_ribo:20200806090300j:plain

36書体が、年間定額

利用するコンピュータは一台のみの制約がありますが、切り替えもシンプルなので年間3,600円はお得です。 フォントを契約する余裕が生まれたのはGitHubスポンサーで支援してくださる方がいたためです。 発表の機会が増えるなら買ってしまえ、と契約に至りました。

カラーパレット

pals::pal.bands(c("#040A0F", "#EB6257", "#E98353", "#F1B952", "#F3CF50", "#F4EDCD"))

f:id:u_ribo:20200806083349p:plain

  • #040A0F
  • #EB6257
  • #E98353
  • #F1B952
  • #F3CF50
  • #F4EDCD

最初は違う色を使っていましたが、どうも馴染まなくて一部を変更しました。そのせいもあってか、普段は4~5色に抑えようとしているので、少し多くなってしまいました。 ドイツっぽいなーと思いますw

Enjoy!

市区町村役所間の距離行列を求める~pdftoolsでのデータ抽出とsfによる算出編~

はじめに

Qiitaのなかで@3tkyさんが、都道府県庁間の距離行列をRで求める記事を書かれていました。

qiita.com

そして、挑戦状のような一文が残されています。

市区町村の役所の緯度経度は 国土地理院 がまとめているので、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()によりテキストを抽出、若干血生臭い文字列処理を行いデータフレーム化します。 これらの処理を一括で実行する関数を書きました。

https://github.com/uribo/kuniezu/blob/22a17ac8ff3567095ca2a934bac9e3d6ccbc0820/data-raw/office_locaiton.R#L5-L59

リンク先のコード(関数部分)をコピーペーストで読み込んでください。 最終的には各都道府県の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()

f:id:u_ribo:20200509170948p:plain

Enjoy!

kuniezu: 日本の国土地理を扱いやすくするRパッケージをCRANに登録しました

https://repository-images.githubusercontent.com/259226015/c936a500-8e9e-11ea-81e1-8c7635268123

はじめに

kuniezuパッケージ (v0.1.0) をCRANにリリースしました。

github.com

このパッケージは、私が業務や趣味で日本国内の地理空間データを扱う時に作っていた関数を一つのパッケージに整理したものです。 空間的に世界規模のデータを扱うのではなく、日本国内に限った話であれば、日本に即した仕様や座標参照系を利用した方が良いことがあります。 そうした日本の地理空間データを処理する際に利用することがある機能や、あると便利なデータセットを提供できるように努めています。

ゆるゆると開発するつもりでいましたが、Twitterで開発中であることをつぶやいたところ想定以上に反響が大きかったのでCRAN登録を急ぐことにしました。

この記事では、そんなkuniezuパッケージの機能紹介を行います。

まずはみなさんパッケージがインストールされていないと思うので、利用されたい方は以下のコマンドを実行してインストールを済ませてください。

install.packages("kuniezu")
library(kuniezu)
# 使い方を説明するために以下のパッケージも読み込んでおきます
library(sf)
library(ggplot2)
library(leaflet)

使い方

最初のリリースでは以下の機能を実装しました。

  • DMS表記を十進数表記に変換
  • 日本測地系2011における平面直角座標系の特定
  • 南西諸島・小笠原諸島を移動した日本地図の描画
  • 地理院タイルをleafletで簡単に利用できるように
  • 国土地理に関するいくつかのデータセット

順に解説します。

parse_*_dohunbyo(): DMS表記を十進数表記に変換

地球上の任意の座標を示す際に使われる緯度と経度の表記には、十進数で表される場合と「度(Degree)・分(Minute)・秒 (Second)」を用いて表記される場合があります。 後者はDMS表記と呼ばれ、十進数で139.741435.6581と表される座標に対してE139°44′28.8869″N35°39′29.1572″のように表現します。 更に日本では度分秒の記号に対して漢字使って東経139度44分28.8869秒北緯35度39分29.1572秒等の表記が用いられることがしばしばあります。

この日本独自のDMS表記を十進数に変換する関数がparse_*_dohunbyo()です(アスタリスクの部分にはlonlatがそれぞれ指定できます)。 lonは東経、西経を示す(180~-180)経度、latが北緯、南緯(90~-90)の緯度を扱うために使い分けます。

# 日本経緯度原点(東京都港区麻布台二丁目)の座標
parse_lon_dohunbyo("東経139度44分28.8869秒")
## [1] 139.7414
parse_lat_dohunbyo("北緯35度39分29.1572秒")
## [1] 35.6581

これは過去のブログ記事に書いた関数をパッケージに含めたものです。 加えて、十進数に変換する前段階として漢字を記号に処理するreplace_dohunbyo_kanji()も用意しています。

replace_dohunbyo_kanji("東経139度44分28.8869秒")
## [1] "E139°44’28.8869."

日本測地系2011における平面直角座標系の特定

日本の現行の測地系である日本測地系2011には、19の平面直角座標系の区分がなされています。

平面直角座標系は狭い範囲の測量や距離の算出に用いるのに適しています。 任意の経緯度を平面直角座標に変換するには、まず座標が平面直角座標系のどの系番号に含まれるのかを知っておく必要が生じます。 これを求めるのがst_detect_jgd2011()です。

sfパッケージでおなじみのst_*()が関数名につくように、関数にはsfcオブジェクトを与えて実行します。

# 日本経緯度原点の座標
st_detect_jgd2011(st_sfc(sf::st_point(c(139.7414, 35.6581)), crs = 4326))
## [1] 6677

move_jpn_rs(): 南西諸島の一部・小笠原諸島を移動した日本地図の描画

日本だけを地図上に描画する際、表示する画面の都合からか、南西諸島の一部や小笠原諸島の位置を海上に移動することがあります(例えばNHKの気象情報のページ白地図)。 これを実現するには、予め移動対象の地物に対する操作が必要となります。

move_jpn_rs()はこれらの島々の位置を変更するための関数です。 地図描画のために奄美群島沖縄県小笠原諸島の位置を変更する処理が実行されます。 また、この地図化もパッケージの関数で行えるようにggplot2ベースのgeom_jpsegment()を用意しました。

使い方は次のように、まず都道府県や市区町村レベルでポリゴンが分かれた地理空間データに対してmove_jpn_rs()を適用、その後マッピングを行います。

move_jpn_rs(jgd2011_bbox) %>%
  ggplot() +
  geom_sf() +
  geom_jpsegment()

デフォルトの設定で移動した地物がわかるように境界線が引かれます。

この関数の元も、かつて記事に書いたものです。

少し前に @shoei05 さんがブラッシュアップしてくれたのをきっかけにパッケージへ含めました。

一方で、まだこの関数はさらなる磨き上げが必要で、解決策がわかる人はぜひ教えてください。

github.com

地理院タイルをleafletで簡単に利用できるように

インタラクティブなウェブ地図を利用したことがある方ならわかるかと思いますが、leafletでは任意の地図タイルを背景画像として扱うことができます。 国土地理院が発行する地理院タイルもまた、この背景画像として利用可能で、 一覧ページに示されるように多種多様な地図タイルが整備されています。

kuniezuパッケージではこの地理院タイルをleafletで簡単に導入できるようになっています。 これらはgsi_tilesオブジェクトを通して呼び出します。

現在のバージョンでは「近年の災害」をはじめとして、すべての地理院タイルに対応しているわけではありません。 kuniezuで利用可能な地理院タイルの一覧は以下のコードで確認できます。

names(gsi_tiles)
##  [1] "standard"             "pale"                 "english"             
##  [4] "lcm25k_2012"          "lcm25k"               "ccm1"                
##  [7] "ccm2"                 "vbm"                  "vbmd_bm"             
## [10] "vbmd_colorrel"        "vbmd_pm"              "vlcd"                
## [13] "lum200k"              "lake1"                "lakedata"            
## [16] "blank"                "seamlessphoto"        "ortho"               
## [19] "airphoto"             "gazo4"                "gazo3"               
## [22] "gazo2"                "gazo1"                "ort_old10"           
## [25] "ort_USA10"            "ort_riku10"           "relief"              
## [28] "anaglyphmap_color"    "anaglyphmap_gray"     "hillshademap"        
## [31] "earthhillshade"       "slopemap"             "slopezone1map"       
## [34] "afm"                  "lcmfc2"               "lcmfc1"              
## [37] "swale"                "cp"                   "jikizu2015_chijiki_d"
## [40] "jikizu2015_chijiki_i" "jikizu2015_chijiki_f" "jikizu2015_chijiki_h"
## [43] "jikizu2015_chijiki_z" "jikizu_chijikid"      "jikizu_chijikii"     
## [46] "jikizu_chijikif"      "jikizu_chijikih"      "jikizu_chijikiz"

使い方は次のようになります。leafletおよびmapviewでの例を紹介します。

# 出力結果は省略します
gsi_tiles$standard %>%
  addCircles(
    data = sf::st_transform(extreme_points %>%
      purrr::reduce(c),
      crs = 4326))

mapview::mapview(st_sfc(sf::st_point(c(139.7414, 35.6581)), crs = 4326), 
                 map = gsi_tiles$pale)

なお、これらの地理院タイル利用する際はタイルごとに設定されている利用規約に従ってください。

国土地理に関するデータセット

これまでの例で出てきたものがありますが、日本の国土地理を扱う上であると便利なデータ等をパッケージで提供するようにしました。 現在、次のようなデータがあります。

  • extreme_points… 東西南北端点の経度緯度データ
  • jgd2011_bbox日本測地系2011の平面直角座標系の区域を示すポリゴンデータ
  • jp47prefectural_offices… 47都道府県の庁舎の位置データ

こんなデータや機能が欲しいという方はぜひ issue に書き込んでください。 バグ報告もありがたいです。 GitHub以外ではTwitterでコメントいただいても対応します。

Enjoy!