jmastatsパッケージによる気象庁の気象データの取得〜過去100年間の気温のヒートマップを作成するまで〜

はじめに

9月に入っても気温の高い日が続いています。「近年の夏の暑さは過去と比べてどうなのか」、気象庁の過去の気象データを使って分析・可視化する例が、 @kaz_ogiwara さんによるSNSの投稿（元は東洋経済ONLINEに掲載された記事）をきっかけに、各新聞社の記事としても取り上げられています。

【朝日新聞】4月の暑い日、やっぱり増えてる？　ヒートマップで確認 (東京)
【朝日新聞】今年8月も過去126年で最も暑かった　2010年の記録を大幅更新 (千葉など)
【秋田魁新報】 140年で一番暑い夏　秋田市の8月平均気温30度、観測史上最高 (秋田)
【西日本新聞】暑くなった福岡、一目瞭然　134年で4度上昇…8月の気温を可視化【14地点ヒートマップ】 (福岡)

このような分析を行うのには、気象庁が公開するデータの取得が必要です。データは誰でも利用可能ですが、データ分析のためにウェブサイトから収集するのは面倒です。データダウンロード機能も用意されていますが、データ容量の制限があったり、手作業が伴う煩わしさからは解放されません。そこで、気象庁の各種データをRから取得・操作するのに役立つRパッケージを作成しました。パッケージ名はjmastatsです。

この記事では、過去の気象データをはじめ、jmastatsパッケージの使い方を紹介するとともに、「徳島」を例として、1923年から2022年の過去100年間の8月の日ごとの観測値をヒートマップで可視化してみます。

パッケージの読み込み

jmastatsパッケージはCRANからインストールできます。次のコマンドを実行してインストールが行われます。

install.packages("jmastats") # v.0.2.0

気温ヒートマップの作成のため、主役となるjmastatsパッケージのほか、いくつかのパッケージを読み込みます。

library(jmastats)
library(tidyverse) # dplyr, lubridateなどのデータ操作用
library(units) # 
library(patchwork) # 複数の図を結合してよしなに配置するためのパッケージ
theme_set(theme_bw()) # ggplot2のテーマ指定

気象データの取得

過去の気象データを取得するために、まずは対象の観測所を選択しましょう。今回は徳島県の徳島市に位置する「徳島」観測所を対象とします。

パッケージに付属のデータセット stations には全国の気象観測地点（AMeDAS）の情報が含まれます。観測地点の配置された標高や観測種目について記録されていますが、この中に気象データの取得に必要な”block_no”が存在します。

# 「徳島」のblock_noを調べます
stations |>
  filter(station_name == "徳島") |>
  distinct(block_no)

# A tibble: 1 × 1
  block_no
  <chr>   
1 47895

block_noがわかれば、jma_collect()関数を使って気象データの取得が行えます。他に必要なのは観測データの種類と期間についてです。今回は各年の8月の日ごとの平均気温を取得するため、jma_collect()関数の引数に次のように指定します。

jma_collect(
  # 日ごとの観測値を指定
  item = "daily",
  # 先ほど確認した「徳島」のblock_no
  block_no = "47895",
  # 参照したい年と月をそれぞれ指定
  year = 2022,
  month = 8)

今回は1922年から2022年にかけての値が対象となります。ここで重要なのは、年以外の引数は変更する必要がないことです。つまり、取得対象は日ごとの観測値、地点は「徳島」、また、期間は8月中であることは変わりません。このような場合、年の引数を変えたコードを繰り返し記述するのも良いですが、繰り返しの処理はプログラムが得意とするところです。以下のようにすると、year以外の引数を固定し、yearに任意の値を指定した結果を一度に得ることが可能となります。

df_tksm <-
  # year引数に与えられる年を指定 ... 1923年から2022年まで
  seq.int(1923, 2022) |>
  purrr::set_names() |>
  purrr::map(
    \(x) jma_collect(item = "daily",
                     block_no = "47895",
                     year = x,
                     month = 7)) |>
  # 一つのデータフレームにまとめる
  # 取得を行った年をyear変数に記録する
  purrr::list_rbind(names_to = "year")

