Ｒで気候変動に備える

Who

会場の近くで農業分野での気象データ利用に関する研究をしています
元々は植物いじりタイプで気象データいじりのビギナーです
生命科学系と比較して気象系はRユーザーが少なくて肩身が狭いです
.Rに初参加、お手柔らかにおねがいします

農業 = funct(天気)

一部の例外を除けば農業は天気に強く依存する

気象データはとっつきにくい

いい感じのやつ、あります

農研機構メッシュ農業気象データシステム
- 気象庁のアメダス観測値と格子状の予報値を統合
- 同一の空間分解能で過去と未来の横断的に解析
- 研究・開発・教育・試用なら無料
- Rインターフェイスあり (要認証)

デモ | パッケージ読み込み

suppressPackageStartupMessages(library(tidyverse))
suppressPackageStartupMessages(library(lubridate))

# devtools::install_github("KeachMurakami/agrmesh")
library(agrmesh)

ℹ WELCOME to R-AMGSDS interface (ver.0.0.1.4) 
ℹ 農研機構は、農業分野や他の分野における研究・開発・教育・試用を目的とする者に、審査に基づきメッシュ農業気象データ（以下、「このデータ」と呼ぶ。）の利用を許可します。
ℹ 特に許可されない限り、このデータを他に転載したり第三者に提供したりすることはできません。
ℹ このデータを利用して作成した情報を販売することはできません。
ℹ 利用者は、利用期間の終了後、速やかに利用報告をすることとします。
ℹ 農研機構は、利用者がこのデータの利用によって生じた結果、ならびに、このデータが利用できないことによって生じた結果について、いかなる責任も負いません。
ℹ このデータを利用して得た成果等を公表する場合は、「農研機構メッシュ農業気象データ（The Agro-Meteorological Grid Square Data, NARO）」を利用した旨を明記してください。
✔ USERID and PASSWORD -> verified

デモ | 積雪深マップ

# 会場の緯度経度
site_lat <-  43.05
site_lon <- 141.37

# 会場周辺の積雪深 (SD) の分布を取得
snowdepth <-
  fetch_amgsds(times = ymd("2024-02-17"), elements = "SD", mode = "area",
               lats = site_lat + c(-0.5, 0.5), lons = site_lon + c(-0.4, 0.4))

print(snowdepth)

# A tibble: 6,130 × 4
   time         lat   lon    SD
   <date>     <dbl> <dbl> <dbl>
 1 2024-02-17  42.6  141.    45
 2 2024-02-17  42.6  141.    43
 3 2024-02-17  42.6  141.    44
 4 2024-02-17  42.6  141.    41
 5 2024-02-17  42.6  141.    43
 6 2024-02-17  42.6  141.    47
 7 2024-02-17  42.6  141.    41
 8 2024-02-17  42.6  141.    41
 9 2024-02-17  42.6  141.    44
10 2024-02-17  42.6  141.    39
# ℹ 6,120 more rows

デモ | 積雪深マップ

# 会場周辺の積雪深 (SD) の分布を可視化
snowdepth |>
  ggplot(aes(lon, lat, fill = SD)) +
  geom_tile() +
  coord_map() +
  scale_fill_viridis_c()

デモ | 将来気温

50年後のこのあたりの気温推移を見る
- 気候モデルMIROC5 x 高排出シナリオ (RCP8.5)

# 2071年のモデルベースの気温を取得
temp_future <-
  fetch_amgsds(times = ymd("2071-05-01", "2071-09-30"), elements = c("TMP_mea", "TMP_max"),
               lats = site_lat, lons = site_lon,
               source = "scenario", model = "MIROC5", RCP = "RCP8.5")

print(temp_future)

# A tibble: 153 × 6
   time         lat   lon site_id TMP_mea TMP_max
   <date>     <dbl> <dbl> <chr>     <dbl>   <dbl>
 1 2071-05-01  43.1  141. 1          13.8    20.2
 2 2071-05-02  43.1  141. 1          16.6    23.8
 3 2071-05-03  43.1  141. 1          16.5    21.5
 4 2071-05-04  43.1  141. 1          14.5    18.8
 5 2071-05-05  43.1  141. 1          12.3    18.2
 6 2071-05-06  43.1  141. 1          12.6    17.2
 7 2071-05-07  43.1  141. 1          16.3    21.5
 8 2071-05-08  43.1  141. 1          15.7    22.2
 9 2071-05-09  43.1  141. 1          15.1    18.6
