cucumber flesh

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

tabularmaps: カラム地図で行政区の大まかな配置を可視化する

カラム地図と呼ばれるものがあります。 これは日本の47都道府県をはじめとした行政区の配置を表上に圧縮表示することで、それぞれの位置関係をわかりやすく伝えるためのプロジェクトです。

f:id:u_ribo:20200422175004p:plain

最近では、カラム地図の開発者の一人、福野泰介さん (@taisukef) による新型コロナウイルス対策ダッシュボードのページでの「現在患者数 / 対策病床数」の可視化で「カラム地図」を見る機会が増えた人もいるのではないでしょうか。

さて、この「カラム地図」をR言語の中で扱うべく、tabularmapsパッケージの開発を進めています。

github.com

以下、この記事ではtabularmapsの使い方を紹介していきます。

パッケージはCRANには登録されていないため、GitHub経由でインストールを行います。以下のコマンドを実行することでパッケージが利用可能になります。

install.packages("remotes")
remotes::install_github("uribo/tabularmaps")
library(tabularmaps)
library(ggplot2)

可視化の際にggplot2パッケージを利用しているので、そちらも読み込んでおく必要があります。

使い方

現在tabularmapsでは

  • 都道府県 (jpn77)
  • 東京23区 (tky23)
  • ISO-3166による国名 (iso3166)

の3種類に対応しています。括弧内の文字列がパッケージに含まれるデータフレームオブジェクトで、県や区などの対象行政名の位置関係を記録しています。

冒頭の都道府県のカラム地図は次のコマンドにより描画されます。

tabularmap(jpn77, 
           fill = region_kanji, 
           label = prefecture_kanji, 
           size = 3,
           family = "IPAexGothic") +
  theme_tabularmap(base_family = "IPAexGothic") +
  scale_fill_jpregion(lang = "jp",
                      name = "八地方区分")

tabularmap()がカラム地図を表示させる関数として機能します。この関数では第一引数に対象のデータフレームを与えて実行します。 そのほかの引数としてggplot2オブジェクトの審美的要素を決定するfill (塗り分けの色基準)、label (表示するラベル文字列)の変数名を与えます。 これに加え、ggplot2::geom_text()に渡す値を指定可能です。 これはラベルの大きさを調整する size やフォントのための family を使うことを想定しています。

上記のコードでは、tabularmap()に加えてtheme_tabularmap()scale_fill_jpregion()を実行しています。これらの関数はオプションとして利用可能なもので、 よりカラム地図っぽい見た目に調整するためのものです。 都道府県を表示させるカラム地図では地方別に塗り分ける色が決まっているため、それに応じた塗り分けを行うようになります。

続いて東京23区の表示をしてみます。先のコードからtabularmap()に与えるデータとラベル用の変数名を変更するだけです。

tabularmap(tky23,
           label = ward_kanji,
           family = "IPAexGothic") +
  theme_tabularmap(base_family = "IPAexGothic") +
  guides(fill = FALSE)

現実の配置や面積がぼやけてしまう欠点もありますが、 行政区域のデータを用意する必要がないのはカラム地図の大きな利点です。 (カラム地図のレイアウト自体はオープンライセンスで開発されています。)

また本来想定されていたダッシュボードでの利用を考える場合、 表示する際のスペースが有効活用できるのも嬉しいです。 日本の場合、細長くなるのでどうしても隙間ができてしまいますが、カラム地図だと基本的にグリッドになるので無駄が少ないです。

今後の機能

fukuno.jig.jp

本家のカラム地図では、市区町村を対象にしたレイアウトが全国47都道府県分定まったとのことです。 これらの情報を追いながら、tabularmapsパッケージでも市区町村版を追加する作業を進めていくつもりです。

全国のデータを追加するのは時間がかかりそうなので、皆さんからの協力(Pull request)をお待ちしています。

github.com

余談

SNSを見ていると、時々ヘンテコな地図が回ってきます。

白地図上に47都道府県の名前を書き込んでください」

中学生の地理で習うような問題です。 回答を見ていると、有名な東京都や神奈川県、形のわかりやすい北海道などは簡単な様子。 一方で北関東の群馬、栃木、茨城の配置、島根と鳥取の位置関係などを間違えている例をしばしば見かけます(私も苦手でした)。 またそれ以上にトンデモな解答も…。

