まだ厨二病

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

🌎全球規模での気候データをRからサクッと得たい(可視化もするよ)

昨日書いた記事で使ったRパッケージの{rWBclimate}、これは良いものだ、と思えたので別途こうして使い方をメモしておく。

近い将来に予測される大規模な気候変動や乾燥化などを扱った論文が増えているように、気象・気候データは各種の研究やデータ解析において重要となる。局所的なものは自前で気象ステーションとかを設置したりする必要があるのだろうけど、傾向を把握するためにざっくりとした広範囲のデータが欲しい時もある。あるいは日本🇯🇵とか世界全体🌎とかを対象にして議論をしたい時などなど...。

{rWBclimate}はそんな世界中の気象・気候データをR上で取得し、扱うためのパッケージである。開発は信頼と実績のROpenSciが行っている。世界銀行 http://www.worldbank.org が提供している開発者向けのウェブAPIの一種であるClimate Data APIを利用することで、気候変動に関する政府間パネルで利用されている全球的な大気循環モデル global circulation model (GCM)の15のシナリオからなるモデリングの値や過去の気象データがサクッと扱えるようになっている。以降の可視化やモデリングはRでやることになるが、こうしたデータが簡単に手に入るのは魅力的だ。

なお、こうした時系列データの解析にはちょっと注意が必要なので次の記事も参考にしてもらいたい

uribo.hatenablog.com

🎓 GCMsの前知識

ちょっと寄り道して{rWBclimate}が利用する大気循環モデルの概要について整理しておく。すでに理解のある人はすっ飛ばしてもらって構わない。詳細はGitHubREADMEにもある。またAPIの仕様はドキュメントを見るのがてっとり早い。

  • モデルのデータ
  • 大気循環モデルではそのシナリオ、対象の国や地域ごとに予測値が異なってくる。
  • データは気温と降水量について
    • データには月および年ごとの平均値と過去の値から予測されうる異常値がある
  • A2とB1シナリオ(温室効果ガスなどの排出量を抑えた生態学的に好ましいシナリオ)
  • 時間スケール
    • 1920年からの過去と2099年までの未来について、20年間ごとに与えられる
    • 過去 1920-1939, 1940-1959, 1960-1979, 1980-1999
    • 将来 2020-2039, 2040-2059, 2060-2079, 2080-2099
  • 空間スケール
    • 2種類の空間スケールを扱う
    • 国単位... 大文字のアルファベッド3文字からなるISOコード。ここを参考にすると良い
    • 大陸ごとの流域単位

🔰 使い方

パッケージはRのパッケージ(拡張機能)などを管理するCRANに登録されている。もしインストールされていなければinstall.packages()関数を使ってローカルパッケージに追加してくれば良い。

# install.packages('rWBclimate')
library(rWBclimate)

モデルデータのダウンロード

get_model_precip()get_model_temp()関数によってモデルデータを取得する。引数で対象の地域 locatorとデータの種類type、取得年の範囲を指定する。

# 日本の降水量のデータ(過去と将来について)を得る
df_jp_precip <- get_model_precip(locator = "JPN", 
    type = "mavg", start = 1980, end = 2020)
# 中身はこんな感じ
df_jp_precip %>% glimpse()
## Observations: 516
## Variables: 7
## $ scenario (fctr) past, past, past, past, past, past, past, past, past...
## $ fromYear (dbl) 1980, 1980, 1980, 1980, 1980, 1980, 1980, 1980, 1980,...
## $ toYear   (dbl) 1999, 1999, 1999, 1999, 1999, 1999, 1999, 1999, 1999,...
## $ gcm      (fctr) bccr_bcm2_0, bccr_bcm2_0, bccr_bcm2_0, bccr_bcm2_0, ...
## $ data     (dbl) 147.78167, 127.11897, 125.91752, 128.32521, 129.09497...
## $ month    (int) 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 1, 2, 3, 4, 5,...
## $ locator  (chr) "JPN", "JPN", "JPN", "JPN", "JPN", "JPN", "JPN", "JPN...
df_jp_precip$gcm %>% levels()
##  [1] "bccr_bcm2_0"     "cccma_cgcm3_1"   "cnrm_cm3"       
##  [4] "csiro_mk3_5"     "gfdl_cm2_0"      "gfdl_cm2_1"     
##  [7] "ingv_echam4"     "inmcm3_0"        "ipsl_cm4"       
## [10] "miroc3_2_medres" "miub_echo_g"     "mpi_echam5"     
## [13] "mri_cgcm2_3_2a"  "ukmo_hadcm3"     "ukmo_hadgem1"
# 各種データフレーム操作のためのパッケージによる処理が行える
# 気象庁の気象研究所によるMRI-CGCM2.3.2シナリオに絞る
df_jp_precip %<>% dplyr::filter(gcm %in% 
    c("mri_cgcm2_3_2a")) %>% dplyr::select(-c(locator, 
    gcm))
