気候予測データの解析環境を構築する　番外編― ”wxbcgribx.py”を動かす

データ解析環境について

はじめに

「気象ビジネス推進コンソーシアム（略称：WXBC）」により、 気象庁のGPVデータを処理し活用する方法を学ぶため、気象庁GPVデータ分析チャレンジ！入門（2023年10月26日）が開催されています。（現在受講中）

www.wxbc.jp

事前準備として、” WXBCオリジナルPythonライブラリwxbcgribX”の実行環境を作っておく必要があるとのこと。

www.wxbc.jp

試行錯誤してWSL2-Docker下で動くようにしましたので、紹介します。

wgrib2のコマンドパスの設定

wxbcgribx.pyの36行目に次のパスを設定します

wgrib2 = "docker exec -w /mnt/c/Users/hoge/{作業ディレクトリ} wgrib2_dock wgrib2"
#wgrib2 = Path("C:/wgrib2/wgrib2.exe") # wgrib2のインストール場所　Windowsの場合
#wgrib2 = Path("/usr/local/bin/wgrib2.exe") # wgrib2のインストール場所　Linuxの場合
#wgrib2 = Path("~/work/grib2/wgrib2/wgrib2") # wgrib2のインストール場所　Macの場合

一時フォルダの設定

私の環境では、tempfile.TemporaryDirectory()が動かないので、作業ディレクトリに”tmp"フォルダを作成して、そこを使うように設定します。

def from_grb(grbpath,matchopt,verbose):
    """
    gpvデータをGRIB2ファイルから取得する関数
    matchopt: データを選別するオプション文字列
    """
    #指定した気象要素のデータを取り出して一時ファイルに格納
    with tempfile.TemporaryDirectory() as td:
        Path("./tmp").mkdir(exist_ok = True)  #作業フォルダの作成
        #dumppath = Path(td)/"extracting.nc" #作業ファイルの指定をコメントアウト
        dumppath = "./tmp/extracting.nc" #作業ファイルの指定を修正
        opt = f'{matchopt} -netcdf {dumppath}'
        rc = wg2(grbpath, opt)
        if verbose:
            print(rc.stdout)#.splitlines())
        #一時ファイルからデータセットをロード
        with xr.open_dataset(dumppath) as ds:
            ds.load()
        #dumppath.unlink(missing_ok=True) #後片付けをコメントアウト
    return ds

事前テストのスクリプトが動きました。

これまで構築した、「気候予測データの解析環境」で”wxbcgribx.py”を動かす方法を紹介しました。試してみてくださいね～。

追記：GoogleColabでwgrib2を使うには

講義で出ていたこの質問、

質問です。どうしても、社内では供与PCへのインストールが出来ない環境にあります。GoogleColabでwgrib2を使うにはどうすればよいでしょうか？

調べてみたらこんな記事がありました。

qiita.com

↑この場合、wxbcgribx.pyにおけるwgrib2のコマンドパスは

wgrib2 = "wgrib2"　

でOKだと思います。

ご参考に～

追記2: 「気象庁GPVデータ分析チャレンジ！基礎編」の”wxbcgribx.py”を動かす

本日（2023.11.30）開催の「気象庁GPVデータ分析チャレンジ！基礎編」で配布された”wxbcgribx.py”は、バージョンアップされているようでうまく動かず・・・

73行目からの"wg2"で実行されているエラーチェックをコメントアウトすると動きました。

def wg2(src,opt=''):
    # wgrib2を実行する関数
    #if not wgrib2.exists():
    #    print("wgrib2 をインストールしてください。")
    #    sys.exit()
    check_path(src)
    rc = subprocess.run(f'{wgrib2} {opt} {src}', 
                        shell=True, text=True, capture_output=True)
#    print(rc.returncode)
#    print(rc.stdout)
#    print(rc.stderr)
    return rc

気候予測データの解析環境を構築する。その7―MCMCサンプラーの実行環境を構築する。

データ解析環境について

はじめに

気候予測データは、将来の平均的な「気候」状態を予測するためのもので、ある地域の気温・降水量の平均値や変動の標準偏差などの統計量がどう変化するか？を予測しています。(ココが知りたい地球温暖化より)

www.cger.nies.go.jp

気温など、連続した数値でかつ正規分布に従うことがわかっていれば、平均値や分散を求めることは簡単です。しかし、ガンマ分布に従う降水量や離散値である猛暑日日数の変化を単純な算術平均で表すことには問題があります。そのようなデータを扱うには"正規分布以外の確率分布を扱える統計モデル"が有用です。

“階層ベイズモデル”は、様々な確率分布を扱い複雑な統計モデルを構築することができ、マルコフ連鎖モンテカルロ法を用いてサンプリングをすることでパラメータを推計することができます。

Pythonでマルコフ連鎖モンテカルロ法を使うためのライブラリ（MCMCサンプラー）には、PyMC、Pyro, NumPyro, TensorFlowProbability, PyStan, PyJAGSなどがありますが、本研究室では"PyMC"を使っています。

www.pymc.io

ということで、PyMCの実行環境を構築していきます。

PyMCのインストール

Python標準の"venv"を使って、MCMC用の環境を作成し、PyMCをインストールします。

# 最新の状態に更新　
sudo apt update
sudo apt upgrade

# venvのインストール
sudo apt install python3.10-venv

# 仮想環境の構築
python3 -m venv v-MCMC

# v-MCMCに入る →プロンプトの頭に（v-MCMC）と表示される
source v-MCMC/bin/activate

# ライブラリの追加
pip3 install scipy pandas
pip3 install plotly kaleido matplotlib
pip3 install ipykernel
pip3 install  cloudpickle cachetools arviz
pip3 install fastprogress typing-extensions
pip3 install pytensor

# PyMCのインストール
pip3 install pymc
pip3 install git+https://github.com/pymc-devs/pymc-experimental.git

VScodeからv-MCMC環境を呼び出す

↓を参考に、VSCodeにvenv環境を認識してもらいます。

qiita.com

GUIでは、右下の「管理」から「設定」を呼び出し、

"venv"で検索して「項目の追加」をクリック

"v-MCMC"と入力してOK

追加される

「表示」→「コマンドパレット」からJupyterインタラクティブシェルを起動

右上のカーネル名をクリック→「別のカーネルを選択」

「Python環境」を選択

「v-MCMC」を選択

カーネルが変更される

一般化線形モデルを動かしてみる