行政境界の複雑な幾何学模様を認識し記憶することに対する人間の限界を感じます。 それが特に自分とは関わりのない地域であれば尚更です。

また、行政区は人間が勝手に定めたもの。時間の経過とともに境界も変化していきます。

情報を伝えるとき、これらの行政区をどこまで厳密に再現する必要があるでしょうか。

地図を表示するにはポリゴンないしラインデータを用意しなくてはいけない。そう思っていた時期が私にもありました。

ですがそうではなかったのです。

必要なのは、情報を正しく、わかりやすく伝えること。

カラム地図では、現実に対して幾らかの嘘をつくことになりますが、それは許容されうる範囲なのではと思います。

参考

気象庁提供の潮位表、過去の気象データをRで読み込む

最近、気象庁が提供する各種データをR上で取り扱うパッケージ、jmastatsに新しい関数を追加しました。 このパッケージはちまちまと更新しており、現在のところCRANへの登録は目指していません。 ですが便利なパッケージなのでこの場で宣伝しておきます。

https://gitlab.com/uribo/jmastats

gitlab.com

今回追加したのは以下の関数です。いずれも気象庁のウェブサイトでダウンロード可能なファイルを読み込むために使います。

  • read_tide_level(), pivot_tide_level(): 潮汐観測資料から潮位表データの読み込み
  • read_jma_weather(): 過去の地点気象データを読み込む

この機能を試したい場合、パッケージをGitLab経由(GitHubではありませんので注意)でインストールする必要があります。次のコマンドを実行するとインストールが行われます。

install.packages("remotes")
remotes::install_gitlab("uribo/jmastats")
library(jmastats)

次にこれらの関数の使い方を紹介します。

潮汐観測資料

潮汐観測資料および各種潮位データ・品質管理情報のページよりダウンロード可能な、 毎時潮位満潮・干潮の時刻と潮位が記録されたテキストファイルを与えて実行します。

ダウンロード済みのローカルファイルのパスを指定しても良いですし、気象庁のURLを直接与えても良いです。また引数による年月、地点の指定も可能です。次の例は、2020年2月の「東京」のデータを対象に実行するものです。

# ファイルが置かれているパスがわかる場合は直接パスを指定すると良いです。
#d <- 
#  read_tide_level("https://www.data.jma.go.jp/gmd/kaiyou/data/db/tide/suisan/txt/2020/TK.txt")
# path引数以外の引数でデータを指定する場合(こちらは一ヶ月単位のデータになります。)
d <- 
  read_tide_level(.year = 2020, .month = 2, .stn = "TK")

d
#> # A tibble: 29 x 42
#> hry_00   hry_01   hry_02   hry_03   hry_04   hry_05   hry_06   hry_07   hry_08   hry_09   hry_10
#> [cm]     [cm]     [cm]     [cm]     [cm]     [cm]     [cm]     [cm]     [cm]     [cm]     [cm]
#> 1      176      156      143      143      153      172      194      214      229      237      236
#> 2      195      180      166      160      162      172      188      204      219      230      232
#> 3      197      190      181      174      173      179      191      204      216      224      228
#> 4      204      206      204      201      197      194      194      196      203      211      218
#> 5      197      212      220      223      221      217      212      210      210      215      223
#> 6      172      196      213      222      225      221      211      199      188      182      183
#> 7      130      164      197      222      235      234      224      208      192      180      176
#> 8      112      148      187      224      251      262      258      241      217      196      182
#> 9       87      120      162      204      240      261      264      250      223      194      172
#> 10       69       93      134      183      229      263      279      276      256      226      195
#> # … with 19 more rows, and 31 more variables: hry_11 [cm], hry_12 [cm], hry_13 [cm], hry_14 [cm],
#> #   hry_15 [cm], hry_16 [cm], hry_17 [cm], hry_18 [cm], hry_19 [cm], hry_20 [cm], hry_21 [cm],
#> #   hry_22 [cm], hry_23 [cm], date <date>, stn <chr>, low_tide_hm_obs1 <time>, low_tide_level_obs1 [cm],
#> #   high_tide_hm_obs1 <time>, high_tide_level_obs1 [cm], low_tide_hm_obs2 <time>,
#> #   low_tide_level_obs2 [cm], high_tide_hm_obs2 <time>, high_tide_level_obs2 [cm],
#> #   low_tide_hm_obs3 <time>, low_tide_level_obs3 [cm], high_tide_hm_obs3 <time>,
#> #   high_tide_level_obs3 [cm], low_tide_hm_obs4 <time>, low_tide_level_obs4 [cm],
#> #   high_tide_hm_obs4 <time>, high_tide_level_obs4 [cm]