df_jp_precip %>% ggplot(aes(x = as.factor(month), 
    y = data, colour = scenario, group = scenario)) + 
    geom_point() + geom_path() + ylab("Average precipitation in degrees C \n between 2020 and 2040") + 
    xlab("Month") + theme_bw()

f:id:u_ribo:20160203102430p:plain

過去の気候データ(モデル値)

過去の気候データは、月、年単位でも得ることができる。シナリオによって与えられた値で実際の観測値ではないということに注意が必要。

# 日本、バヌアツ、ドイツの年平均気温を取得
df_temp <- get_historical_temp(locator = c("JPN", 
    "VUT", "DEU"), time_scale = "year") %>% 
    dplyr::filter(year >= 1980)
df_temp %>% ggplot(aes(x = year, y = data, 
    group = locator, colour = locator)) + 
    geom_point() + geom_path() + ylab("Average historical temperature (cent.)") + 
    xlab("Year") + theme_bw()

f:id:u_ribo:20160203102439p:plain

地図へのマッピング

{rWBclimate}では取得した気候データを地図上にマッピングするためのKMLファイルをダウンロードしてRに読み込む関数があり、地図を描くために必要なkmlファイルをダウンロードするためのパスを設定しておく必要がある。保存したくない場合はtempdir()関数とかを使うと良いだろう。ここではoptions()関数によって引数kmlpathの値がそれになる。都度入力が面倒な場合は.Rprofileに環境変数として書いておけば良い。環境変数の設定はこちらを参考にして欲しい。

uribo.hatenablog.com

# KMLファイル用のパスを与えておく
options(kmlpath = paste(Sys.getenv("HOME"), 
    "Dropbox/maps", "kml", sep = "/"))
# get_ensemble_temp()で気温データを取得する
# アジア地域の年平均気温(ちょっと時間がかかる)
df_asia_precip <- get_ensemble_temp(locator = Asia_country, 
    type = "annualavg", start = 2020, end = 2040) %>% 
    dplyr::filter(scenario == "b1", fromYear == 
        2020, percentile == 50)
# KMLファイルをRオブジェクトとして読み込む。必要なKMLファイルがあればダウンロードしてくる
df_asia_contry <- create_map_df(Asia_country)
df_asia_contry %>% glimpse()
## Observations: 46,744
## Variables: 8
## $ long  (dbl) 67.77988, 67.77422, 67.78693, 68.00139, 68.05801, 68.182...
## $ lat   (dbl) 37.18582, 37.11555, 37.09221, 36.93610, 36.93253, 37.017...
## $ order (int) 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 1...
## $ hole  (lgl) FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, FALSE, ...
## $ piece (fctr) 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1...
## $ id    (chr) "0", "0", "0", "0", "0", "0", "0", "0", "0", "0", "0", "...
## $ group (chr) "0.1-AFG", "0.1-AFG", "0.1-AFG", "0.1-AFG", "0.1-AFG", "...
## $ ID    (chr) "AFG", "AFG", "AFG", "AFG", "AFG", "AFG", "AFG", "AFG", ...
# 地図を描画するためのデータフレームを作成する。引数return_map
# = TRUEでプロットを行う
climate_map(map_df = df_asia_contry, data_df = df_asia_precip, 
    return_map = FALSE)

マッピングの自由度を上げるためには{ggplot2}や関連のパッケージの関数を使ったほうが良い。すんなりとオブジェクトの引き渡しをしてくれる。

climate_map(map_df = df_asia_contry, data_df = df_asia_precip) + 
    ggtitle("Predict Annual Temperature\n between 2020 and 2040") + 
    scale_fill_continuous("Temperature", 
        low = "Yellow", high = "red") + coord_map(projection = "mercator") + 
    ggthemes::theme_map()

f:id:u_ribo:20160203102520p:plain

📜 参考

💻 実行環境

package * version date source
dplyr * 0.4.3.9000 2015-10-28 Github (hadley/dplyr@52d10f6)
ggplot2 * 2.0.0 2015-12-18 CRAN (R 3.2.3)
magrittr * 1.5 2016-01-13 Github (smbache/magrittr@00a1fe3)
rWBclimate * 0.1.3 2014-04-19 CRAN (R 3.2.0)