PyMCの公式には様々なチュートリアルが公開されています。動作確認に“GLM: Linear regression”を実行してみます。

https://www.pymc.io/projects/docs/en/latest/learn/core_notebooks/GLM_linear.html

まずは環境設定と擬似データを作成

# ライブラリの読み込み
import arviz as az
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import pymc as pm
import xarray as xr
from pymc import HalfCauchy, Model, Normal, sample

# 環境設定
RANDOM_SEED = 8927
rng = np.random.default_rng(RANDOM_SEED)

%config InlineBackend.figure_format = 'retina'
az.style.use("arviz-darkgrid")

# 擬似データ作成
size = 200
true_intercept = 1
true_slope = 2

x = np.linspace(0, 1, size)
# y = a + b*x の線形回帰モデル
true_regression_line = true_intercept + true_slope * x
# データにノイズを加える
y = true_regression_line + rng.normal(scale=0.5, size=size)

data = pd.DataFrame(dict(x=x, y=y))

# データのグラフを描いてみる
fig = plt.figure(figsize=(7, 7))
ax = fig.add_subplot(111, xlabel="x", ylabel="y", title="Generated data and underlying model")
ax.plot(x, y, "x", label="sampled data")
ax.plot(x, true_regression_line, label="true regression line", lw=2.0)
plt.legend(loc=0);

次にモデルを組み立てます。チュートリアルでは一気にサンプリングまでしていますが、ここではまず、モデルを組み立てるところまで。

with Model() as model:  # model specifications in PyMC are wrapped in a with-statement
    # Define priors
    sigma = HalfCauchy("sigma", beta=10)
    intercept = Normal("Intercept", 0, sigma=20)
    slope = Normal("slope", 0, sigma=20)

    # Define likelihood
    likelihood = Normal("y", mu=intercept + slope * x, sigma=sigma, observed=y)

チュートリアルにはありませんが、以下のコマンドでモデルの内容を表示できます。

# 定義した各パラメータの情報
display(model.model)

# モデル構造
modeldag = pm.model_to_graphviz(model)
display(modeldag)

サンプリングを実行

with model:
    idata = sample(3000)

結果を表示します。 ※チュートリアルでは"bambi"などなるものを使って、Rで使うようなチルダ式を使ってモデルを作成する方法が紹介されていますが、ここでは割愛します。

# サンプリング結果
az.plot_trace(idata, figsize=(10, 7));

# 予測結果を計算
idata.posterior["y_model"] = idata.posterior["Intercept"] + idata.posterior["slope"] * xr.DataArray(x)

# 予測結果を描く
_, ax = plt.subplots(figsize=(7, 7))
az.plot_lm(idata=idata, y="y", num_samples=100, axes=ax, y_model="y_model")
ax.set_title("Posterior predictive regression lines")
ax.set_xlabel("x");


動きましたか？

次回はこのPyMCを使って、 「日本域CMIP6データ」を解析したいと思います。

「日本域CMIP6データ」を可視化する　その１－データの入手とマップ化

気候予測データについて

はじめに

「日本域CMIP6データ」は統計的ダウンスケーリング手法を用いて作成された、20世紀初頭から21世紀末までの日単位で全国1kmメッシュの気候予測情報で、5種類の最新の全球気候モデル、3種類の温室効果ガス排出想定に基づいた将来予測データが公開されています。

www.nies.go.jp

利用可能な 8 変数（日最低・最高・平均気温、降水量、全天日射量、下向き長波放射、風速、相対湿度）のうち日最高気温データについて、まずはマップ化してみましょう。

  1. 気候予測データセット（DS2022）のホームページにアクセス
    diasjp.net

  2. データセット紹介→日本域CMIP6データに進む

  3. データダウンロードをクリック

  4. "データをダウンロード"をクリック

  5. CMIP6をベースにしたCDFDM手法による日本域バイアス補正気候シナリオデータのダウンロード一覧ページにログイン

  6. ”MRI-ESM2-0”の1950年代の日最高気温データを検索し、全ファイルにチェックを入れる
    ディレクトリ指定：/NIES2020_jpnCDFDM_CMIP6/tasmax/day/MRI-ESM2-0/historical キーワード指定：-195(部分一致)

  7. 一括ダウンロードする

データを解凍する

# ユーザホームディレクトリにCMIP6フォルダを作成し、
# ダウンロードした”data***.tgz”を移動しておく。
cd /mnt/c/Users/hoge/CMIP6/  
find ./ -type f -name "*.tgz" | xargs -n 1 tar zxf

次のディレクトリにデータが解凍されます。 C:\Users\hoge\CMIP6\dias\data\NIES2020_jpnCDFDM_CMIP6\tasmax\day\MRI-ESM2-0

QGISに読み込む

入手したデータはnetcdf形式ですので、そのままQGISに読み込むことができます。

メッシュレイヤとして読み込む

「レイヤ」→「レイヤを追加」→「メッシュレイヤを追加」

ソース：メッシュデータセットに"tasmax_day_MRI-ESM2-0_historical_r1i1p1f1_19500101-19501231_cdfdm.nc"を指定して「追加」

CRSの設定

緯度経度のグリッドですので、EPSG4326をつかいます。

時間設定

レイヤの時間を19500101-19501231になるように設定します。

日最高気温のマップが表示されます 。

一応アニメーションも動きます。

「日本域CMIP6データ」の日最高気温データを入手し、マップ化するまでを紹介しました。

実は、月平均値なら気候変動適応情報プラットフォームのWebGISで簡単に可視化できます。

adaptation-platform.nies.go.jp

簡単ですので、まずはこちら↑を試してみてください。

次回は、入手したデータからピンポイントの情報を取り出し、集計する方法を紹介したいと思います。

※本記事では以下のデータを利用させていただきました。 石崎 紀子, 2021: CMIP6をベースにしたCDFDM手法による日本域バイアス補正気候シナリオデータ, Ver.1, 国立環境研究所, doi:10.17595/20210501.001. (参照: 2023/09/28~10/01)

「全球及び日本域150年連続実験データ」を可視化する　その7－札幌の平均気温のシミュレーション結果をヒートマップで表現する。

気候予測データについて

はじめに

今年の夏は北海道でもうだるような暑さが続きました。気象庁の”夏の天候のまとめ”によれば、2023年（令和5年）夏（6～8月）の日本の平均気温は1898年以降で夏として最も高くなったとのこと。

www.jma.go.jp