10 2071-05-10  43.1  141. 1          14.8    17.1
# ℹ 143 more rows

デモ | 将来気温

観測値ベースの2011年の値と比べると…

# 2011年の観測値ベースの気温を取得
temp_present <-
  fetch_amgsds(times = ymd("2011-05-01", "2011-09-30"), elements = c("TMP_mea", "TMP_max"),
               lats = site_lat, lons = site_lon)

# データを結合
temps <-
  bind_rows(present = temp_present, future = temp_future, .id = "period") |>
  mutate(date = yday(time) + ymd("2099-12-31")) # 年の異なるデータを並べて表示するための処理

print(temps)

# A tibble: 306 × 8
   period  time         lat   lon site_id TMP_mea TMP_max date      
   <chr>   <date>     <dbl> <dbl> <chr>     <dbl>   <dbl> <date>    
 1 present 2011-05-01  43.1  141. 1         10.8    13.5  2100-05-01
 2 present 2011-05-02  43.1  141. 1          6.27   10.4  2100-05-02
 3 present 2011-05-03  43.1  141. 1          6.69   10.0  2100-05-03
 4 present 2011-05-04  43.1  141. 1          5.35    6.88 2100-05-04
 5 present 2011-05-05  43.1  141. 1          7.15   11.3  2100-05-05
 6 present 2011-05-06  43.1  141. 1          9.30   15.4  2100-05-06
 7 present 2011-05-07  43.1  141. 1         10.8    13.5  2100-05-07
 8 present 2011-05-08  43.1  141. 1          9.44   12.9  2100-05-08
 9 present 2011-05-09  43.1  141. 1         10.3    14.0  2100-05-09
10 present 2011-05-10  43.1  141. 1          9.27   13.6  2100-05-10
# ℹ 296 more rows

デモ | 将来気温

観測値ベースの2011年の値と比べると…

temps |>
  ggplot(aes(date, TMP_max, col = period)) +
  geom_line() +
  geom_smooth()

`geom_smooth()` using method = 'loess' and formula = 'y ~ x'

デモ | イネっぽい作物の収穫期

定植 (5月1日) からの有効積算気温が2200 °C・日で収穫

※ 日付・数字は計算の一例で実態には即しません

\[有効積算気温 = \sum_i{\max(日平均気温_{\textrm{day } i}, 0)}\]

デモ | イネっぽい作物の収穫期

temps |>
  group_by(period) |>
  mutate(TMP_cumsum = cumsum(pmax(0, TMP_mea)),
         harvest = sum(TMP_cumsum < 2200) + min(date)) |>
  ggplot(aes(date, TMP_cumsum, col = period)) +
  geom_line() +
  geom_hline(yintercept = 2200, linetype = 2) +
  geom_vline(aes(xintercept = harvest), linetype = 2)

不確実な未来に備える

来週の天気もわからないのに50年後の天気なんて…

アンサンブル予報

気象庁のwebページを引用

数値予報結果の誤差の原因は大きく2つに分けられ、ひとつは初期値に含まれる誤差が拡大すること、もうひとつは数値予報モデルが完全ではないことです。数値予報では、わずかに異なる2つの初期値から予報した2つの予報結果は、初めのうち互いによく似ていますが、その差は時間の経過とともに拡大します。これは、大気の運動にある特徴的な性質「初期値の小さな差が将来大きく増大する」というカオス（混沌）的な振る舞いのひとつです。実際の数値予報では、観測データの誤差や解析手法の限界から、初期値に含まれる誤差をゼロにすることはできず、時間とともに誤差が拡大することを避けることはできません。また、数値予報では、計算機の性能の限界により、ある大きさの格子を用いた近似式で気温や風等の予測計算を行います。このように近似式を使っていることからも、予報結果に誤差が生じます。

このような誤差の拡大を事前に把握するため、「アンサンブル（集団）予報」という数値予報の手法を利用しています。この手法では、ある時刻に少しずつ異なる初期値を多数用意するなどして多数の予報を行い、平均やばらつきの程度といった統計的な情報を用いて気象現象の発生を確率的に捉えることが可能となります。

