geofacetで日本のデータが利用可能になりました
もう数ヶ月も前の話ですが、CRANに登録されているgeofacetパッケージが更新されて、バージョンが0.1.9となりました。geofacetパッケージは、ggplot2ベースの作図パッケージの一種ですが、グラフを地理的な空間関係と対応付けて描画するという特徴があります。例えば、次のような図を見たことがある人がいるのではないでしょうか。
この図はアメリカでの2000年から2016年にかけての失業率を州に分けて表示しています。ただ図を並べるだけでなく、各州の空間関係を考慮した配置をすることで各州の関係がわかりやすく可視化されていると思います。この図はgeofacetのサンプルコードとして掲載されているもので、下記を実行することで描画されます。
library(ggplot2) library(geofacet) ggplot(state_unemp, aes(year, rate)) + geom_line() + facet_geo(~ state, grid = "us_state_grid2") + scale_x_continuous(labels = function(x) paste0("'", substr(x, 3, 4))) + ylab("Unemployment Rate (%)")
geofacetには、こうした図を書くために、いくつかの国や地域に対応していますが、バージョン0.1.9で日本にも対応するようになりました！追加した当人として、geofacetで日本のデータを利用する例を紹介します。アメリカに比べると少し不恰好で、現実の大きさや配置を再現できていない部分もあるかと思いますが、ぜひお使いください！
現在のgeofacetの日本データは、都道府県ごとに図を分割して描画するようになっています。なので都道府県ごとのデータを用意する必要があります。まずはグラフに描画するデータを用意しましょう。ここでは気象庁が提供する過去の気象データを用います。このデータをRから取得するパッケージを作成中ですが、その話はまたの機会にするとします。
全国47の各都道府県からそれぞれ1地点の観測所を選び、昨年6月の気象データを取得してきます。これらの操作の解説は省略しますが、jma_collect()という関数を利用して、気象庁の過去データを取得します。引数に、取得対象の項目（daily: 日毎の観測値）、対象の観測所、年月をそれぞれ指定して実行します。
library(dplyr) library(purrr) library(jpndistrict) library(jmastats) # 自作パッケージ data("station", package = "jmastats")
df_target_stations <- stations %>% # 観測装置の種類が「有線ロボット気象計」であるものに制限 # (いわゆるアメダス) filter(station_type %in% c("四")) %>% group_by(pref_code) %>% slice(1L) %>% ungroup() %>% select(pref_code, station_name, block_no)
df_weather <- df_target_stations %>% mutate(weather = pmap(., ~ jmastats::jma_collect("daily", block_no = ..3, year = 2017, month = 6) %>% dplyr::select(1, 2, 5))) %>% tidyr::unnest() %>% filter(!is.na(pref_code)) %>% sf::st_set_geometry(NULL) %>% parse_unit()
データを確認しておきます。気温と降水量についてのデータを持っていることがわかります。
glimpse(df_weather)
## Observations: 1,410
## Variables: 6
## $ pref_code <chr> "01", "01", "01", "01", "01", "01", "01",...
## $ station_name <chr> "宗谷岬", "宗谷岬", "宗谷岬", "宗谷岬", "宗谷岬", "宗谷岬",...
## $ block_no <chr> "1284", "1284", "1284", "1284", "1284", "...
## $ date <date> 2017-06-01, 2017-06-02, 2017-06-03, 2017...
## $ precipitation_sum_mm <S3: units> 7.5 mm, 7.0 mm, 0.0 mm, 0.0 mm, 0.0...
## $ temperature_average <S3: units> 8.1 °C, 5.6 °C, 4.8 °C, 4.6 °C, 7.2...
ではgeofacetを使い、全国の6月の気象データを可視化してみましょう。まずはベースとなるグラフを作成しますが、これはggplot2で普段作成するやり方と変わりません。
library(ggforce) # unitsオブジェクトをグラフに描画するのに使います p <- ggplot(df_weather %>% left_join(jp_prefs_grid1 %>% select(code_pref_jis, name), by = c("pref_code" = "code_pref_jis")), aes(date, temperature_average)) + geom_line(color = "steelblue")
ggplot2で、特定のグループ、ここでは都道府県別の観測所 (station_name)ごとに図を分けて表示したいときはfecet_wrap()を使います。
p + facet_wrap(~ station_name)
パネルに分割して、個々のデータが見やすくなりましたが、観測所の名前を見ても、見知らぬ土地の観測所では、北の寒い地域のデータなのか、南の暖かい地域のものなのかといったことがわかりにくいです。こうした時にgeofacetが役立ちます。
それでは先ほどの観測所ごとの表示を、空間配置を考慮して作成し直しましょう。geofacetの関数facet_geo()もggplot2::facet_wrap()同様に操作します。分割する変数を指定し、地図のデータセットをgrid引数に与えて実行するだけです。
この、地図のデータセットはgeofacetが持っているもので、ここに都道府県の名前や配置が定義されています。中身はただのデータフレームですが、一応確認しておきましょう。
glimpse(jp_prefs_grid1)
## Observations: 47
## Variables: 7
## $ code <int> 1, 2, 3, 5, 4, 6, 7, 8, 15, 9, 10, 11, 16, 17, 1...
## $ code_pref_jis <chr> "01", "02", "03", "05", "04", "06", "07", "08", ...
## $ name <chr> "Hokkaido", "Aomori", "Iwate", "Akita", "Miyagi"...
## $ name_abb <chr> "HKD", "AOM", "IWT", "AKT", "MYG", "YGT", "FKS",...
## $ name_region <chr> "Hokkaido", "Tohoku", "Tohoku", "Tohoku", "Tohok...
## $ col <int> 15, 14, 14, 13, 14, 13, 14, 13, 12, 12, 14, 13, ...
## $ row <int> 1, 2, 3, 3, 4, 4, 5, 5, 5, 6, 6, 6, 6, 6, 7, 7, ...
col, rowの変数が、geo_facet()を実行した際の配置となります。パネルの配置とデータを紐付ける変数が必要となりますが、今回は2桁の都道府県コード(code_pref_jis)を利用します。また、今回は観測所の名前をパネルのラベルとして利用したいので、jp_prefs_grid1の方にも都道府県コードに対応する観測所の名前を与えておく必要があります。
地図のデータセットの方に対応する変数がない場合は次のようにエラーとなります。
p + facet_geo(~ station_name, grid = "jp_prefs_grid1") # Error: The values of the specified facet_geo column 'station_name' do not match any column of the specified grid.
共通のラベルを用意するには以下のコードを実行します。
selected_station_names <- df_weather %>% distinct(pref_code, station_name) %>% arrange(pref_code) %>% magrittr::use_series(station_name) # 北海道(北)から5つの観測所の名前を出力 head(selected_station_names)
## [1] "宗谷岬" "小田野沢" "種市" "駒ノ湯" "休屋" "差首鍋"
# 変数を追加 jp_prefs_grid1 <- jp_prefs_grid1 %>% arrange(code_pref_jis) %>% mutate(station_name = selected_station_names)
p <- p + facet_geo(~ station_name, grid = "jp_prefs_grid1") p
できました！もう少し出力を整えていきましょう。
# remotes::install_github("yutannihilation/gghighlight") library(gghighlight) ggplot(df_weather %>% left_join(jp_prefs_grid1 %>% select(code_pref_jis, name), by = c("pref_code" = "code_pref_jis")), aes(date, temperature_average, color = pref_code)) + geom_line() + gghighlight(use_direct_label = FALSE) + facet_geo(~ station_name, grid = "jp_prefs_grid1") + scale_x_date(date_breaks = "1 week", date_labels = "%W") + guides(color = FALSE) + theme(axis.text = element_text(size = 5.2), strip.text.x = element_text(size = 7.6)) + labs(title = "June 2017", caption = "Data Source: Japan Meteorological Agency website\n (http://www.data.jma.go.jp/obd/stats/etrn/index.php)", x = "# Weeks", y = "Temperature")
また、流行りのgghighlightを使って次のようにしても綺麗だと思います (ラベルの文字サイズが大きすぎますね...修正する方法がわかりませんでした)。
追記: 早速 ラベルの文字サイズを調整する方法を教えていただきました(
label_params = list(size = ...)。更に、ラベルを非表示にする方法まで(use_direct_label = FALSE)！！ id: yutannihilation さん、ありがとうございます。
ラベルのサイズ変えるのはlabel_params = list(size = 0.5)という感じでいけるはず。ちゃんとドキュメント書いてなくてすいません。。
— Hiroaki Yutani (@yutannihilation) 2018年6月22日
# remotes::install_github("yutannihilation/gghighlight") library(gghighlight) ggplot(df_weather %>% left_join(jp_prefs_grid1 %>% select(code_pref_jis, name), by = c("pref_code" = "code_pref_jis")), aes(date, temperature_average, color = pref_code)) + geom_line() + gghighlight(use_direct_label = FALSE) + facet_geo(~ name, grid = "jp_prefs_grid1") + scale_x_date(date_breaks = "1 week", date_labels = "%W") + guides(color = FALSE) + theme(axis.text = element_text(size = 5.2), strip.text.x = element_text(size = 7.6)) + labs(title = "June 2017", caption = "Data Source: Japan Meteorological Agency website\n (http://www.data.jma.go.jp/obd/stats/etrn/index.php)", x = "# Weeks", y = "Temperature")
標準の配置が気に入らなかったり、自作の配置を適用したいときは「からだにいいもの」さんの
を参考にデータを編集すると良いです。
それでは。Enjoy!
このページで掲載している気象データおよび図は、気象庁ホームページ (http://www.data.jma.go.jp/obd/stats/etrn/index.php) のデータをもとに瓜生真也が加工・編集したものです。