別紙（概況、統計値等）によれば、北海道の平均気温は平年より3℃も高かったらしい。

https://www.data.jma.go.jp/obd/stats/data/stat/tenko2023jja_besshi.pdf

実際、これまでに経験したことのない暑さだったことは間違いないのですが、気温の変化をよりわかりすく示すヒートマップというものがあるようです

www.sakigake.jp

www.nishinippon.co.jp

これはわかりやすい！ということで、１５０連続実験データでもこの表現手法を使ってヒートマップを作成してみたいと思います。

データの切り出し

ここでは１例として札幌管区気象台に最も近い格子点における８月の平均気温のヒートマップを作成します。

地理座標は北緯43度3.6分、東経141度19.7分ですので、対応する格子点の位置を特定します。

https://www.jma.go.jp/jma/kishou/know/amedas/ame_master.pdf

# ライブラリをインポート
import os
import numpy as np
import pandas as pd
import geopandas as gpd
from shapely.geometry import Polygon

os.chdir("/mnt/c/Users/hoge/ONEFIFTY")

# 札幌管区気象台の緯度経度
src_lat = 43 + 3.6/60
src_lon = 141 + 19.7/60

# 日本域20㎞150年連続実験データの投影法での座標を計算
wgs_SRS = osr.SpatialReference()
wgs_SRS.ImportFromEPSG(4326)

lcc_SRS = osr.SpatialReference()
lcc_SRS.ImportFromProj4(
"+proj=lcc +lat_1=30.0 +lat_2=60.0 +lon_0=135. +lat_0=0 +x_0=0 +y_0=0 +a=6371000 +b=6371000 +units=m +no_defs")

trans = osr.CoordinateTransformation(wgs_SRS,lcc_SRS)
out_x,out_y,z = trans.TransformPoint(src_lat,src_lon)

#最も近い格子点の行番号、列番号を特定する
x_vec = (-1930000.1287495417 + 10000) + np.arange(0,191) * 20000
y_vec = ((5833346.794147097 - 10000) - np.arange(0,155) * 20000)[::-1]

x_ind = np.argmin(np.abs(x_vec - out_x))
y_ind = np.argmin(np.abs(y_vec - out_y))

# 確認のためのベクタを作成する
x = x_vec[x_ind]
y = y_vec[y_ind]

df = pd.DataFrame(["Sapporo"])
df.columns=["name"]
gm = [Polygon([(x-10000,y-10000),(x-10000,y+10000),(x+10000,y+10000),(x+10000,y-10000),(x-10000,y-10000)])]

sapporo_gdf = gpd.GeoDataFrame(df,geometry=gm,crs=lcc_SRS.ExportToWkt())
sapporo_gdf.to_crs("EPSG:4326").to_file("sappro_grid.fgb", driver="FlatGeobuf")
print(x_ind,y_ind) # 121 121

格子点を中心とした20㎞メッシュに札幌が含まれていることが確認できます。

150年連続実験データの処理

「全球及び日本域150年連続実験データ」の20km150年連続実験データから、（121,121）に位置するデータを切り出し、地上気温の日平均値を計算します。

# ライブラリをインポート
import os
import pathlib
import numpy as np
import pandas as pd
import xarray as xr
import rioxarray

# 処理するファイルの一覧を作成する
os.chdir("/mnt/c/Users/hoge/ONEFIFTY")

griblist_df = pd.DataFrame({"grib_path":list(pathlib.Path("./").glob("**/*.grib"))})
griblist_df["grib_name"] = griblist_df["grib_path"].map(lambda x: x.name)
griblist_df["year"] = griblist_df["grib_name"].str.extract(".*_(\d{4})\d{2}.grib").astype(int)
griblist_df["month"] = griblist_df["grib_name"].str.extract(".*_\d{4}(\d{2}).grib").astype(int)
griblist_df = griblist_df.sort_values("year").reset_index(drop=True)

# 年ごとの日平均値を収納する配列を作成
TA_daily_mean_arr = np.zeros((griblist_df.shape[0],31))

# 年ごとに日平均値を計算
for h_index,h_row in griblist_df.iterrows():
    ds = xr.open_dataset(h_row["grib_path"], engine="cfgrib")
    ds["time"] = ds["time"] - np.timedelta64(1, 'h')
    da = ds["t"][:,121,121].resample(time='1D').mean()
    TA_daily_mean_arr[h_index,:] = da.data

np.save("TA_daily_mean.npy",TA_daily_mean_arr)

ヒートマップの作成

plotlyでヒートマップを描きます。

# ライブラリをインポート
import os
import numpy as np
import plotly.express as px
import matplotlib.cm as cm

os.chdir("/mnt/c/Users/hoge/ONEFIFTY")

#　データをロード
TA_daily_mean_arr = np.load("TA_daily_mean.npy")

#　平年値（統計期間1991～2020年）を計算
TA_normal = TA_daily_mean_arr[40:70,:].mean()

# ヒートマップの描画
fig = px.imshow(TA_daily_mean_arr, aspect="auto",width=550,height=800, 
color_continuous_scale='RdBu_r',
color_continuous_midpoint=TA_normal)

# レイアウト編集
fig.update_layout(
    title=dict(text="札幌周辺の８月の日平均気温シミュレーション（1951~2099）",
                y=0.98),
    margin=dict(t=40, b=45, l=25, r=120)
)   

#脚注１を作成
fig.add_annotation(
    font_size=10,
    font_family="sans-serif",
    x=1.04,
    y=0.83,
    xref="paper",
    yref="paper",
    xanchor="left",
    yanchor="top",
    borderwidth=1,
    bordercolor='black',
    align="left",
    text="「日本域150年連続<br>実験」による、札幌<br>に最も近い格子点に<br>
　　　　おける1951年から<br>2099年までの８月の<br>日平均気温シミュレ<br>
　　　　ーション結果。縦軸<br>が '年' , 横軸が '月日' 。<br>8月の日平均気温の<br>
　　　　平年値（1991~2020<br>年）を白とし、高温<br>を赤、低温を青とし<br>
　　　　た色分けにより気温<br>の違いを表現した。<br>結果は特定の日の「<br>
　　　　気象」の状態（何月<br>何日に気温は何度か）<br>を予測しているもの<br>
　　　　ではないことに注意！"
)

