「昔はこんなに暑くなかった」をR言語で可視化する

2023年の8月もそろそろ終わります。
しかし、まだまだ暑くて秋の気配はまだまだ来ないようです。

さてここ数年、7月に入ったころから「昔はこんなに暑くなかった」「いや、そんなことはない」といった話題がSNSをにぎわせています。
私も数年前にこんなグラフを作って、周りではそこそこ評判良かったです。

東京の8月の最低気温の分布の推移。いまさらだけど、ちゃんと書き直したのもあげておく。最低気温の上限が90年代以降上昇しているという結論は変わらないけど。boxplotもおまけで付けておいた。 pic.twitter.com/wfxgfr7I3O
— ボブさん (@bob3bob3) 2018年7月18日

今年もいろんな人が気温の可視化をしていて、非常に興味深かったです。

いくつか挙げてみましょう。

1) 東京における夏（6月～9月）の気温、過去148年分のヒートマップ。

東京における夏の気温を過去148年分ヒートマップで表現。
夏の気温は戦後徐々に上昇傾向にあり、現在では厳しい暑さが7月から9月中旬まで続くことがわかります。
データは気象庁「過去の気象データ検索」から取得。 pic.twitter.com/AdIm92rLKD
— 荻原和樹＠『データ思考入門』講談社現代新書 (@kaz_ogiwara) 2023年7月16日

これは数年前に評判になった記事のアップデート版ですね。
toyokeizai.net

2) 猛暑日の日数のカウント。
www.facebook.com
これはNHKの関東ローカル番組首都圏ネットワークで紹介されていたグラフですがシンプルでいて非常に分かりやすいです。

3) 基準値との偏差のグラフ
www.data.jma.go.jp
こちらは本家本元気象庁謹製のグラフ。強いです。

松本健太郎さんもこれをベースに独自のグラフを描いています。
comemo.nikkei.com

さて、どの可視化も素晴らしいですが、やっぱり自分で可視化してみたいよね、ということでこれらのグラフをR言語で再現してみます。
方針として、完コピではなくR言語らしい（ggplot2らしい）可視化を目指します。

0) 気温のデータを取得する。
1) 猛暑日の日数
2) 平均気温のヒートマップ
3)偏差のグラフ

0) 気温のデータを取得する。

可視化をする前に、大元となるデータを取得しなくてはいけません。

過去の気温のデータは気象庁が公開してくれているので、ありがたくスクレイピングさせてもらいます。

www.data.jma.go.jp

WebAPIで取得できれば良かったんですが、天気予報のWebAPIは見つけられたものの、過去データのものは発見できなくてスクレイピングすることにしました、

まずは使いたいパッケージを一通り呼び出す。
インストールされていなければ自動的にインストールされます。

#パッケージ管理用パッケージ
if (!require(pacman)) {install.packages("pacman")} 

#必要なパッケージがインストールされていなければインストールする
pacman::p_load( 
  conflicted, #関数の衝突防止
  tidyverse,  #モダンなデータ操作
  rvest,      #Webスクレイピング
  polite,     #紳士的なWebスクレイピング
  janitor,    #データクリーニング
  RcppRoll,   #移動平均算出
  ggridges    #リッジグラフ
  )

気象庁の過去の気温のページは地域ごと月ごとに分かれているので、それらを連続して取得し、さらに分析しやすい形に整形まで行う関数を作ります。

# スクレイピングから基本的なクリーニング、データ整形まで一連の流れを関数化
get_temp_data <- function(URL) {
  URL |>
    bow() |> # Introduce yourself to the host
    scrape(content = "text/html; charset=UTF-8") |> # 実際にhtmlファイルを取得
    html_table() |> # htmlのテーブルをdf形式に変換
    pluck(6) |> # 6個目のテーブルを抜き出す
    clean_names() |> # 列名を扱いやすい形に整える
    select(ri, qi_wen_c, qi_wen_c_2, qi_wen_c_3, shi_du) |> # 日付、平均気温、最高気温、最低気温、平均湿度
    slice(-(1:3)) |> # 余計な行を取り除く
    mutate( # 型の変換
      day = ri,
      mean_temp = qi_wen_c |> as.double(),
      max_temp = qi_wen_c_2 |> as.double(),
      min_temp = qi_wen_c_3 |> as.double(),
      mean_humidity = shi_du |> as.integer(),
      DI = 0.81 * mean_temp + 0.01 * mean_humidity * (0.99 * mean_temp - 14.3) + 46.3, #不快指数
      .keep = "none"
    )
}