気象庁のサーバーへの負荷を軽減するため、jma_collect()関数では取得したデータをローカルストレージにキャッシュとして保存します。そのため、一度取得したデータについては、再度取得する必要はありません。実行間隔についても同様に、気象庁サーバーへ過剰な負荷を与えないように、一定時間の間隔を空けてデータを取得するようになっています。そのため、対象の年数が多い（かつキャッシュが存在しない）場合は、取得に時間を要します。今回も対象の年数が多いのでコードを実行したらしばらく休憩しましょう。

キャッシュ機能は便利ですが、データの更新が行われた場合は、キャッシュを削除する必要があります。また、キャッシュを利用せずにデータを取得する引数オプションcache = FALSEも利用できますので覚えておくと良いでしょう。

さて、取得したデータを確認してみましょう。

glimpse(df_tksm)

Rows: 3,100
Columns: 10
$ year          <chr> "1923", "1923", "1923", "1923", "1923", "1923", "1923", …
$ date          <date> 1923-07-01, 1923-07-02, 1923-07-03, 1923-07-04, 1923-07…
$ pressure      <tibble[,2]> <tbl_df[26 x 2]>
$ precipitation <tibble[,3]> <tbl_df[26 x 3]>
$ temperature   <tibble[,3]> <tbl_df[26 x 3]>
$ humidity      <tibble[,2]> <tbl_df[26 x 2]>
$ wind          <tibble[,5]> <tbl_df[26 x 5]>
$ sunshine      <tibble[,1]> <tbl_df[26 x 1]>
$ snow          <tibble[,2]> <tbl_df[26 x 2]>
$ weather_time  <tibble[,2]> <tbl_df[26 x 2]>

pressure 気圧
precipitation 降水量
temperature 気温
humidity 湿度

のように観測種目ごとに変数が分かれています。jma_collect()関数の返り値はtibbleデータフレームですが、格納される変数は別のデータフレームを内包しています。つまり以下のような階層構造のあるデータフレームとなっています。

気象庁の過去の気象データでは、一つの観測種目ごとに複数の項目が含まれるという特徴があります。例えば気温については、平均気温、最高気温、最低気温の3項目が含まれます。どの気象データを含んでいるかは対象の観測所によりますが、こうした複数の項目を通常の横並びの変数として扱うと列数が増えて処理が面倒になることがあります。

そこでjmastatsでは、観測種目ごとにデータフレームとしてまとめることで、関心のある観測種目のみを容易に参照できるようになっています。今回は気温についてのみを対象とするので他の変数は必要ありません。そのためデータフレームの列選択を行うdplyr::select()関数を使って簡単に気温のデータだけを取り出すことができます。

df_tksm <-
  df_tksm |>
  # 日付と気温についての変数を選択
  select(year,
         date,
         temperature) |>
　 # ネストされたデータフレームを展開する
  tidyr::unnest(cols = temperature)

glimpse(df_tksm)

Rows: 3,100
Columns: 5
$ year         <chr> "1923", "1923", "1923", "1923", "1923", "1923", "1923", "…
$ date         <date> 1923-07-01, 1923-07-02, 1923-07-03, 1923-07-04, 1923-07-…
$ `average(℃)` <dbl> 22.3, 22.4, 20.7, 22.6, 20.5, 22.9, 22.9, 21.7, 25.1, 23.…
$ `max(℃)`     <dbl> 26.0, 25.3, 24.3, 27.8, 23.2, 28.5, 28.5, 24.2, 31.0, 27.…
$ `min(℃)`     <dbl> 20.3, 20.3, 18.0, 16.4, 18.6, 18.3, 17.9, 19.9, 20.4, 21.…

tidyr::unnest()関数で気温について記録したtemperature変数を展開すると、平均気温、最高気温と最低気温の値が記録されていることがわかります。