# 脚注２を作成
fig.add_annotation(
    font_size=10,
    font_family="sans-serif",
    x=1.03,
    y=0.12,
    xref="paper",
    yref="paper",
    xanchor="left",
    yanchor="top",
    align="left",
    showarrow=False,
    text="本ヒートマップは、<br>文部科学省「統合<br>的気候モデル高度<br>
　　　　化研究プログラム」<br>において、地球シ<br>ミュレータを用い<br>
　　　　て作成されたデー<br>タを使用させてい<br>ただいています。"
)

# y軸を設定
fig.update_yaxes(
    tickvals=[0,9,19,29,39,49,59,69,79,89,99,109,119,129,139,148],
    ticktext=["1951","1960","1970","1980","1990","2000","2010","2020","2030","2040","2050","2060","2070","2080","2090","2099"]
    )

# x軸を設定
fig.update_xaxes(
    tickvals=[0,4,9,14,19,24,29],
    ticktext=["08/01","08/05","08/10","08/15","08/20","08/25","08/30"],
    tickangle=-50
)

# カラーバーを設定
fig.update_coloraxes(
    colorbar_len=0.4,
    colorbar_y=0.32,
    colorbar_yanchor="middle",
    colorbar_ticklabeloverflow="allow",
    colorbar_tickvals=[284,310],
    colorbar_ticktext=["low","high"],
    colorbar_title="日平均気温",
    colorbar_title_side="top"
    )

# ロゴを描画
fig.add_layout_image(
    sizex=0.2,sizey=0.2,
    source="https://pbs.twimg.com/profile_images/1664190624701222912/8gaaJj-l_400x400.png",
    xref="paper",
    yref="paper",
    x=1.05,
    y=0.995)

fig.show()

ヒートマップができました。

皆さんのお住まいの地域でも同様のヒートマップを作成することができます。

「全球及び日本域150年連続実験データ」の日本域20km150年連続実験データの日平均気温をヒートマップで表現する方法を紹介しました。試してみてくださいね～。

※本記事では文部科学省「統合的気候モデル高度化研究プログラム」において、地球シミュレータを用いて作成されたデータを使用しました。

「全球及び日本域確率的気候予測データ（d4PDF シリーズ）」を可視化する － 日本域ダウンスケーリングデータを用いて大雨の将来予測を可視化する。

気候予測データについて

はじめに

「全球及び日本域確率的気候予測データ（d4PDF シリーズ）」は、「大規模アンサンブル気候シミュレーション出力をまとめたデータベースの総称」で、「社会的関心の高い豪雨や台風などの発現頻度が低い顕著な大気現象の将来変化を評価・推定することに重点を置いて設計されたデータベース」です。（気候予測データセット解説書第１章．P1-9より抜粋）

大雨といえば、北海道でも8月はじめにこんなニュースがありました。

www3.nhk.or.jp

将来このような気象が、どのくらいの頻度で発生すると予測されるか？d4PDFデータを用いて可視化してみましょう。

データ入手

20kmの地域気候モデルのデータのうち、北海道が含まれる範囲を切り出してダウンロードします。

  1. 気候予測データセット(DS2022)のホームページにアクセス
    diasjp.net

  2. データセット紹介→全球及び日本域確率的気候予測データ（d4PDF シリーズ）に進む

  3. 日本域ダウンスケーリングをクリック

  4. データ取得/閲覧をクリック

  5. d4PDF_RCM(Subset)タブを選択

  6. 切り出し条件を設定（北海道エリア）

    • 実験：過去実験：HPB
    • 期間：1951 08 - 1951 08
    • 変数：地上大気データ：surf → rain：Precipitation
    • アンサンブル：Select All
    • 領域：Xleft 112, Xright 140, Ybottom 110, Ytop 135


  7. 切り出しをクリック

"subset_tar"というファイルがダウンロードされます。

データ切り出しスクリプトの利用

HPBだけでも1951年から2011年まで60年分マウスクリックするのはつらいので、ダウンロードスクリプトを使用します。

スクリプトはページの最後にありますので、ダウンロードしてホームディレクトリに作成した"d4PDF"フォルダに保存し、WSL2コンソールで解凍します。

cd　/mnt/c/Users/hoge/d4PDF
tar -xvf d4pdf-subset-script-3.tgz

先ほどマウスで行った切り出し処理のコマンドはこのようになります。

python3 d4pdf-sb.py -o subset.tar -f 1951-08 -t 1951-08 -X 112,140 -Y 110,135 RCM/HPB/Surf rain

実行時にユーザのホームディレクトリにある.netrcを参照するようなので、作成しておきます。

nano ~/.netrc

DIASのユーザー設定を記入します。

machine d4pdf.diasjp.net
login {DIASのユーザー名}
password {DIASのパスワード}

ファイルのパーミッションを変更しておきます。

chmod og-rw ~/.netrc

Pythonのsubprocessを用いて1年づつ実行します。

# ライブラリのインポート
import os
import numpy as np
import pandas as pd
import subprocess
os.chdir("/mnt/c/Users/hoge/d4PDF")

# HPB切り出し条件の設定
subset_df = pd.DataFrame({"year":np.arange(1951,2012)})
subset_df["output"] = subset_df["year"].map(lambda x: "HPB_{}.tar".format(x))
subset_df["from"] = subset_df["year"].map(lambda x: "{}-08".format(x))
subset_df["to"] = subset_df["from"]

for _,hrow in subset_df.iterrows():
    cmd = ["python3","d4pdf-sb.py","-o",hrow["output"],
　　　　"-f",hrow["from"],"-t",hrow["to"],"-X","112,140","-Y","110,135",
　　　　"RCM/HPB/surf","rain"]
    subprocess.run(cmd)

# HFB_4K_CC切り出し条件の設定
subset_df = pd.DataFrame({"year":np.arange(2051,2111)})
subset_df["output"] = subset_df["year"].map(lambda x: "HFB_4K_CC_{}.tar".format(x))
subset_df["from"] = subset_df["year"].map(lambda x: "{}-08".format(x))
subset_df["to"] = subset_df["from"]

# 切り出し実行
for _,hrow in subset_df.iterrows():
    cmd = ["python3","d4pdf-sb.py","-o",hrow["output"],
　　　　"-f",hrow["from"],"-t",hrow["to"],"-X","112,140","-Y","110,135",
　　　　"RCM/HFB_4K_CC/surf","rain"]
    subprocess.run(cmd)

# 'HFB_4K_GF' 'HFB_4K_HA' 'HFB_4K_MI' 'HFB_4K_MP' 'HFB_4K_MR' も同様にダウンロードします。
tarの解凍

