「全球及び日本域150年連続実験データ」を可視化する その4―統合プログラム全球60km150年連続実験データを入手しQGIS上でアニメーション表示する
はじめに
「全球及び日本域150年連続実験データ」は1950 年から 2099 年までの連続した 150 年のシミュレーション実験の出力データで、60km格子の全球気候情報と20km格子の日本域気候情報が公開されています。
60km格子の全球気候情報を用いることにより、世界的な温暖化傾向を可視化することが可能です。
日本域20㎞と同様に、QGIS上でアニメーション表示してみましょう。
netCDFの出力にはXarrayを使います。インストールしておいてくだいさい。
データ入手
ここでは、RCP8.5シナリオの月平均値データを入手します。
統合プログラム全球60km領域20km150年連続実験データセットのダウンロード一覧ページにログイン
過去実験、1950年代、地上月平均データを検索する。
全てのファイルにチェックを入れ一括ダウンロードする →2010年代まで繰り返す。
rcp8.5、2010年代、地上月平均データを検索する。
全てのファイルにチェックを入れ一括ダウンロードする →2090年代まで繰り返す。
ctlファイルを検索する
HPDとHFD_HighResMIPのsfc_avr_mon.ctlにチェックを入れ、一括ダウンロードする.
データを解凍する
# ユーザホームディレクトリのONEFIFTYフォルダに、 # ダウンロードした”data***.tgz”を移動しておく。 cd /mnt/c/Users/hoge/ONEFIFTY/ find ./ -type f -name "*.tgz" | xargs -n 1 tar zxf
次のディレクトリにデータが解凍されます。
C:\Users\hoge\ONEFIFTY\dias\data\GCM60_NHRCM20_150yr_TOUGOU\GCM60
GrADSで表示
GrADSでデータを読み込む際に必要な CTL ファイルをダウンロードしたので、まずはGrADSで表示してみます。
#データフォルダに移動 cd /mnt/d/Users/kazuh/ONEFIFTY/dias/data/GCM60_NHRCM20_150yr_TOUGOU/GCM60/HPD #GrDASを立ち上げ grads
#データの読み込み open sfc_avr_mon.ctl #コントロールファイルの内容を表示 q file 1 #timestepを指定 set t 3 #グリッドを塗り潰して色別で表示。 set gxout grfill #メルカトル図法 set mproj latlon #地表気温マップを表示 d ta
GrADSでアニメーションしてみる
set mproj robinson set t 1 780 set loopdim t d ta
実際はアニメーションするハズ
netCDFファイルの作成
GrADS Scriptを使えばもっと手の込んだアニメーションを作れるのですが、QGISで扱うことにより直感的な操作が可能になります。
そこで、地上気温の時系列データを切り出し、QGISの理解できるnetCDFファイルを作成します。
コントロールファイルから情報を抽出
sfc_avr_mon.ctlの内容は
DSET ^%y4%m2/sfc_avr_mon_HPD_%y4%m2.dr
OPTIONS TEMPLATE BIG_ENDIAN TEMPLATE
TITLE sfc_avr_mon_HPD_195001 AVR ( Time average data )
UNDEF -9.99E33
XDEF 640 LINEAR 0.00000 0.56250
YDEF 320 LEVELS
-89.57009 -89.01318 -88.45297 -87.89203 -87.33080
-86.76944 -86.20800 -85.64651 -85.08499 -84.52345
-83.96190 -83.40033 -82.83876 -82.27718 -81.71559
-81.15400 -80.59240 -80.03080 -79.46920 -78.90760
-78.34600 -77.78439 -77.22278 -76.66117 -76.09956
-75.53795 -74.97634 -74.41473 -73.85311 -73.29150
-72.72988 -72.16827 -71.60665 -71.04504 -70.48342
-69.92181 -69.36019 -68.79857 -68.23695 -67.67534
-67.11372 -66.55210 -65.99048 -65.42886 -64.86725
-64.30563 -63.74401 -63.18239 -62.62077 -62.05915
-61.49753 -60.93591 -60.37429 -59.81267 -59.25105
-58.68943 -58.12781 -57.56619 -57.00457 -56.44295
-55.88133 -55.31971 -54.75809 -54.19647 -53.63485
-53.07323 -52.51161 -51.94999 -51.38837 -50.82675
-50.26513 -49.70351 -49.14189 -48.58026 -48.01864
-47.45702 -46.89540 -46.33378 -45.77216 -45.21054
-44.64892 -44.08730 -43.52567 -42.96405 -42.40243
-41.84081 -41.27919 -40.71757 -40.15595 -39.59433
-39.03270 -38.47108 -37.90946 -37.34784 -36.78622
-36.22460 -35.66298 -35.10136 -34.53973 -33.97811
-33.41649 -32.85487 -32.29325 -31.73163 -31.17000
-30.60838 -30.04676 -29.48514 -28.92352 -28.36190
-27.80028 -27.23865 -26.67703 -26.11541 -25.55379
-24.99217 -24.43055 -23.86892 -23.30730 -22.74568
-22.18406 -21.62244 -21.06082 -20.49919 -19.93757
-19.37595 -18.81433 -18.25271 -17.69109 -17.12946
-16.56784 -16.00622 -15.44460 -14.88298 -14.32136
-13.75973 -13.19811 -12.63649 -12.07487 -11.51325
-10.95162 -10.39000 -9.82838 -9.26676 -8.70514
-8.14352 -7.58189 -7.02027 -6.45865 -5.89703
-5.33541 -4.77379 -4.21216 -3.65054 -3.08892
-2.52730 -1.96568 -1.40405 -0.84243 -0.28081
0.28081 0.84243 1.40405 1.96568 2.52730
3.08892 3.65054 4.21216 4.77379 5.33541
5.89703 6.45865 7.02027 7.58189 8.14352
8.70514 9.26676 9.82838 10.39000 10.95162
11.51325 12.07487 12.63649 13.19811 13.75973
14.32136 14.88298 15.44460 16.00622 16.56784
17.12946 17.69109 18.25271 18.81433 19.37595
19.93757 20.49919 21.06082 21.62244 22.18406
22.74568 23.30730 23.86892 24.43055 24.99217
25.55379 26.11541 26.67703 27.23865 27.80028
28.36190 28.92352 29.48514 30.04676 30.60838
31.17000 31.73163 32.29325 32.85487 33.41649
33.97811 34.53973 35.10136 35.66298 36.22460
36.78622 37.34784 37.90946 38.47108 39.03270
39.59433 40.15595 40.71757 41.27919 41.84081
42.40243 42.96405 43.52567 44.08730 44.64892
45.21054 45.77216 46.33378 46.89540 47.45702
48.01864 48.58026 49.14189 49.70351 50.26513
50.82675 51.38837 51.94999 52.51161 53.07323
53.63485 54.19647 54.75809 55.31971 55.88133
56.44295 57.00457 57.56619 58.12781 58.68943
59.25105 59.81267 60.37429 60.93591 61.49753
62.05915 62.62077 63.18239 63.74401 64.30563
64.86725 65.42886 65.99048 66.55210 67.11372
67.67534 68.23695 68.79857 69.36019 69.92181
70.48342 71.04504 71.60665 72.16827 72.72988
73.29150 73.85311 74.41473 74.97634 75.53795
76.09956 76.66117 77.22278 77.78439 78.34600
78.90760 79.46920 80.03080 80.59240 81.15400
81.71559 82.27718 82.83876 83.40033 83.96190
84.52345 85.08499 85.64651 86.20800 86.76944
87.33080 87.89203 88.45297 89.01318 89.57009
ZDEF 1 LEVELS 1000.00000
TDEF 780 LINEAR 1JAN1950 1MO
VARS 59
TA 0 99 Surface Air Temperature at 2m K 1
TGEF 0 99 Effective Ground temperature (Radiation) K 1
SLP 0 99 SEA LEVEL PRESSSURE Pa 1
PS 0 99 SURFACE PRESSURE Pa 1
UA 0 99 Surface Zonal Velocity at 10m m/s 1
VA 0 99 Surface Merid. Velocity at 10m m/s 1
WIND 0 99 Surface Air Wind Speed at 10m m/s 1
RHA 0 99 Surface Air Relative Humidity at 2m % 1
QA 0 99 Surface Air Specific Humidity at 2m kg/kg 1
PRECIPI 0 99 Total Precipitation kg/m**2/s 1
SNP 0 99 Snow Precipitation kg/m**2/s 1
PPCI 0 99 Convective precipitation kg/m**2/s 1
EVSPS 0 99 water vapor flux kg/m**2/s 1
UMOM 0 99 momentum flux (X) N/m**2 1
VMOM 0 99 momentum flux (Y) N/m**2 1
FLLH 0 99 latent heat flux W/m**2 1
FLSH 0 99 sensible heat flux W/m**2 1
DLWB 0 99 Downard Longwave Radiation at the Bot W/m**2 1
ULWB 0 99 Upward Longwave Radiation at the Bot W/m**2 1
DSWB 0 99 Downward Shortwave Radiation at the Bot W/m**2 1
USWB 0 99 Upward Shortwave Radiation at the Bot W/m**2 1
CSDSWB 0 99 Downward Shortwave Radiation at the Bot (Clear Sky) W/m**2 1
CSUSWB 0 99 Upward Shortwave Radiation at the Bot (Clear Sky) W/m**2 1
CSDLWB 0 99 Downard Longwave Radiation at the Bot (Clear Sky) W/m**2 1
DSWT 0 99 Downward Shortwave Radiation at the Top W/m**2 1
USWT 0 99 Upward Shortwave Radiation at the Top W/m**2 1
ULWT 0 99 OLR (Upward Longwave Radiation at the Top) W/m**2 1
CSULWT 0 99 OLR clear sky (Upward Longwave Radiation at the Top) W/m**2 1
CSUSWT 0 99 Upward Shortwave Radiation at the Top (Clear Sky) W/m**2 1
PWATER 0 99 Precipitable Water kg/m**2 1
TCLOUD 0 99 Total cloud amount % 1
TCWC 0 99 Total cloud water content kg/m**2 1
VINTQU 0 99 column total of QU kg/kg m/s Pa 1
VINTQV 0 99 column total of QV kg/kg m/s Pa 1
AICE 0 99 Area fraction of sea ice % 1
WSL010 0 99 H2O SOIL upper 10cm kg/m**2 0
H2OSLT 0 99 H2O SOIL (total) kg/m**2 0
ROFS 0 99 surface runoff kg/m**2/s 1
ROF 0 99 Total Runoff kg/m**2/s 1
EVDWVEG 0 99 evap/dew on leaf (downward) kg/m**2/s 1
EVDWSL 0 99 evap/dew on soil (downward) kg/m**2/s 1
TRNSL 0 99 transpiration from soil (downward) kg/m**2/s 1
TMPSL1 0 99 TMP SOIL L1 K 0
TMPSL2 0 99 TMP SOIL L2 K 0
TMPSL3 0 99 TMP SOIL L3 K 0
TMPSL4 0 99 TMP SOIL L4 K 0
H2OSL1 0 99 H2O SOIL L1 kg/m**2 0
H2OSL2 0 99 H2O SOIL L2 kg/m**2 0
H2OSL3 0 99 H2O SOIL L3 kg/m**2 0
CVRSNWA 0 99 Snow Coverage 0-1 0
SWE 0 99 Snow water equivalent kg/m**2 1
DEPSNW 0 99 Snow Depth * CVRSNWA M 0
TMPSNW 0 99 TMP SNOW (vertical ave)* CVRSNWA K 0
EVDWSN 0 99 evap/dew on snow (downward) kg/m**2/s 1
SN2SL 0 99 snow melt water from snow to soil (downward) kg/m**2/s 1
YICE 0 99 Mass of Sea Ice kg/m**2 1
YSNW 0 99 Mass of snow on sea ice kg/m**2 1
TOTALHP 0 99 Total heating due to physics processes W/m**2 1
TOTALHM 0 99 Total Heating due to moist process J/m**2 1
ENDVARS
xgrads.CtlDescriptor(file="")
ここから、
- x方向のグリッド数：640
- y方向のグリッド数：320
- 要素数：59
であることが読み取れます。
この情報から、「緯度方向の格子座標の配列」と「変数の内容をまとめたデータフレーム」を作成しておきます。
# ライブラリをインポート import os import pandas as pd # ctlファイルのYDEFとVARSが記載されている行を特定 os.chdir("/mnt/c/Users/hoge/ONEFIFTY") substr_df = pd.read_fwf("./dias/data/GCM60_NHRCM20_150yr_TOUGOU/GCM60/HPD/sfc_avr_mon.ctl", widths=[10],header=None,skip_blank_lines=False).fillna("#") YDEF_ind = substr_df.loc[substr_df[0].str.startswith("YDEF")].index.values[0] YDEF_nrow = substr_df.loc[substr_df[0].str.startswith("ZDEF")].index.values[0] - YDEF_ind - 1 VARS_ind = substr_df.loc[substr_df[0].str.startswith("VARS")].index.values[0] VARS_nrow = substr_df.loc[substr_df[0].str.startswith("ENDVARS")].index.values[0] - VARS_ind - 1 # 格子の緯度情報を読み込み FLAT_vec = pd.read_fwf("./dias/data/GCM60_NHRCM20_150yr_TOUGOU/GCM60/HPD/sfc_avr_mon.ctl", header=None,skiprows=YDEF_ind+1,nrows=YDEF_nrow, colspecs=[(6,16),(17,27),(28,38),(39,49),(50,60)]).stack().values # 要素情報を読み込み parm_df = pd.read_fwf("./dias/data/GCM60_NHRCM20_150yr_TOUGOU/GCM60/HPD/sfc_avr_mon.ctl", header=None,skiprows=VARS_ind+1,nrows=VARS_nrow, colspecs=[(0,9),(25,100),(90,104),(107,108)],names=["varname","description","unit","flag"])
地上気温の月平均値の配列の作成
まずは地上気温の月平均値の3次元配列（time,y,x）を作成します。
# ライブラリをインポート 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") drlist_df = pd.DataFrame({"dr_path":list(pathlib.Path("./").glob("**/sfc_avr_mon_*.dr"))}) drlist_df["dr_name"] = drlist_df["dr_path"].map(lambda x: x.name) drlist_df["date"] = pd.to_datetime(drlist_df["dr_name"].str.extract(".*_(\d{6}).dr")[0],format="%Y%m") drlist_df = drlist_df.sort_values("date").reset_index(drop=True) # 空のNumpyArrayを用意する TA_arr = np.zeros([1800,320,640],np.float32) # 地上気温を切り出し、NumpyArrayに代入 for h_index,h_row in drlist_df.iterrows(): mon_dat = np.fromfile(str(h_row["dr_path"]),dtype='>f4') mon_arr = mon_dat.reshape(59,320,640) TA_arr[h_index,:,:] = mon_arr[0,:,:] # 子午線が中心に来るように変換 TA_arr = np.concatenate([TA_arr[:,:,320:],TA_arr[:,:,:320]],axis=2)
Xarrayの作成
ひとまず、Xarray形式に変換し、netCDF形式で保存します。
# lon軸の設定 FLON_vec = (np.arange(640) -320 ) * 0.56250 # xarray.Datasetを作成 TA_ds = xr.Dataset( { "temperature": (["time", "y", "x" ],TA_arr) }, coords={ "y": ("y", FLAT_vec), "x": ("x", FLON_vec), "time": drlist_df["date"], "reference_time": pd.Timestamp("1950-01-01"), }, ) # netCDFへ保存 TA_ds.to_netcdf("GCM60_RCP85_tmp_sfc_avr_mon.nc")
年平均値の時系列を作成
保存した配列は延べ1800ヵ月分、1.44ギガあり、（私の環境では）大きすぎて、そのままQGISで扱えないため、年平均値を計算して、150年分に縮小します。
# ライブラリをインポート import os import pathlib import numpy as np import pandas as pd import xarray as xr import affine import rioxarray # データの読み込み TA_ds = xr.open_dataset("GCM60_RCP85_tmp_sfc_avr_mon.nc") # 年平均値を計算 year_mean_ds = TA_ds.resample(time="A-DEC").mean(dim="time") # 投影情報を設定 year_mean_ds = year_mean_ds.rio.write_crs("+proj=longlat +a=6371000 +b=6371000 +over +no_defs") year_mean_ds = year_mean_ds.rio.write_transform(affine.Affine(0.5625,0,-180.2815,0,0.5616063492063492,-89.85701587301588)) # netCDFへ保存 year_mean_ds.to_netcdf("GCM60_RCP85_tmp_sfc_avr_year.nc")
QGISでアニメーション
QGISのメッシュレイヤに読み込み、アニメーションしてみます。
投影法の登録
このデータで用いられている投映法は、これまでに使っていたカスタム投影法と中央子午線が異なるため、新たに登録します。
QGISを立ち上げ、「設定」→「カスタム投影法」で”オプション－ユーザ定義CRS”を開きます。
＋ボタン（CRSを追加）をクリックし、次のように設定します
- 名前：GCM60_longlat
- 形式：Proj文字列（非推奨）
- パラメータ：+proj=longlat +a=6371000 +b=6371000 +over +no_defs
「検証」をクリックして確認してから、「OK」で登録します。
QGISに読み込む
「レイヤ」→「レイヤを追加」→「メッシュレイヤを追加」
ソース：メッシュデータセットに"GCM60_RCP85_tmp_sfc_avr_year.nc"を指定して「追加」
QGIS上で時系列アニメーションできます
他の気象要素についても同様に作成することができます。
「全球及び日本域150年連続実験データ」の全球60km150年連続実験データを入手し、QGISでアニメーション表示する方法を紹介しました。試してみてくださいね～。
※本記事では文部科学省「統合的気候モデル高度化研究プログラム」において、地球シミュレータを用いて作成されたデータを使用しました。