# 複数年分のデータを一括取得し整形する関数
# base_urlにはhttpsからyear=までのurlを与える。
# year_rangeには取得対象の西暦をベクトルで与える。 
# month_rangeには取得対象月をベクトルで与える。
get_multiple_years_data <- function(base_url, year_range, month_range) {
  year_month_range <- crossing(year_range, month_range)
 
  paste0(base_url, year_month_range$year_range, "&month=", year_month_range$month_range, "&day=&view=p1") |>
    map(\(url) get_temp_data(url), .progress = TRUE) |>
    list_rbind(names_to = "rowid") |> # 年ごとに分かれたデータを一つにまとめる
    left_join( # 年を列を追加
      year_month_range |>
        rowid_to_column(),
      by = join_by(rowid)
    ) |>
    select(!rowid) |> # 不要になった列を削除
    mutate(age = paste0(floor(year_range / 10) * 10, "s")) |> #年代を追加
    rename(year = year_range, month =month_range)   #列名の調整
  }

これで準備できたので実際にデータを取得します。
今回は、札幌、秋田、東京、那覇の4地点の6月から9月までの気温などデータを取得してみます。
サーバーに負荷がかからないように時間をかけてデータ取得するため、1地点あたり2時間程度かかります。

念のため、データが取得できたらすぐに保存します。

東京：

#1875年から2022年まで
data_heatmap_tokyo1 <- get_multiple_years_data(
  base_url = "https://www.data.jma.go.jp/obd/stats/etrn/view/daily_s1.php?prec_no=44&block_no=47662&year=",
  year_range = 1875:2022,
  month_range = 6:9
  )

# 今年の8月まで
data_heatmap_tokyo2 <- get_multiple_years_data(
  base_url = "https://www.data.jma.go.jp/obd/stats/etrn/view/daily_s1.php?prec_no=44&block_no=47662&year=",
  year_range = 2023,
  month_range = 6:8
)

# がっちゃんこ
data_heatmap_tokyo <- bind_rows(
  data_heatmap_tokyo1,
  data_heatmap_tokyo2
) |>
  mutate(
  year,
  month = formatC(month, width = 2, flag = "0"),
  day = formatC(as.integer(day), width = 2, flag = "0")
  )

save(data_heatmap_tokyo, file = "data_heatmap_tokyo.RData")

札幌：

data_heatmap_sapporo1 <- get_multiple_years_data(
  base_url = "https://www.data.jma.go.jp/obd/stats/etrn/view/daily_s1.php?prec_no=14&block_no=47412&year=",
  year_range = 1879:2022,
  month_range = 6:9
)

data_heatmap_sapporo2 <- get_multiple_years_data(
  base_url = "https://www.data.jma.go.jp/obd/stats/etrn/view/daily_s1.php?prec_no=14&block_no=47412&year=",
  year_range = 2023,
  month_range = 6:8
)

# がっちゃんこ
data_heatmap_sapporo <- bind_rows(
  data_heatmap_sapporo1,
  data_heatmap_sapporo2
) |>
  mutate(
    year,
    month = formatC(month, width = 2, flag = "0"),
    day = formatC(as.integer(day), width = 2, flag = "0")
  )

save(data_heatmap_sapporo, file = "data_heatmap_sapporo.RData")

那覇：