WSL2コンソールから解凍します

cd　/mnt/c/Users/hoge/d4PDF
find ./ -type f -name "*.tar" | xargs -n 1 tar xf

「平年8月1か月分の総雨量」を求める

過去実験（HPB）のすべての年、すべてのアンサンブルについて8月の総雨量を計算し、その平均値を「平年8月1か月分の総雨量」とします。

ここでは、時間雨量0.5㎜未満を無降水として扱いました。

# ライブラリのインポート
import os
import pathlib
import numpy as np
import pandas as pd

os.chdir("/mnt/c/Users/hoge/d4PDF")

# 処理ファイルのリストを作成する
files_df = pd.DataFrame({"f_path":list(pathlib.Path("/mnt/c/Users/hoge/d4PDF/d4PDF_RCM/HPB").glob("**/*rain_subset.dat"))})

# 月総雨量を収納する配列を作成しておく
monthly_arr = np.zeros((files_df.shape[0],26,29))

# 1ファイルずつ処理
for ind,hrow in files_df.iterrows():

    # バイナリファイルを読み込み
    dvec = np.fromfile(hrow["f_path"],dtype='>f4')

    # 4次元の配列に変換
    darr = dvec.reshape(744,26,29)

    # 総雨量を計算して、空の配列に格納
    monthly_arr[ind,:,:] = darr.sum(axis=0)

# 平均年値を計算して保存しておく
climate_normal = monthly_arr.mean(axis=0)
np.save("./aug_rain/climate_normal_sumrain.npy",climate_normal)

結果を図にしてみると、こんな感じです。

メッシュ平年値図と比べて、値の範囲的には一致しているようです。

https://www.data.jma.go.jp/obd/stats/etrn/view/atlas/docs2020/png_large/prec/precipitation_08.png

https://www.data.jma.go.jp/obd/stats/etrn/view/atlas/docs2020/png_large/prec/precipitation_08.png

降雨イベントの検出と統計処理

雨の降り始めから降り終わりまでを１つの降雨イベント(ひと雨)とし、降り始めからの積算雨量が８月の総雨量の平年値を超える回数を計算します。ここでは、24時間以上の無降水が継続したときを降雨イベントの区切りとしました。

# ライブラリのインポート
import os
import pathlib
import numpy as np
import pandas as pd

os.chdir("/mnt/c/Users/hoge/d4PDF")

def res_counts_total_storm_prcip(ts_vec,threshold):
    #24hr雨量の計算
    mvsum_24hr_vec = np.concatenate(
        [np.zeros(23),
        np.convolve(ts_vec,np.ones(24),mode="valid")]
    )

    #降雨イベントのインデックスを作成
    pr_event_no =np.concatenate(
        [np.zeros(1),
        (np.convolve(mvsum_24hr_vec == 0,(-1,1),mode="valid") == 1).cumsum()]
    ).astype(int)

    #イベントごとの積算雨量を計算
    TP_vec = np.bincount(pr_event_no, weights=ts_vec)

    #閾値を超える回数を集計して返す
    return((TP_vec > threshold).sum())

# 8月の総雨量の平年値を読み込んでおく
climate_normal = np.load("./aug_rain/climate_normal_sumrain.npy")

# ファイルリストの作成
files_df = pd.DataFrame({"f_path":list(pathlib.Path("./d4PDF_RCM/").glob("**/*rain_subset.dat"))})
files_df["f_name"] = files_df["f_path"].map(lambda x:x.name)
files_df["ens_name"] = files_df["f_name"].str.extract("surf_(.*)_\d{6}_rain.*")

# アンサンブルごとに超過回数を計算する
for e_name,hgrp in files_df.groupby("ens_name"):
    # 返り値の入れ物を作成
    counts_arr = np.zeros((hgrp.shape[0],26,29))

    # インデックスをリセット
    hgrp = hgrp.reset_index(drop=True)

    for Tind,hrow in hgrp.iterrows():
        # ファイルの読み込み
        dvec = np.fromfile(hrow["f_path"],dtype=">f4")
        dvec = dvec.astype(float)

        # 4次元の配列に変換
        darr = dvec.reshape(744,26,29)

        # 時間雨量0.5㎜を無降水に
        darr[darr<0.5] = 0

        # 次元を入れ替えて、Y,X,timeにする
        darr = darr.transpose(1,2,0) 

        # 格子点ごとに処理
        for Yind,X_darr in enumerate(darr):
            for Xind,ts_vec in enumerate(X_darr):
                thresh = climate_normal[Yind,Xind]
                counts_arr[Tind,Yind,Xind] = res_counts_total_storm_prcip(ts_vec,thresh)

    # npy形式で保存
  np.save("./aug_rain/counts_over_threshold_rain_total_storm_prec_{}.npy".format(e_name),counts_arr)

大雨の発生確率の計算

「降り始めからの積算雨量が閾値を超えるような降雨イベント」が発生する回数はポアソン分布に従うと仮定し、年に1回以上発生する確率を計算します。

# ライブラリ呼び出し
import os
import pathlib
import numpy as np
import pandas as pd
from osgeo import gdal, gdalconst, gdal_array 
from osgeo import osr
from scipy.stats import poisson

# 地理座標系のパラメータを設定
ul_x = -1930000.1287495417 + 112*20000
ul_y = 5833346.794147097 - (155-135-1)*20000
x_res = 20000
y_res = -20000
x_size = 29
y_size = 26

outRasterSRS = osr.SpatialReference()
outRasterSRS.ImportFromProj4(
"+proj=lcc +lat_1=30.0 +lat_2=60.0 +lon_0=135. +lat_0=0 +x_0=0 +y_0=0 +a=6371000 +b=6371000 +units=m +no_defs")

# ファイルリストの作成
files_df = pd.DataFrame({"f_path":list(pathlib.Path("./aug_rain").glob("counts*.npy"))})
files_df["f_name"] = files_df["f_path"].map(lambda x: x.name)
files_df["exper"] = files_df["f_name"].str.extract(".*_prec_([A-Z]*)_.*")