潮位の単位がcmで記録されているので、unitsパッケージを使って単位を与えているのが特徴です。 また、元データは毎時潮位と満潮・干潮の潮位データが一つのファイルになっているので列が多いです(行数は日数に対応)。そのため、jmastatsではこのデータを分離し、縦長のデータフレームに変換する関数 pivot_tide_level()を用意しました。

read_tide_level()で読み込んだデータを引数に与えて実行します。この関数の返り値は2つのデータフレームを含んだリストです。hourlyが毎時潮位、tideが満潮・干潮の潮位となります。

d_long <- 
  d %>% 
  pivot_tide_level()
d_long
#> $hourly
#> # A tibble: 696 x 3
#> datetime            stn   tide_value
#> <dttm>              <chr>       [cm]
#> 1 2020-02-01 00:00:00 TK           176
#> 2 2020-02-01 01:00:00 TK           156
#> 3 2020-02-01 02:00:00 TK           143
#> 4 2020-02-01 03:00:00 TK           143
#> 5 2020-02-01 04:00:00 TK           153
#> 6 2020-02-01 05:00:00 TK           172
#> 7 2020-02-01 06:00:00 TK           194
#> 8 2020-02-01 07:00:00 TK           214
#> 9 2020-02-01 08:00:00 TK           229
#> 10 2020-02-01 09:00:00 TK           237
#> # … with 686 more rows
#> 
#> $tide
#> # A tibble: 112 x 6
#> date       stn   tide_level count time   tide_value
#> <date>     <chr> <chr>      <chr> <time>       [cm]
#> 1 2020-02-01 TK    low        1     09:23         238
#> 2 2020-02-01 TK    high       1     02:33         141
#> 3 2020-02-01 TK    low        2     21:35         210
#> 4 2020-02-01 TK    high       2     15:23         155
#> 5 2020-02-02 TK    low        1     09:45         232
#> 6 2020-02-02 TK    high       1     03:16         160
#> 7 2020-02-02 TK    low        2     23:18         198
#> 8 2020-02-02 TK    high       2     16:23         151
#> 9 2020-02-03 TK    low        1     10:25         228
#> 10 2020-02-03 TK    high       1     03:41         173
#> # … with 102 more rows

せっかくなので可視化してみましょう。

library(ggplot2)
library(ggforce)
d_long %>% 
  purrr::pluck("hourly") %>% 
  ggplot(aes(datetime, tide_value)) + 
  geom_line(color = "red") +
  scale_x_datetime(date_labels = "%m/%d") +
  theme_light(base_family = "IPAexGothic") +
  labs(title = "毎時潮位グラフ 2020年2月",
       subtitle = "東京")

f:id:u_ribo:20200420231912p:plain

地点名がわからないときは tide_stationデータセットで見つけることができます。なお、年によって観測地点が変わるのでそこは注意です。

tide_station
#> Simple feature collection with 1670 features and 6 fields
#> geometry type:  POINT
#> dimension:      XY
#> bbox:           xmin: 122.95 ymin: 24.28333 xmax: 153.9833 ymax: 45.4
#> CRS:            EPSG:4326
#> # A tibble: 1,670 x 7
#> year  id    stn   station_name address                  type                  geometry
#> <chr> <chr> <chr> <chr>        <chr>                    <chr>              <POINT [°]>
#>   1 1997  1     WN    稚内         北海道 稚内市 新港町     フロート式     (141.6833 45.4)
#> 2 1997  2     AS    網走         北海道 網走市 港町       フロート式 (144.2833 44.01667)
#> 3 1997  3     HN    花咲         北海道 根室市 花咲港     フロート式 (145.5667 43.28333)
#> 4 1997  4     KR    釧路         北海道 釧路市 港町       フロート式 (144.3833 42.96667)
#> 5 1997  5     HK    函館         北海道 函館市 海岸町     フロート式 (140.7333 41.78333)
#> 6 1997  6     B3    小樽         北海道 小樽市 色内3丁目 音波式              (141 43.2)
#> 7 1997  7     SH    下北         青森県 むつ市 大字関根   音波式       (141.25 41.36666)
#> 8 1997  8     HC    八戸         青森県 八戸市 新湊3丁目 フロート式 (141.5333 40.53333)
#> 9 1997  9     MY    宮古         岩手県 宮古市 日立浜町   フロート式 (141.9833 39.63334)
#> 10 1997  10    OF    大船渡       岩手県 大船渡市 赤崎町   フロート式   (141.75 39.01667)
#> # … with 1,660 more rows