アンサンブル予報

気象庁のwebページを引用

数値予報結果の誤差の原因は大きく2つに分けられ、ひとつは初期値に含まれる誤差が拡大すること、もうひとつは数値予報モデルが完全ではないことです。数値予報では、わずかに異なる2つの初期値から予報した2つの予報結果は、初めのうち互いによく似ていますが、その差は時間の経過とともに拡大します。これは、大気の運動にある特徴的な性質「初期値の小さな差が将来大きく増大する」というカオス（混沌）的な振る舞いのひとつです。実際の数値予報では、観測データの誤差や解析手法の限界から、初期値に含まれる誤差をゼロにすることはできず、時間とともに誤差が拡大することを避けることはできません。また、数値予報では、計算機の性能の限界により、ある大きさの格子を用いた近似式で気温や風等の予測計算を行います。このように近似式を使っていることからも、予報結果に誤差が生じます。

このような誤差の拡大を事前に把握するため、「アンサンブル（集団）予報」という数値予報の手法を利用しています。この手法では、ある時刻に少しずつ異なる初期値を多数用意するなどして多数の予報を行い、平均やばらつきの程度といった統計的な情報を用いて気象現象の発生を確率的に捉えることが可能となります。

アンサンブル予報

\[y_t = 1.1 \times y_{t-1} + \epsilon\]

sim_chaos <-
  function(y0, times = 30){
    noise <- rnorm(times, 0, .1)
    y <- y0
    
    for(i in 2:times){
      y[i] <- 1.1 * y[i-1] + noise[i-1]
    }
    
    tibble(time = 1:times, y)
  }

set.seed(20240217)

runif(10, min = 0.99, max = 1.01) |>
  map_dfr(.id = "id", sim_chaos) |>
  ggplot(aes(time, y, col = id, group = id)) +
  geom_line(show.legend = FALSE)

アンサンブル予報

個々のメンバーは

「将来気候のもと生起しうる気象」からの抽出
(十分なスピンアップを経れば) 互いに独立

「50年に一度の…」のような確率的な影響評価

アンサンブルに触れる

地球温暖化対策に資するアンサンブル気候予報データベース (d4PDF)
- 気候モデル別・CO2排出シナリオ別に最大100メンバで60年分の計算
- データはgrib形式 (GRIdded Binary) またはGrADS向けの形式 (.ctl + .dat)
  - grib形式: rNOMADSパッケージで扱うことができる
  - GrADS形式: Rから直接扱うことはできない (たぶん)

TL;DRな処理

過去実験・昇温実験の日平均気温をGrADS向けの形式でダウンロードする
Rで読み込みやすいNetCDF形式に変換する
tidyncパッケージで読み込み
下処理して保存

# https://d4pdf.diasjp.net/d4PDF.cgiより全球実験の日平均気温 (高度 2 m; sfc_avr_day TA) のアンサンブルデータをダウンロード
# 同ページで提供されている切り出し用Pythonスクリプト (d4pdf-sb.py) を使用
# 昇温実験 (HFB_4K_**) については6気候モデルの15アンサンブルの各1年分

YOUR_DIAS_PW <- "XXX"
YOUR_DIAS_ID <- "XXX"