for exper,hgrp in files_df.groupby("exper"):
    # 発生回数の平均を計算
    sum_arr = np.zeros((26,29))
    sum_counts = 0

    for _,hrow in hgrp.iterrows():
        counts_arr = np.load(hrow["f_path"])
        counts_arr[counts_arr > 1] = 1
        sum_arr += counts_arr.sum(axis=0)
        sum_counts += counts_arr.shape[0]

    lambda_arr = sum_arr/sum_counts

    # 1回以上発生する確率
    PO_arr = np.vectorize(lambda x: 1-poisson.cdf(0,x))(lambda_arr)
    
    # geotiffに保存
    name = "Probability_occurrence_heavy_rain{}.gtiff".format(exper)
    output = gdal.GetDriverByName('GTiff').Create(name,x_size,y_size,1,gdal.GDT_Float64)
    output.SetProjection(outRasterSRS.ExportToWkt())
    output.SetGeoTransform((ul_x, x_res, 0, ul_y, 0, y_res))

    outband = output.GetRasterBand(1)
    outband.WriteArray(PO_arr[::-1,:].astype("float64"))
    outband.FlushCache()
    output = None

QGISに読み込む

「レイヤ」→「レイヤを追加」→「ラスタレイヤを追加」

ソース：ラスタデータセットに”Probability_occurrence_heavy_rainHPB.gtiff”を指定して「追加」

過去実験：HPBにおける「降り始めからの積算雨量が閾値を超えるような降雨イベント」が発生する確率のマップが表示されます。

将来４度昇温実験：HFBのマップ

変化率の計算

QGISのラスタ計算機を使えば、過去実験と将来４度昇温実験における大雨発生確率の変化倍率を計算できます。

「ラスタ」→「ラスタ計算機」

計算条件を設定

  • ディスクに書き込まないオンザフライ・ラスタを作成にチェック
  • レイヤ名：HPB_HFB_ratio
  • 式：”Probability_occurence_heavy_rainHFB@1"/”Probability_occurence_heavy_rainHPB@1"

変化倍率のマップが計算されます。

と、この記事を書いていたらこんなニュースが。

mainichi.jp

気象庁気象研究所のプレスリリースを読むと、

www.mri-jma.go.jp

  • d4PDFは領域版でもメッシュ間隔が20kmであるため、狭い範囲に短時間に多量の雨をもたらす線状降水帯や、複雑な地形の影響を受けた大雨を再現することはできませんでした。　
  • 今回実施した5kmメッシュ（水平格子間隔）のシミュレーションは、d4PDFの20kmメッシュのシミュレーションに比べて、発生頻度の低い大雨の再現性が向上し線状降水帯の検出も可能とな りました。
  • 本データセットはデータ統合・解析システム（DIAS）を通じて公開する予定です。

とのことですので、公開を待ちたいと思います。

「全球及び日本域確率的気候予測データ（d4PDF シリーズ）」の 日本域ダウンスケーリングデータを用いて大雨の将来予測を可視化する方法を紹介しました。試してみてくださいね～。

※本記事では、文部科学省による複数の学術研究プログラム（「創生」、「統合」、SI-CAT、DIAS）間連携および地球シミュレータにより作成された d4PDF を使用させていただきました。

「全球及び日本域150年連続実験データ」を可視化する その6－回る地球上に地上気温変化を描く。

気候予測データについて

はじめに

前回は、QGIS上で全球の地上気温変化アニメーションを作成しました。

cci-labo.hateblo.jp

ただやはり、本物の地球を観察しているような視点で描きたい！ ということで（画面上で）回る地球儀にデータを貼り付けてアニメーションを作ってみます。

温度変化マップの作成

前回作成した年平均気温の変化データからcartpyを用いて正距円筒図法のマップを作成します。

# ライブラリのインポート
import os
import pandas as pd
import xarray as xr
import matplotlib.pyplot as plt
import cartopy.crs as ccrs
from matplotlib.colors import LinearSegmentedColormap, Normalize
from matplotlib.colorbar import ColorbarBase

# カラーマップの作成
colors = ["blue", "#ffffe1", "red", "black"]
nodes = [0.0, 0.2, 0.6, 1.0]
TC_cmap =  LinearSegmentedColormap.from_list("tccmap", list(zip(nodes, colors)))

# データの読み込み
os.chdir("/mnt/c/Users/hoge/ONEFIFTY")
ds = xr.open_dataset("GCM60_RCP85_tmpchange_sfc_avr_year.nc")
relief = plt.imread("../natural_earth/MSR_50M/MSR_50M.tif")

# 正距円筒図法のマップを作成
for h_date,h_arr in ds.groupby("time"):
    fig = plt.figure(figsize=(24, 12))
    ax = fig.add_subplot(1,1,1,projection=ccrs.PlateCarree())
    h_arr["temperature"].plot.imshow(ax=ax, cmap=TC_cmap, vmin=-5.0,vmax=20.0,
                                     transform=ccrs.PlateCarree(),add_colorbar=False,add_labels=False)
    ax.imshow(relief,alpha=0.2, cmap="gray",extent=(-180,180,-90,90))
    ax.coastlines(resolution='50m',linewidth=.1)
    ax.gridlines(crs=ccrs.PlateCarree(), draw_labels=False, linewidth=1, alpha=0.8)

    plt.axis('tight')
    plt.axis('off')
    fig.subplots_adjust(left=0, right=1, bottom=0, top=1)

    #保存
    plt.savefig(pd.to_datetime(h_date).strftime("./TCanu_latlon_png/anu_%Y.png"))
    plt.clf()
    plt.close()

# カラーバーを描く
fig = plt.figure(figsize=(1,7))
ax = fig.add_axes([0.05, 0.05, 0.9, 0.9])
norm = Normalize(vmin=-5.0, vmax=20.0)
cb = ColorbarBase(ax, orientation='vertical', 
                           cmap=TC_cmap,norm=norm, extend='both')
cb.ax.yaxis.set_tick_params(color="white")
cb.ax.set_yticklabels(["Cold","0","","","","Warm"])
cb.outline.set_edgecolor("white")
plt.setp(plt.getp(cb.ax.axes, 'yticklabels'), color="white",fontsize=20)

fig.patch.set_facecolor("black")   

plt.savefig('TC_cb_var.png', bbox_inches='tight', facecolor="black")

アニメーション作成

アニメーションロゴの作成ではGIMP上でGIFを作成しましたが、GIFを使うと256階調のインデックスカラーに変換されてしまい、グラデーションがうまく出ない・・・

qiita.com

そこで1フレームずつpngを作成し、連番画像からffmpegを使ってmp4を作ることにしました。

タイトル−キャプション画像の作成

GIMPを用いて、アニメーション動画の背景とタイトル・キャプションを作成し、pngで保存しておきます。

フレーム画像の作成