data_heatmap_naha1 <- get_multiple_years_data(
  base_url = "https://www.data.jma.go.jp/obd/stats/etrn/view/daily_s1.php?prec_no=91&block_no=47936&year=",
  year_range = 1890:2022,
  month_range = 6:9
) |>
  mutate(
    year,
    month = formatC(month, width = 2, flag = "0"),
    day = formatC(as.integer(day), width = 2, flag = "0")
  )

data_heatmap_naha2 <- get_multiple_years_data(
  base_url = "https://www.data.jma.go.jp/obd/stats/etrn/view/daily_s1.php?prec_no=91&block_no=47936&year=",
  year_range = 2023,
  month_range = 6:8
) |>
  mutate(
    year,
    month = formatC(month, width = 2, flag = "0"),
    day = formatC(as.integer(day), width = 2, flag = "0")
  )

data_heatmap_naha <- bind_rows(
  data_heatmap_naha1,
  data_heatmap_naha2
) |>
  mutate(
    year,
    month = formatC(month, width = 2, flag = "0"),
    day = formatC(as.integer(day), width = 2, flag = "0")
  )

save(data_heatmap_naha, file = "data_heatmap_naha.RData")

秋田：

data_heatmap_akita1 <- get_multiple_years_data(
  base_url = "https://www.data.jma.go.jp/obd/stats/etrn/view/daily_s1.php?prec_no=32&block_no=47582&year=",
  year_range = 1883:2022,
  month_range = 6:9
) |>
  mutate(
    year,
    month = formatC(month, width = 2, flag = "0"),
    day = formatC(as.integer(day), width = 2, flag = "0")
  )

data_heatmap_akita2 <- get_multiple_years_data(
  base_url = "https://www.data.jma.go.jp/obd/stats/etrn/view/daily_s1.php?prec_no=32&block_no=47582&year=",
  year_range = 2023,
  month_range = 6:8
) |>
  mutate(
    year,
    month = formatC(month, width = 2, flag = "0"),
    day = formatC(as.integer(day), width = 2, flag = "0")
  )

data_heatmap_akita <- bind_rows(
  data_heatmap_akita1,
  data_heatmap_akita2
) |>
  mutate(
    year,
    month = formatC(month, width = 2, flag = "0"),
    day = formatC(as.integer(day), width = 2, flag = "0")
  )

save(data_heatmap_akita, file = "data_heatmap_akita.RData")

これでデータ取得完了です。

1) 猛暑日の日数

まずは地点ごとに猛暑日の日数をカウントしてみましょう。

# 保存したデータを呼び出す
attach("data_heatmap_tokyo.RData")
attach("data_heatmap_sapporo.RData")
attach("data_heatmap_naha.RData")
attach("data_heatmap_akita.RData")

# 猛暑日（最高気温35度以上）の日数を年ごとにカウントし棒グラフにする
bind_rows(
  data_heatmap_tokyo |> mutate(city = "東京"),
  data_heatmap_sapporo |> mutate(city = "札幌"),
  data_heatmap_naha |> mutate(city = "那覇"),
  data_heatmap_akita |> mutate(city = "秋田")
) |>
  mutate(
    extremely_hot_day = if_else(max_temp >= 35, 1, 0),
    city = city |>
      ordered(levels = c("東京", "札幌", "那覇", "秋田"))
  ) |>
  select(year, extremely_hot_day, city) |>
  summarise(extremely_hot_day = sum(extremely_hot_day, na.rm = TRUE), .by = c("year", "city")) |>
  ggplot(aes(x=year, y=extremely_hot_day, fill = extremely_hot_day)) +
  geom_col() +
  scale_fill_gradientn(colors = c("blue", "red")) +
  labs(
    title = "6月から9月の猛暑日（最高気温35度以上）の日数",
    subtitle = "2023年は8月20日までのデータ",
    x = "年",
    y = "猛暑日の日数",
    fill = "日数"
    ) +
  facet_wrap(vars(city), ncol=2) +
  theme(
    axis.text.x = element_text(angle = 270, hjust = 1)
  )