過去の地点気象データ

jmastatsには、すでにウェブスクレイピングにより気象データを取得する関数 jma_collect() がありました。一方でこの関数では地点数が多い・期間が長い場合に 気象庁のページへの負荷が多くなってしまう問題がありました。

今回追加した read_jma_weather() ではデータの用意こそユーザが各自で行う必要がありますが、任意の地点・期間のデータをRで読み込むには便利な関数となっています。 気象庁からダウンロードしたファイルのフォーマットは https://www.data.jma.go.jp/gmd/risk/obsdl/top/help3.html にあるようにやや煩雑となっていて、すぐに使える形式ではありません。 read_jma_weather()はRでの分析、可視化が素早くできるよう、自動的にデータを整形した形式で読み込みます。

関数を実行する前に対象のデータをダウンロードしておきましょう。 過去の気象データ・ダウンロードのページから関心の地点データを取得します。ここでは「つくば(舘野)」の「日平均気温」と「降水量の日合計」について、最近一ヶ月分のデータを用意しました。

d <- 
  read_jma_weather("~/Downloads/data.csv")
#> Selected station: つくば(館野)

#> Parsed with column specification:
#> cols(
#>   date = col_character(),
#>   `つくば(館野)_平均気温(℃)` = col_double(),
#>   `つくば(館野)_平均気温(℃)_品質情報` = col_double(),
#>   `つくば(館野)_平均気温(℃)_均質番号` = col_double(),
#>   `つくば(館野)_降水量の合計(mm)` = col_double(),
#>   `つくば(館野)_降水量の合計(mm)_現象なし情報` = col_double(),
#>   `つくば(館野)_降水量の合計(mm)_品質情報` = col_double(),
#>   `つくば(館野)_降水量の合計(mm)_均質番号` = col_double()
#> )

dplyr::glimpse(d)
#> Rows: 31
#> Columns: 8
#> $ date                                            <date> 2020-03-20, 2020-03…
#> $ `つくば(館野)_平均気温(℃)`                   <dbl> 12.9, 11.2, 14.7, 8.5, 6.0, 7.0,…
#> $ `つくば(館野)_平均気温(℃)_品質情報`          <dbl> 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, …
#> $ `つくば(館野)_平均気温(℃)_均質番号`          <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, …
#> $ `つくば(館野)_降水量の合計(mm)`               <dbl> 0.0, 0.0, 0.0, 0.0, 1.5, 0.0, 0.0…
#> $ `つくば(館野)_降水量の合計(mm)_現象なし情報`  <dbl> 0, 1, 1, 0, 0, 1, 1, 0, 0, 0, 1, 0, 0, …
#> $ `つくば(館野)_降水量の合計(mm)_品質情報`      <dbl> 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8, 8…
#> $ `つくば(館野)_降水量の合計(mm)_均質番号`      <dbl> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1…

お試しいただき、改善や機能追加のリクエストをいただけると助かります。 また、こうしたRパッケージの開発を支援してくださる方も引き続き募集しています。

https://uribo.hatenablog.com/entry/2020/03/27/170257 uribo.hatenablog.com

どうぞよろしくお願いします。

jpmeshバージョン1.2.0をリリース

統計調査などで使われる標準地域メッシュをRで扱うjpmeshパッケージのバージョン1.2.0 をCRANにリリースしました。以前のバージョンは1.1.3でした。マイナーアップデートですがいくつかの変更点・新機能がありますので紹介します。

f:id:u_ribo:20200328093427p:plain

インストール