作成したマップをGIMPを使って地球儀に変換し、タイトル−キャプション画像に貼り付けてからpngに書き出します。ファイル名等の文字列の扱いはPythonでやったほうが楽なので、一連の処理をするスプリクトを作成し、Pythonで実行します。

処理スクリプト

1枚のマップを読み込んで地球儀に貼り付け、回転角を調節してから背景画像に貼り付け、pngに書き出すScript-fu

(define (latlonmap-to-flame base_png inFile outDir yr_str rotate1 rotate2 sq-width)
    (let* 
        (
            ;変数の設定
            (base_img (car (gimp-file-load RUN-NONINTERACTIVE  base_png  base_png)))
            (globe_img (car (gimp-file-load RUN-NONINTERACTIVE inFile inFile)))
            (src-layer (car (gimp-image-get-active-layer globe_img)))
            (globe_layers (cadr (gimp-image-get-layers globe_img)))
            (rotate_y rotate1)
            (new_globelayer 0)
            (l 0)
            (l_name "fl_")
            (layer_id 0)           
        ) 

        ;レイヤー名（fl-{year}_{number}）を返す関数
        (define (res_layer_name yr_str layer_id)
            (let* (
                (s (string-append "0000" (number->string layer_id)))
                (n (string-length s))
                )
                (string-append "fl-" yr_str "_" 
                    (substring s (- n 2) n)
                )
            )
        )

        ;球体に変換したレイヤを背景に結合し、pngで書き出す関数
        (define (create_flame l)
            (let* (
                (new_flame (car (gimp-image-duplicate base_img )))
                (new_flamelayer (car (gimp-layer-new-from-drawable l new_flame)))
                (l_name (car (gimp-layer-get-name l)))
                (outFile (string-append outDir "/" l_name ".png"))
                (tx_layer 0) 
                (d_able 0)
                )
                (gimp-image-add-layer new_flame new_flamelayer 0)
                (gimp-layer-set-offsets new_flamelayer 50 0)
                (set! tx_layer (car (gimp-text-fontname new_flame new_flamelayer 10 0 
　　　　　　　　　　　　　yr_str -1 TRUE 100 POINTS "Rounded-L M+ 2p Heavy")))
                (gimp-floating-sel-to-layer tx_layer)
                (gimp-text-layer-set-color tx_layer '(255 255 255))
                (gimp-image-merge-visible-layers new_flame CLIP-TO-BOTTOM-LAYER )
                (set! d_able (car (gimp-image-get-active-drawable new_flame) ))
                (file-png-save-defaults RUN-NONINTERACTIVE new_flame d_able outFile outFile)
                (gimp-image-delete new_flame)
            )
        )

        ;マップを縦横同じサイズの画像に変換
        (gimp-image-scale globe_img sq-width sq-width)

       ;球体に変換し0.5刻みで回転させたレイヤを作成
        (while (< rotate_y rotate2) 
            (set! new_globelayer (car(gimp-layer-copy src-layer TRUE)))
            (set! l_name (res_layer_name yr_str layer_id))
            (gimp-layer-set-name new_globelayer l_name)
            (gimp-image-add-layer globe_img new_globelayer 0)
            (plug-in-map-object 1 globe_img new_globelayer 1 0.5 0.5 20.0 0.5 0.5 0.0 
　　　　　　　　　　　　0 0 0 0 0 0 0 rotate_y 35 1 '(255 255 255) 0.5 0.5 10 0 0 1 
　　　　　　　　　　　　0.4 0.9 0.4 0 27.0 TRUE FALSE FALSE TRUE 0.48 0.5 0.5 0.5 
　　　　　　　　　　　　1.0 -1 -1 -1 -1 -1 -1 -1 -1)
            (set! rotate_y (+ rotate_y 0.5))
            (set! layer_id (+ layer_id 1))
        )
        (gimp-image-remove-layer globe_img src-layer)

        ;レイヤごとに背景と結合し、保存する。
        (set! globe_layers (cadr (gimp-image-get-layers globe_img)))
        (for-each 
            (lambda (l)
                (create_flame l)
            )
            (vector->list globe_layers)
        )

        (gimp-image-delete base_img)
        (gimp-image-delete globe_img)

    )
)

;GIMPへの登録設定
(script-fu-register
 "latlonmap-to-flame"                        ;func name
 "latlonmap-to-flame"                                  ;menu label
 "月平均値アニメーションのフレームを作成" ;description
 "ccilobo"                             ;author
 "copyright 2023, ccilabo"        ;copyright notice
 "Sep 14, 2023"                          ;date created
 ""                     ;image type that the script works on
 )

(script-fu-menu-register "latlonmap-to-flame" "<Image>/Script-Fu/CCILabo")
1年づつ処理する

pythonのsubprocessを用いて 1ファイルづつ処理し、最後にファイル名に連番を振ります。

# ライブラリの読み込み
import os
import pathlib
import numpy as np
import pandas as pd
import subprocess

# 処理ファイルのリストを作成
maps_df = pd.DataFrame({"f_path":list(pathlib.Path("/mnt/c/Users/hoge/ONEFIFTY/TCanu_latlon_png").glob("*.png"))})
maps_df["f_name"] = maps_df["f_path"].map(lambda x: x.name)
maps_df["yr_str"] = maps_df["f_name"].str.extract("anu_(\d{4}).png")

# 30年で1回転（1年あたり24フレーム）するように開始角：終了角を設定する
maps_df["rotate1"] = maps_df.index.values*24%720/2 - 180
maps_df["rotate2"] = maps_df.index.values*24%720/2 - 168

# 1ファイルずつ処理
for ind,hrow in maps_df.iterrows():
    cmd_str = 'gimp-2.10 -d -i -b \'(latlonmap-to-flame \"/mnt/c/Users/hoge/ONEFIFTY/annually_animebase.png\" \"{0}\" \"/mnt/c/Users/hoge/ONEFIFTY/TCanu_flame_png\" \"{1}\" {2} {3} 918)\' -b \'(gimp-quit 0)\''.format(str(hrow["f_path"]),hrow["yr_str"],hrow["rotate1"],hrow["rotate2"])
    subprocess.run(cmd_str,shell=True)

#　ファイル名を連番に変更
flame_df = pd.DataFrame({"f_path":list(pathlib.Path("/mnt/c/Users/hoge/ONEFIFTY/TCanu_flame_png").glob("*.png"))})
flame_df["f_stem"] = flame_df["f_path"].map(lambda x: x.stem)
flame_df = flame_df.sort_values("f_stem").reset_index(drop=True)
flame_df["f_serial"] =flame_df.apply(lambda x: "fl-{:04}.png".format(x.name),axis=1)

for ind,hrow in flame_df.iterrows():
    hrow["f_path"].rename(hrow["f_path"].parent / hrow["f_serial"])
動画に変換

ffmpegを使ってmp4動画に変換します

qiita.com

cd /mnt/d/Users/kazuh/ONEFIFTY/
ffmpeg -framerate 30 -i ./TCanu_flame_png/fl-%04d.png -vcodec libx264 -pix_fmt yuv420p -r 60 TC_globe_annualy_150yr.mp4

mp4動画ができました。



youtu.be

「全球及び日本域150年連続実験データ」の全球60km150年連続実験データの温度変化を”回る地球儀”上でアニメーションする方法を紹介しました。試してみてくださいね～。

※本記事では文部科学省「統合的気候モデル高度化研究プログラム」において、地球シミュレータを用いて作成されたデータを使用しました。

「全球及び日本域150年連続実験データ」を可視化する その5ー統合プログラム全球60km150年連続実験データから温度変化のアニメーションを作成・表示する。

気候予測データについて

はじめに

前回は60km格子の全球気候情報を使って世界的な温暖化傾予測データをアニメーション化しました。

ただ、全球で同じ気温スケールを使用しているため、高緯度域と低緯度域で色の変化域が異なり、気温がどう上がっていくかを直感的に把握するのはちょっと難しい。

同様のアニメーションを検索すると、NASAの公開する1880~2022年の全球温暖化アニメーションが引っかかります。

youtu.be

このアニメーションでは、1951-1980年を基準（ゼロ）とし、そこからどの程度、温度が変化するかを描いています。

これを真似させてもらって、150年連続実験データの温度変化のアニメーションを作成してみましょう。

温度変化データの作成

全球60km150年連続実験データの1951-1980年の平均気温マップを作成し、一律に引くことで温度変化を計算します。

# ライブラリのインポート
import os
import pathlib
import numpy as np
import pandas as pd
import xarray as xr
import affine
import rioxarray

# 年平均データの読み込み
os.chdir("/mnt/c/Users/hoge/ONEFIFTY")
year_mean_ds = xr.open_dataset("GCM60_RCP85_tmp_sfc_avr_year.nc")

# 1951~1980の平均を計算
mean_da = year_mean_ds["temperature"].sel(time=slice('1951','1980')).mean(dim="time")

# 温度変化を計算
temprature_change_ds = year_mean_ds - mean_da

このままマップ化するとグリニッジを中心とした世界地図になり、一番描きたい日本周辺が東端に追いやられてしまうので、東経135度を中心としたデータに作り直してからnetCDF形式で出力します。

# 温度変化データの配列を取り出す
tc_arr = temprature_change_ds["temperature"].data

# 西経45度で切り、左右を入れ替える
JP_tc_arr = np.concatenate([tc_arr[:,:,240:],tc_arr[:,:,:240]],axis=2)

# 格子の緯度情報、経度情報を作成
FLON_vec = (np.arange(640) -79.5) * 0.56250 
FLAT_vec = (np.arange(320) -159.5 ) * 0.5616063492063492

# Xarrayのデータセットを作成
JP_tc_ds = xr.Dataset(
    {
        "temperature": (["time", "y","x" ],JP_tc_arr)
    },
    coords={
        "y": ("y", FLAT_vec),
        "x": ("x", FLON_vec),
        "time": temprature_change_ds["time"].data,
        "reference_time": pd.Timestamp("1950-01-01"),
    },
)

# 投影情報を設定
JP_tc_ds = JP_tc_ds.rio.write_crs("+proj=longlat +over +a=6371000 +b=6371000 +no_defs")
JP_tc_ds = JP_tc_ds.rio.write_transform()

# netCDFへ保存
JP_tc_ds.to_netcdf("GCM60_RCP85_tmpchange_sfc_avr_year.nc")

QGISでアニメーション

QGISのメッシュレイヤに読み込み、アニメーションしてみます。

QGISに読み込む

「レイヤ」→「レイヤを追加」→「メッシュレイヤを追加」

ソース：メッシュデータセットに"GGCM60_RCP85_tmpchange_sfc_avr_year.nc"を指定して「追加」


カラーランプを編集

レイヤパネルの「GGCM60_RCP85_tmpchange_sfc_avr_year.nc」を右クリックし、レイヤプロパティを立ち上げます。

「シンボロジ」タブ→「等高線」タブを選択し、以下の項目を設定し、適用します。

  • 最小： -5.0
  • 最大： 20.0
  • リサンプリング方法：隣接平均
  • 補間方法：線形（Linear）色が連続的に変化

カラーランプのプルダウンから「カラーランプを新規作成」を選択します。

「勾配グリッド（Gradient）」を選択してOK

色１に青（RGB:0,0,255）、色２に黒（RGB:0,0,0）を設定。

カラーバーをダブルクリックしてグラディエーションストップを追加し、相対位置と色を設定。

  • 相対位置：20％　色：淡黄（RGB:255,255:225）
  • 相対位置：60％　色：赤(RGB:255,0,0)


投影法の変更

NASAのアニメーションはロビンソン図法で描かれているので、そのように変更します。

投影法の登録

「設定」→「カスタム投影法」で”オプション－ユーザ定義CRS”を開き、以下の投影法を追加します。

  • 名前：JP_robin
  • 形式：Proj文字列（非推奨）
  • パラメータ：+proj=robin +lon_0=135 +x_0=0 +y_0=0 +datum=WGS84 +units=m +no_defs
プロジェクトの投影法の変更

右下の「現在のCRS」をクリック

「あらかじめ定義されたCRS」エリアの「ユーザー定義の座標系」から「JP_robin」を選択して「OK」

ロビンソン図法で表示されます。

アニメーション表示

時系列コントローラを利用してアニメーションできます。

年によって暑くなったり、寒くなったり、変動しながらだんだん世界的に気温が上がっていく様子がアニメーション化できました。

他の気象要素についても同様に作成することができます。

「全球及び日本域150年連続実験データ」の全球60km150年連続実験データの温度変化をアニメーションする方法を紹介しました。試してみてくださいね～。

※本記事では文部科学省「統合的気候モデル高度化研究プログラム」において、地球シミュレータを用いて作成されたデータを使用しました。