# 過去実験 (HPB) の90アンサンブルから1年ずつの気温推移を取得
walk(1:90, function(i){
  system(
    str_glue("
      echo {YOUR_DIAS_PW} | python3 d4pdf-sb.py -u {YOUR_DIAS_ID} -o hpb_{i}.tar -f 2010-01 -t 2010-12 -e {i} -W 141 -E 142 -S 42 -N 43 GCM/HPB/sfc_avr_day TA
    ")
  )
})

# 6気候モデルでの昇温実験 (HFB_4K_**) の各15アンサンブルから1年ずつの気温推移を取得
walk(101:115, function(i){
  system(
    str_glue("
      echo {YOUR_DIAS_PW} | python3 d4pdf-sb.py -u {YOUR_DIAS_ID} -o hfb_cc_{i}.tar -f 2071-01 -t 2071-12 -e {i} -W 141 -E 142 -S 42 -N 43 GCM/HFB_4K_CC/sfc_avr_day TA
      echo {YOUR_DIAS_PW} | python3 d4pdf-sb.py -u {YOUR_DIAS_ID} -o hfb_ha_{i}.tar -f 2071-01 -t 2071-12 -e {i} -W 141 -E 142 -S 42 -N 43 GCM/HFB_4K_HA/sfc_avr_day TA
      echo {YOUR_DIAS_PW} | python3 d4pdf-sb.py -u {YOUR_DIAS_ID} -o hfb_gf_{i}.tar -f 2071-01 -t 2071-12 -e {i} -W 141 -E 142 -S 42 -N 43 GCM/HFB_4K_GF/sfc_avr_day TA
      echo {YOUR_DIAS_PW} | python3 d4pdf-sb.py -u {YOUR_DIAS_ID} -o hfb_mi_{i}.tar -f 2071-01 -t 2071-12 -e {i} -W 141 -E 142 -S 42 -N 43 GCM/HFB_4K_MI/sfc_avr_day TA
      echo {YOUR_DIAS_PW} | python3 d4pdf-sb.py -u {YOUR_DIAS_ID} -o hfb_mp_{i}.tar -f 2071-01 -t 2071-12 -e {i} -W 141 -E 142 -S 42 -N 43 GCM/HFB_4K_MP/sfc_avr_day TA
      echo {YOUR_DIAS_PW} | python3 d4pdf-sb.py -u {YOUR_DIAS_ID} -o hfb_mr_{i}.tar -f 2071-01 -t 2071-12 -e {i} -W 141 -E 142 -S 42 -N 43 GCM/HFB_4K_MR/sfc_avr_day TA
    ")
  )
})


# ダウンロードされたtarファイルを解凍してまとめる
dir(pattern = ".tar$") |>
  walk(untar, exdir = "ensembles")

# GrADS -> NetCDF変換関数
# r-bloggersから拝借、元のブログはリンク切れに...
# https://www.r-bloggers.com/2012/01/converting-grads-files-to-netcdf-from-within-r-using-cdo/
# CDO (Climate Data Operators) というコマンドラインツールのインストールが必要
grads2nc <-
  function(ctl_input, nc_output, verbose) {  
    if(Sys.which("cdo")[[1]] == ""){
      stop("cdo executable not found. Check PATH or install cdo.") 
    }
    if(verbose){
      ext = ""
    } else {
      ext = "> /dev/null"
    }
    cmd = sprintf("cdo -f nc import_binary %s %s %s", ctl_input, nc_output, ext)
    t = system.time(ret_code <- system(cmd))
    if(verbose & (ret_code != 127)){
      message(sprintf("Processing \"%s\" to \"%s\" is done, time spent:", ctl_input, nc_output))
      print(t)
    } else if(ret_code == 127){
      stop("Conversion of grads file to nc failed")
    }
    return(invisible(ret_code))
  }

# NetCDF -> tibble変換関数
# grads2ncで出力されたNetCDFを一時ファイルとして受け取り、tidyncで読み込む
read_grads <-
  function(input, verbose = FALSE){
    tmp <- tempfile()
    grads2nc(input, nc_output = tmp, verbose = verbose)
    
    # NetCDFを開いて緯度経度範囲指定で１グリッドの値を抽出する
    output <-
      tidync::tidync(tmp) |>
      tidync::hyper_tibble(lon = between(lon, 141.6, 142),
                           lat = between(lat, 42.8, 43.2))
    
    file.remove(tmp)
    return(output)
  }

# GrADSのファイルパスを取得し、文字列から各種情報を抜き出す
ctl_files <- dir("ensembles/", pattern = "ctl", recursive = TRUE, full.names = TRUE)
model <- str_extract(ctl_files, "HPB|HFB_4K_..")
member <- str_extract(ctl_files, "m[:digit:]{3}")
yyyymm <- str_extract(ctl_files, "[:digit:]{6}")

# GrADSファイルを一括で読み込んでtibble形式で整理する
temp_ensembles <-
  ctl_files |>
  map(read_grads) |>
  list(model, member, yyyymm) |>
  pmap_dfr(function(x1, x2, x3, x4){
    x1 |>
    transmute(
      model = x2,
      member = x3,
      date = ymd(paste0(x4, time %/% 24 + 1)), # 時単位 -> 日単位
      year = year(date),
      TMP_mea = ta - 273.15, # 絶対温度 -> 摂氏
      period = if_else(model == "HPB", "present", "future"),
      dummy_date = yday(date) + ymd("2999-12-31") # 異なる年のデータを並べて表示するための処理)
    )
  })

# データ保存@dropbox
write_csv(temp_ensembles, "temp_ensembles.csv")

アンサンブルに触れる

過去実験: 90年分 (90メンバー)
RCP8.5昇温実験: 90年分 (6気候モデル x 15メンバー)

temp_ensembles <- read_csv("https://www.dropbox.com/scl/fi/wc91b6aif5fxn8gvl2kik/temp_ensembles.csv?rlkey=e0c3f4qgkcwacplm0ewzn6r9y&dl=1")

Rows: 65700 Columns: 7
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr  (3): model, member, period
dbl  (2): year, TMP_mea
date (2): date, dummy_date

ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.

head(temp_ensembles)

# A tibble: 6 × 7
  model     member date        year TMP_mea period dummy_date
  <chr>     <chr>  <date>     <dbl>   <dbl> <chr>  <date>    
1 HFB_4K_CC m101   2071-01-01  2071  -3.15  future 3000-01-01
2 HFB_4K_CC m101   2071-01-02  2071  -2.23  future 3000-01-02
3 HFB_4K_CC m101   2071-01-03  2071  -5.68  future 3000-01-03
4 HFB_4K_CC m101   2071-01-04  2071  -4.34  future 3000-01-04
5 HFB_4K_CC m101   2071-01-05  2071   1.23  future 3000-01-05
6 HFB_4K_CC m101   2071-01-06  2071   0.603 future 3000-01-06

temp_ensembles |>
  select(model, period) |>
  table()

           period
model       future present
  HFB_4K_CC   5475       0
  HFB_4K_GF   5475       0
  HFB_4K_HA   5475       0
  HFB_4K_MI   5475       0
  HFB_4K_MP   5475       0
  HFB_4K_MR   5475       0
  HPB            0   32850

アンサンブルに触れる

単年度では確信を持てないトレンドがくっきり

temp_ensembles |>
  ggplot(aes(dummy_date, TMP_mea, col = model, group = paste(model, member, year))) +
  geom_line(linewidth = .1, alpha = .1) +
  geom_line(data = \(x) summarise(x, TMP_mea = mean(TMP_mea), .by = c(dummy_date, model)), aes(group = model)) +
  scale_x_date(date_labels = "%b", date_breaks = "1 month", name = NULL)

アンサンブルに触れる | 作物収穫期

アンサンブルでみると積算温度のばらつきがわかる

cumtemp_ensembles <-
  temp_ensembles |>
  filter(month(date) %in% 5:9) |>
  group_by(model, member, year) |>
  mutate(TMP_cumsum = cumsum(pmax(0, TMP_mea)),
         harvest = sum(TMP_cumsum < 2200) + min(dummy_date)) 

cumtemp_ensembles |>
  group_by(period, member, model, year) |>
  ggplot(aes(dummy_date, TMP_cumsum, col = model, group = paste(model, member, year))) +
  geom_line(linewidth = .5, alpha = .5) +
  geom_hline(yintercept = 2200) +
  annotate("rect", fill = NA, col = "black", linetype = 2,
           xmin = ymd("3000-08-01"), xmax = ymd("3000-09-30"), ymin = 2100, ymax = 2300)

アンサンブルに触れる | 作物収穫期

収穫期の変化を確率的に議論できる
- 「現在の90年に1回の高温年は将来のXパーセンタイルに相当」

cumtemp_ensembles |>
  distinct(period, model, member, year, harvest) |>
  ggplot(aes(harvest)) +
  geom_histogram(aes(y = after_stat(count), fill = model), alpha = .3, col = "white", binwidth = 1) +
  geom_density(aes(y = after_stat(density) * 90, group = period, col = period))

雑感

気象データの周辺環境はＲにあまり優しくない？
- ガチ勢は未だFortran、周辺勢はPython、最近はJuliaも
- 差分方程式的なループ計算だとＲの弱点が目立つ
解析の下処理と本番の比重によっては他言語推奨かも
- 自分はtidyverseの引力から脱出できません

ENJOY!

参考:

謝辞: 本研究はJSPS科研費 JP23K14051の助成を受けたものです。