2020年3月28日現在、WindowsおよびmacOSでのバイナリ版が用意されていない状況です。これらのOSを利用されている方はインストールをもうしばらくお待ちください。

install.packages("jpmesh") # v1.2.0
library(jpmesh)
library(sf)

メッシュサイズの指定を文字列から数値で行うように

これまでのjpmeshの関数では、メッシュサイズを指定する際に 80km のような文字列での指定が必要でした。

meshcode_set(mesh_size = "80km")

coords_to_mesh(141.3468, 43.06462, mesh_size = "500m")

rmesh(1, mesh_size = "1km")

これに対して、今回リリースされたv1.2.0では1kmを1として、すべてのメッシュサイズの指定を数値で行うように変更しました。すなわち、80kmのときは80、500mでは0.5、250mメッシュサイズの場合は0.25です。

上記のコードをv1.2.0で実行する際は以下のようにします。

meshcode_set(mesh_size = 80)

coords_to_mesh(141.3468, 43.06462, mesh_size = 0.5)

rmesh(1, mesh_size = 1)

メッシュサイズの変更 - mesh_convert()

次は追加した関数です。入力に与えたメッシュコードを適当なメッシュサイズに変更する mesh_convert() を用意しました。メッシュサイズのスケールアップ・スケールダウンが簡単に行えます。

set.seed(1)
rmesh(1) %>% 
  mesh_convert(to_mesh_size = 0.25)
##  [1] "4040745711" "4040745712" "4040745713" "4040745714" "4040745721"
##  [6] "4040745722" "4040745723" "4040745724" "4040745731" "4040745732"
## [11] "4040745733" "4040745734" "4040745741" "4040745742" "4040745743"
## [16] "4040745744"

スケールダウンの機能を持つ関数として fine_separate() がありますが、こちらは変換したいメッシュサイズにするまでに繰り返しの処理が必要になる場合がありました。

rmesh(1) %>% 
  fine_separate() %>% 
  purrr::map(fine_separate) %>% 
  purrr::reduce(c) %>% 
  unique()

mesh_convert()では 80km (80)から125m (0.125)までのメッシュサイズに対応しており、柔軟にメッシュサイズの変更ができます。

set.seed(2)
rmesh(1, mesh_size = 1) %>% 
  mesh_convert(to_mesh_size = 80)
## [1] "4931"

メッシュコードを変数に含むデータフレームのsf化 - meshcode_sf()

メッシュコードを含むデータフレームを処理することがあります。

set.seed(3)
df_mesh <- 
  tibble::tibble(
    meshcode = rmesh(1, mesh_size = 1) %>% 
      fine_separate(),
    value = rnorm(4))
df_mesh
## # A tibble: 4 x 2
##   meshcode   value
##   <chr>      <dbl>
## 1 492816651 -0.293
## 2 492816652  0.259
## 3 492816653 -1.15 
## 4 492816654  0.196

ここからメッシュコードのgeometry (polygon) を生成し、sfオブジェクトとして扱うための関数が meshcode_sf() です。引数 mesh_var にメッシュコードが記録されている変数名を与えて実行します。このとき、dplyrやtidyrで変数を指定するように引用符は必要ありません。

sf_mesh <- 
  df_mesh %>% 
  meshcode_sf(mesh_var = meshcode)
class(sf_mesh)
## [1] "sf"         "tbl_df"     "tbl"        "data.frame"
sf_mesh %>% 
  st_geometry() %>% 
  plot()

f:id:u_ribo:20200328092349p:plain

この機能は自分自身がよく使っていたものなので、実装できて満足しています。

次のリリースでは?

jpmeshは現在80kmメッシュから1kmメッシュまでの基準地域メッシュ、500m、250m、125mの分割地域メッシュに対応しています。国内で使われるメッシュコードには他に100mと50mの統合地域メッシュが存在します。これらについてサポートするかが議論されています。

github.com

100mメッシュが250mメッシュと同じ桁数になってしまうのでどうにかしないといけません。

また、対象のメッシュコードが人口集中地区(DID)をはじめとした国土政策関係の区域に含まれるかを判定する関数を実装予定でいます。

github.com

ご意見あるかたはGitHub もしくは Twitter (@u_ribo)まで。

スポンサーの検討もお願いいたします。

uribo.hatenablog.com