うおお、いきなり東京の増え方がヤバいですね。
秋田も北国のイメージですが意外と多い。
札幌と那覇が似たようなもの。
もしかして沖縄は避暑地になるのか？？？

せっかくなので真夏日（最高気温30度以上）も含めてカウントしてみましょう。

bind_rows(
  data_heatmap_tokyo |> mutate(city = "東京"),
  data_heatmap_sapporo |> mutate(city = "札幌"),
  data_heatmap_naha |> mutate(city = "那覇"),
  data_heatmap_akita |> mutate(city = "秋田")
) |>
  mutate(
    hot_day = if_else(max_temp >= 30, 1, 0),
    city = city |>
      ordered(levels = c("東京", "札幌", "那覇", "秋田"))
  ) |>
  select(year, hot_day, city) |>
  summarise(hot_day = sum(hot_day, na.rm = TRUE), .by = c("year", "city")) |>
  ggplot(aes(x=year, y=hot_day, fill = hot_day)) +
  geom_col() +
  scale_fill_gradientn(colors = c("blue", "red")) +
  labs(
    title = "6月から9月の真夏日と猛暑日（最高気温30度以上）の日数",
    subtitle = "2023年は8月20日までのデータ",
    x = "年",
    y = "猛暑日の日数",
    fill = "日数"
    ) +
  facet_wrap(vars(city), ncol=2) +
  theme(
    axis.text.x = element_text(angle = 270, hjust = 1)
    )

あ、やっぱり全体としては沖縄の方が暑いですね……

2) 平均気温のヒートマップ

お次は平均気温のヒートマップです。
これもいちいち長いコードを書かないで済むように一時的な描画用関数を作ります。
本当は関数内で最高気温の変数名を指定したかったんだけど、うまくできなかったのでtempで決め打ちしてます。

# ヒートマップ描画関数
temp_heatmap <- function(df, title) {
  df |>
    drop_na(temp) |>
    ggplot(aes(x = day, y = year, fill = temp)) +
    geom_tile(color="black") + 
    scale_y_reverse() +
    facet_grid(. ~ month) +
    theme(
      axis.text.x = element_text(angle = 90, hjust = 1, size = 3)
      ) +
    scale_fill_gradientn(colors = c("blue", "yellowgreen", "yellow","red")) +
    labs(title = title, x = "日", y = "年", fill = title)
}

東京：

data_heatmap_tokyo |>
  mutate(temp = mean_temp) |>
  temp_heatmap("東京の平均気温")

札幌：

data_heatmap_sapporo |>
  mutate(temp = mean_temp) |>
  temp_heatmap("札幌の平均気温")

秋田：

data_heatmap_akita |>
  mutate(temp = mean_temp) |>
  temp_heatmap("秋田の平均気温")

那覇：

data_heatmap_naha |>
  mutate(temp = mean_temp) |>
  temp_heatmap("那覇の平均気温")

那覇は戦争の期間のデータが欠測値になってますね……

全部まとめちゃいましょう。

data_all <- bind_rows(
  東京 = data_heatmap_tokyo,
  札幌 = data_heatmap_sapporo,
  那覇 = data_heatmap_naha,
  秋田 = data_heatmap_akita,
  .id = "city"
) |>
  mutate(
    temp = mean_temp,
    city = city |> factor(levels = c("札幌", "秋田", "東京", "那覇"))
    )

data_all |>
  drop_na(temp) |>
  ggplot(aes(x = day, y = year, fill = temp)) +
  geom_tile() + 
  scale_y_reverse() +
  facet_grid(city ~ month) +
  theme(
    axis.text.x = element_text(angle = 270, hjust = 1, size = 3)
    ) +
  scale_fill_gradientn(colors = c("blue", "yellowgreen", "yellow","red")) +
  labs(title = "平均気温", x = "日", y = "年", fill = "平均気温")

どこも年々気温が高くなり、夏の期間が長くなっていることが分かりますね。

3)偏差のグラフ