つづいてjmastatsのもう一つの便利な関数を紹介します。気象データには、記録した値の単位が変数名として与えられています。変数に単位が記録されていると、すぐに単位を確認することができますが、データの扱いには不便です。例えば、気温の場合、単位は℃ですが、変数名に単位が含まれていると、変数名を参照するのに単位も入力することになるので手間がかかります。そこでjmastatsでは、変数名から単位を取り除くparse_unit()関数を用意しています。この関数は単純に変数名から単位を除くだけでなく、単位に応じて値をunitsオブジェクトに変換します。unitsオブジェクトとはunitsパッケージにより提供される、測定単位を扱うためのクラスです。unitsオブジェクトの詳細はunitsパッケージのドキュメントを参照してください。

df_tksm <-
  df_tksm |>
  #　変数名から単位を取り除き、変数の値を対応するunitsオブジェクトに変換
  parse_unit()

glimpse(df_tksm)

Rows: 3,100
Columns: 5
$ year    <chr> "1923", "1923", "1923", "1923", "1923", "1923", "1923", "1923"…
$ date    <date> 1923-07-01, 1923-07-02, 1923-07-03, 1923-07-04, 1923-07-05, 1…
$ average [°C] 22.3 [°C], 22.4 [°C], 20.7 [°C], 22.6 [°C], 20.5 [°C], 22.9 [°C…
$ max     [°C] 26.0 [°C], 25.3 [°C], 24.3 [°C], 27.8 [°C], 23.2 [°C], 28.5 [°C…
$ min     [°C] 20.3 [°C], 20.3 [°C], 18.0 [°C], 16.4 [°C], 18.6 [°C], 18.3 [°C…

気温についての3つの項目の単位はいずれも摂氏（℃）でしたが、これらは変数名からはなくなり、各値についていることがわかります。

最後に、可視化に向けてちょっとした処理を追加しておきます。

df_tksm <- 
  df_tksm |>
  mutate(year = as.integer(year),
         day = lubridate::day(date))

気温のヒートマップ

jmastatsパッケージにはヒートマップ作成の機能はありません。ですが取得したデータをもとに、ggplot2パッケージを使って簡単にヒートマップを作成できます。今回はヒートマップを描画するために次のような関数を用意しました。

plot_temp_hm <- function(data, var, ...) {
  pal_relative <-
    jmastats:::jma_pal(palette = "relative", .attribute = TRUE)
  data |>
    mutate(year = forcats::fct_rev(as.character(year)),
           {{ var }} := drop_units({{ var }})) |>
    mutate(z = case_when(
      {{ var }} >= 35 ~ "35~",
      between({{ var }}, 30, 35) ~ "30~35",
      between({{ var }}, 25, 30) ~ "25~30",
      between({{ var }}, 20, 25) ~ "20~25",
      between({{ var }}, 15, 20) ~ "15~20",
      .default = as.character({{ var }})
    )) |>
    ggplot() +
    aes(day, y = year, fill = z) +
    geom_tile(color = "white", linewidth = 0.05, ...) +
    coord_equal() +
    scale_fill_manual(values = rev(pal_relative$colors |>
                                     purrr::set_names(pal_relative$amedas$labels)),
                      drop = TRUE) +
    scale_y_discrete(breaks = as.character(seq.int(1920, 2020, by = 10))) +
    guides(fill = guide_legend(title = NULL,
                             reverse = TRUE))
}

p1_ave <- 
  plot_temp_hm(df_tksm, average, show.legend = TRUE)

p1_ave +
  labs(title = "「徳島」（徳島県徳島市）の8月の気温",
       subtitle = str_glue("平均気温({intToUtf8(c(176, 67))})"),
       caption = "1923年から2022年の8月の各日の観測値\nデータ源: 気象庁")

出力されたヒートマップを見ると、2000年代はほぼ毎年、平均気温が35℃を超える日があることがわかります。また直近では2018年に平均気温35℃を超える日が連続していたことも確認できました。

平均気温のほかに最高気温、最低気温についても同様の図を作成します。作図のためのコードを関数化しているので、同じような図を作成する場合は、対象の変数を変えるだけで簡単に作図が行えます。

p1_ave <- 
  plot_temp_hm(df_tksm, average, show.legend = FALSE) +
  labs(subtitle = str_glue("平均気温({intToUtf8(c(176, 67))})"))
p2_max <-
  plot_temp_hm(df_tksm, max, show.legend = FALSE) +
  labs(subtitle = str_glue("最高気温({intToUtf8(c(176, 67))})"))
p3_min <-
  plot_temp_hm(df_tksm, min, show.legend = FALSE) +
  labs(subtitle = str_glue("最低気温({intToUtf8(c(176, 67))})"))

凡例も用意して。

p_legend <-
  df_tksm |>
  tidyr::pivot_longer(cols = c(average, min, max)) |>
  mutate(value = drop_units(value)) |>
  mutate(z = case_when(
    value >= 35 ~ "35~",
    between(value, 30, 35) ~ "30~35",
    between(value, 25, 30) ~ "25~30",
    between(value, 20, 25) ~ "20~25",
    between(value, 15, 20) ~ "15~20",
    .default = as.character(value)
  )) |>
  filter(!is.na(z)) |>
  ggplot() +
  aes(day, value, fill = z) +
  geom_tile() +
  scale_fill_manual(values = rev(pal_relative$colors |>
                                   purrr::set_names(pal_relative$amedas$labels)),
                    drop = TRUE) +
  guides(fill = guide_legend(title = NULL,
                             reverse = TRUE))

p_legend <- cowplot::get_legend(p_legend)

複数の図をpatchworkパッケージで結合しましょう。

p1_ave + p2_max + p3_min +
  p_legend +
  plot_layout(widths = c(0.3, 0.3, 0.3, 0.1), ncol = 4, guides = "auto") +
  plot_annotation(title = "「徳島」（徳島県徳島市）の8月の気温",
                  caption = "1923年から2022年の8月の各日の観測値\nデータ源: 気象庁")

最高気温が35℃を超える日が目立ちます。1950年代より前にもありますが、これは何でしょうか？最低気温については、平均気温や最高気温よりも近年の気温の上昇の様子がうかがえます。昔は8月中でも最低気温が15〜20℃だった日が数日はあったようですが、2000年代に入ると最低気温も上昇し、20℃を下回る日はほとんどありません。8月中頃からは最低気温が連日25〜30℃となる日が続く傾向が見られます。

おまけ: 気象情報の配色

さて、この記事のなかで使われた配色にピンと来た人はするどい（？）人です。ヒートマップの配色は気象庁が定める「気象情報の配色に関する設定指針(PDF)」に従って行いました。この指針は気象情報に関する注意・警戒の喚起効果を高め、高齢者にも配慮した配色を行うことを目的としています。

例えば、気象庁の指針では、「中間レベルの値は標準的で、値が大きくまたは小さくなるほど注意を促したい情報の配色」としてアメダスおよび天気予報分布の気温を示す際に、次の配色が利用されます。

pal_relative$colors |> 
  set_names(pal_relative$amedas$labels) |> 
  pals::pal.bands()

国内の気象データを可視化する際は、気象庁の指針に従う配色を行うと統一感が出て良いと個人的には思います。ですのでjmastatsでは気象情報の配色を再現するためのオブジェクトを用意しています（エクスポートされていないので、利用の際にはjmastats:::jma_pal()のように:::を使って関数を呼び出す必要があります）。

おわりに

気象庁の過去の気象データをRから取得するためのパッケージjmastatsを紹介しました。気象データの取得は面倒な作業ですが、jmastatsパッケージを使えば、気象データの取得・操作が簡単に行えます。この記事の中で扱ったヒートマップも、関心のある観測所のblock_noを指定するだけで簡単に作成できます。気象データを使った分析・可視化を行う際はjmastatsパッケージを是非ご利用ください。