気象庁の偏差のグラフは1991年〜2020年の30年平均値を基準値としているので、ここではそれに倣って各地点の基準値との偏差を計算したうえで、5年移動平均値とそれに基づく平滑曲線を描いててみます。

これも地点ごとに簡単に書けるよう一時的な関数を作ります。

my_plot <- function(df) {
  reference_value <- df |>
    dplyr::filter(year >= 1991, year <= 2020) |>
    summarise(
      mean_max = mean(max_temp, na.rm = TRUE),
      mean_mean = mean(mean_temp, na.rm = TRUE),
      mean_min = mean(min_temp, na.rm = TRUE)
      )
  
  df |>
    summarise(
      max_temp = max_temp |> mean(na.rm = TRUE),
      mean_temp = mean_temp |> mean(na.rm = TRUE),
      min_temp = min_temp |> mean(na.rm = TRUE),
      .by = year
    ) |>
    mutate(
      year = year,
      diff_max = max_temp - reference_value$mean_max,
      diff_mean = mean_temp - reference_value$mean_mean,
      diff_min = min_temp - reference_value$mean_min,
      `最高気温の移動平均` = roll_mean(diff_max, n = 5, align = "right", fill = NA, na.rm = TRUE),
      `平均気温の移動平均` = roll_mean(diff_mean, n = 5, align = "right", fill = NA, na.rm = TRUE),
      `最低気温の移動平均` = roll_mean(diff_min, n = 5, align = "right", fill = NA, na.rm = TRUE),
      .keep = "none"
    ) |>
    select(year, `最高気温の移動平均`, `平均気温の移動平均`, `最低気温の移動平均`) |>
    pivot_longer(!year) |>
    ggplot(aes(x = year, y = value, color = name, group = name)) +
    geom_hline(yintercept = 0, color = "black", linetype = 2) +  
    geom_line(linewidth = 0.5) +
    geom_smooth(se = FALSE, linewidth = 1) +
    theme(legend.position = "bottom") +
    labs(
      subtitle = "偏差の5年移動平均値と平滑曲線", 
      x = "年", 
      y = "気温の偏差", 
      color = "気温の種類"
      )
}

東京：

data_heatmap_tokyo |>
  my_plot() +
  labs(title = "東京")

札幌：

data_heatmap_sapporo |>
  my_plot() +
  labs(title = "札幌")

秋田：

data_heatmap_akita |>
  my_plot() +
  labs(title = "秋田")

那覇：

data_heatmap_naha  |>
  my_plot() +
  labs(title = "那覇")

いずれの地点も1980年ごろから動きが変わっていることが分かります。

4)10年ごとの分布で見る
最後にオリジナルな可視化をやってみます。
平均値は分かりやすいのですが、やはり分布の形を見てみたいということで、10年区切りで分布の変化が分かるようなグラフを描いてみます。
ここでは8月の最低気温に絞ってみます。

まず、リッジプロット。

data_all |>
  dplyr::filter(month == "08") |>
  ggplot(aes(x = min_temp, y = age, group=age, fill = after_stat(x))) +
  geom_density_ridges_gradient() +
  scale_fill_viridis_c(option = "C")+
  theme_bw(base_size = 9) +
  theme(legend.position = "bottom") +
  scale_x_continuous(limits = c(6, 34)) +
  facet_wrap(vars(city), ncol=4) +
  labs(title = "8月の最低気温の分布", x = "最低気温", y = "年代", fill = "最低気温")

そしてバイオリンプロットと箱ひげ図の組み合わせ。

data_all |>
  dplyr::filter(month == "08") |>
  ggplot(aes(x = age, y = min_temp)) +
  geom_violin(fill = "orange") +
  geom_boxplot(width = 0.3, fill = "lightblue", outlier.shape = NA) +
  theme(
    text = element_text(size = 9),
    axis.text.x = element_text(angle = 270, hjust = 1)
    ) +
  facet_wrap(vars(city), ncol=2) +
  labs(title = "8月の最低気温の分布", x = "年代", y = "最低気温")