「全球及び日本域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でアニメーション表示する方法を紹介しました。試してみてくださいね～。
※本記事では文部科学省「統合的気候モデル高度化研究プログラム」において、地球シミュレータを用いて作成されたデータを使用しました。
「全球及び日本域150年連続実験データ」を可視化する その３ー統合プログラム領域20km150年連続実験データの時系列データをQGISでアニメーション表示する。
はじめに
前回は、150年連続実験データの時系列データをgifアニメーションにする方法を紹介しました。
ただ、やはり、興味のあるエリア（ROI）を拡大したり、動画を止めたり、ほかの情報と重ね合わせたい！
それなら、時間軸を持ったnetCDF形式で出力し、QGISにメッシュレイヤとして読み込むことで実現できそう。
ということで、やってみましょう。
netCDFの出力にはXarrayを使います。なければインストールしておいてください。
時系列netCDFの作成
まずは、各年8月上旬の地上気温の平均値の3次元配列（time,y,x）を作成します。
# ライブラリをインポート 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) # 空のNumpyArrayを用意する temp_arr = np.zeros([149,155,191],np.float32) # 1年ずつ期間平均し、配列に代入 for h_index,h_row in griblist_df.iterrows(): ds = xr.open_dataset(h_row["grib_path"], engine="cfgrib") da = ds["t"] da_eary = da.sel(time=slice("{}-08-01 1:00:00".format(h_row["year"]), "{}-08-11 0:00:00".format(h_row["year"]))) da_eary_mean = da_eary.mean(dim="time") temp_arr[h_index,:,:] = da_eary_mean.data
netCDFファイルの作成
xarray.Datasetを作成し、投影パラメータを設定してnetCDF形式のファイルへ出力します。
# x軸、y軸の設定 xx_vec = np.arange(191) * 20000 + 10000 -1930000.1287495417 yy_vec = np.arange(155) * -20000 - 10000 + 5833346.794147097 # xarray.Datasetを作成 temp_ds = xr.Dataset( { "temperature": (["time", "y","x" ], temp_arr[:,::-1,:]) }, coords={ "y": ("y", yy_vec), "x": ("x", xx_vec), "time": pd.date_range("1951-08-01",freq="AS-AUG",periods=149), "reference_time": pd.Timestamp("1950-08-01"), }, ) # 投影情報を設定 temp_ds = temp_ds.rio.write_crs("+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") temp_ds = temp_ds.rio.write_transform() # netCDFへ保存 temp_ds.to_netcdf("NHRCM20_RCP85_tmp_early_Aug.nc")
QGISに読み込む
QGISのメッシュレイヤに読み込み、アニメーションしてみます。
「レイヤ」→「レイヤを追加」→「メッシュレイヤを追加」
ソース：メッシュデータセットに"NHRCM20_RCP85_tmp_early_Aug.nc"を指定して「追加」
地表気温のマップが表示されます
時系列コントローラーパネルを開く
「アニメーション時系列ナビを表示」→ステップを1年に設定して再生
QGIS上で時系列アニメーションできます
自由に拡大してアニメーションできます
他の月、気象要素についても同様に作成することができます。
「全球及び日本域150年連続実験データ」の日本域20km150年連続実験データの時系列データをQGISでアニメーション表示する方法を紹介しました。試してみてくださいね～。
※本記事では文部科学省「統合的気候モデル高度化研究プログラム」において、地球シミュレータを用いて作成されたデータを使用しました。
気候予測データの解析環境を構築する。その6―Pythonライブラリ：Xarrayを導入する。
はじめに
前回、150年連続データセットのgifアニメーションを（力技で？）作成してみました
が、やはりQGIS上でアニメーションをぱたぱたさせたい！自分の住む地域を拡大してみてみたい！ということで、方法をさがしていたところ、以下の動画を発見。
時間軸が設定された”netCDF形式”のデータをメッシュレイヤとして追加すると、時系列コントローラを使用してパタパタさせることができるらしい。
”GUNMA GIS GEEK”にもwgrib2をnetCDFに変換して、QGIS読み込んでアニメーションさせる方法が載っていました。
Python上で時間平均したマップをnetCDF形式で出力すれば、前回作ったgifアニメーションをQGIS上でぱたぱたさせることができそうです。
で、Python上でnetcdfを扱うライブラリは"netCDF4"と思ってたら、最近は”Xarray"らしい。
さらに、”cfgrib”エンジンがあれば直接gribを扱えるそう。
これは便利！ということで、解析環境にXarrayと関連ライブラリを追加インストールします。
Xarrayと関連するライブラリのインストール
eccodeのライブラリをaptでインストールしておきます。
sudo apt update sudo apt upgrade sudo apt install libeccodes-dev
関連するライブラリを入れてから、Xarrayをインストールします。
sudo pip3 install netcdf4 #recommended if you want to use xarray for reading or writing netCDF files sudo pip3 install cfgrib #Python interface to map GRIB files to the NetCDF Common Data Model following the CF Convention using ecCodes. sudo pip3 install nc-time-axis #for plotting cftime.datetime objects sudo pip3 install pooch #manages your Python library's sample data files sudo pip3 install xarray #N-D labeled arrays and datasets in Python sudo pip3 install rioxarray #geospatial xarray extension powered by rasterio
テスト
GRIB Data Exampleを試してみます。
# ライブラリを呼び出す import xarray as xr import matplotlib.pyplot as plt # テストデータをロード ds = xr.tutorial.load_dataset("era5-2mt-2019-03-uk.grib", engine="cfgrib") # cartpyでマップ化 import cartopy.crs as ccrs import cartopy fig = plt.figure(figsize=(10, 10)) ax = plt.axes(projection=ccrs.Robinson()) ax.coastlines(resolution="10m") plot = ds.t2m[0].plot( cmap=plt.cm.coolwarm, transform=ccrs.PlateCarree(), cbar_kwargs={"shrink": 0.6} ) plt.title("ERA5 - 2m temperature British Isles March 2019")
動きましたか？
これでnetCDF形式のファイルを扱う環境が整いました。 さあ、気象予測データの解析をはじめてみましょう！
※QGISのメッシュレイヤ機能について、ずいぶん前に追加されていたのですが、まったく追えていませんでした・・これからいろいろ試していきたいと思っています。
「全球及び日本域150年連続実験データ」を可視化する その２ー統合プログラム領域20km150年連続実験データの時系列変化を可視化する。
はじめに
「全球及び日本域150年連続実験データ」は気候変動がいつの段階でどのくらいまで進むか？その時間推移を把握することができます。
その変化を視覚的に捉えるため、マップの時系列アニメーションを作成してみましょう。
時系列ラスタの作成
前回入手した、RCP8.5シナリオの8月の149年分のデータを集計し、8月上旬（1日～10日）の気温の時系列GeoTIFFを作成します。
1951/8/1~1951/8/10の平均マップ
まずは1951年8月上旬（1日～10日）について、地上気温の平均マップを作ってみます。
# ライブラリをインポート import os import pathlib import numpy as np import pandas as pd import datetime import docker from osgeo import gdal, gdalconst, gdal_array import matplotlib.pyplot as plt import cartopy.crs as ccrs # "wgrib_docker”コンテナを起動する client = docker.from_env() wgrib_cont =client.containers.get("wgrib_dock_mntd") wgrib_cont.start() # 処理ファイルのパスを設定 ctl_path = pathlib.Path("/mnt/c/Users/hoge/ONEFIFTY/dias/data/GCM60_NHRCM20_150yr_TOUGOU/NHRCM20/RCP8.5/1950/surf_HPD_195108.ctl") grib_path = pathlib.Path("/mnt/c/Users/hoge/ONEFIFTY/dias/data/GCM60_NHRCM20_150yr_TOUGOU/NHRCM20/RCP8.5/1950/surf_HPD_195108.grib") # ctlファイルのVAR行を特定 substr_df = pd.read_fwf(ctl_path,widths=[10],header=None,skip_blank_lines=False).fillna("#") n_skip = substr_df.loc[substr_df[0].str.startswith("var")].index.values[0] + 1 # VARの内容を読み込み parm_df = pd.read_fwf(ctl_path,header=None,skiprows=n_skip,nrows=17, colspecs=[(0,8),(11,13),(21,100)], names=["varname","parm_no","description"], dtype={"varname":str, "parm_no":int, "description":str}) # gribファイルのヘッダ情報を読み込み res_code,res_out = wgrib_cont.exec_run('wgrib -min {}'.format(str(grib_path)), workdir=str(grib_path.parents[2]),stderr=False) # レコード情報を整理 reco_df = pd.DataFrame({"v_info":res_out.decode().split("\n")}) reco_df = reco_df.loc[reco_df["v_info"].str.match("^\d")] reco_df["band_no"] = reco_df["v_info"].str.extract("(^\d*):.*",expand=False).astype(int) reco_df["date_time"] = pd.to_datetime( ("19" + reco_df["v_info"].str.extract(".*:d=(\d*-d*):.*",expand=False)), format="%Y%m%d%H-%M") reco_df["parm_str"] = reco_df["v_info"].str.extract(".*:d=\d*-\d*:([A-Z]*):.*",expand=False) reco_df["parm_no"] = reco_df["v_info"].str.extract(".*kpds5=(\d*)",expand=False).astype(int) reco_df = reco_df.reset_index(drop=True) # parm_dfを結合する reco_df = reco_df.merge(parm_df) # 8/1~8/10の地上気温レコードを絞り込む reco_tmp = reco_df.loc[ (reco_df["varname"] == "tmp") & (reco_df["date_time"] >= datetime.datetime(1951,8,1,0,0)) & (reco_df["date_time"] <= datetime.datetime(1951,8,11,0,0)) ] # gribファイルをGDALデータセットとして読み込む grib_ds = gdal.Open(str(grib_path),gdalconst.GA_ReadOnly) # 1レコードずつ読み込み、平均値を求める tmp_arr_sum = np.zeros((155,191)) for h_ind,h_row in reco_tmp.iterrows(): tmp_arr_sum += grib_ds.GetRasterBand(h_row["band_no"]).ReadAsArray() tmp_arr_mean = tmp_arr_sum/reco_tmp.shape[0] # Cartpyで表示してみる。 fig = plt.figure(figsize=(10, 5)) ax = plt.axes(projection=ccrs.LambertConformal(central_longitude=135.0, central_latitude=0,standard_parallels=(30.0,60.0))) ax.imshow(tmp_arr_mean,origin="upper", transform=ccrs.LambertConformal(central_longitude=135.0, central_latitude=0,standard_parallels=(30.0,60.0)), extent=[-1930000, 1890000, 2733347, 5833347]) ax.coastlines() ax.gridlines()
cartpyの使い方は以下を参考にさせていただいています。
処理の関数化
この処理を149年分繰り返し、149バンドのGeoTIFFを作成します。 ループを回すために、まずはそれぞれの処理を関数化しておきます。
# ライブラリをインポート import os import pathlib import numpy as np import pandas as pd import datetime import docker from osgeo import gdal, gdalconst, gdal_array from osgeo import osr # ctlファイルを読み込みparameter_dfを返す関数 def res_parameter_df(ctl_path): # 行を特定 substr_df = pd.read_fwf(ctl_path,widths=[10],header=None,skip_blank_lines=False).fillna("#") n_skip = substr_df.loc[substr_df[0].str.startswith("var")].index.values[0] + 1 # 読み込み parm_df = pd.read_fwf(ctl_path,header=None,skiprows=n_skip,nrows=17, colspecs=[(0,8),(11,13),(21,100)], names=["varname","parm_no","description"], dtype={"varname":str, "parm_no":int, "description":str}) # parm_dfを戻す return parm_df # gribファイルのヘッダ情報を読み込み、records_dfを返す関数 def res_records_df(grib_path,year,parm_df): # ヘッダ情報を読み込み res_code,res_out = wgrib_cont.exec_run('wgrib -min {}'.format(str(grib_path)), workdir=str(grib_path.parents[2]),stderr=False) # レコード情報を整理 reco_df = pd.DataFrame({"v_info":res_out.decode().split("\n")}) reco_df = reco_df.loc[reco_df["v_info"].str.match("^\d")] reco_df["band_no"] = reco_df["v_info"].str.extract("(^\d*):.*",expand=False).astype(int) reco_df["date_time"] = pd.to_datetime((str(h_row["year"])[:2] + reco_df["v_info"].str.extract(".*:d=(\d*-\d*):.*",expand=False)), format="%Y%m%d%H-%M") reco_df["parm_str"] = reco_df["v_info"].str.extract(".*:d=\d*-\d*:([A-Z]*):.*",expand=False) reco_df["parm_no"] = reco_df["v_info"].str.extract(".*kpds5=(\d*)",expand=False).astype(int) reco_df = reco_df.reset_index(drop=True) #parm_dfを結合してから戻す return reco_df.merge(parm_df) # gribファイルから特定要素、期間のデータを切り出し、平均値の配列を返す関数 def res_mean_arr(grib_ds,reco_df,varname,datetime_start,datetime_end): reco_sub = reco_df.loc[ (reco_df["varname"] == varname) & (reco_df["date_time"] >= datetime_start) & (reco_df["date_time"] <= datetime_end) ] # 1レコードずつ読み込み、平均値を求める arr_sum = np.zeros((155,191)) for h_ind,h_row in reco_sub.iterrows(): arr_sum += grib_ds.GetRasterBand(h_row["band_no"]).ReadAsArray() arr_mean = arr_sum/reco_sub.shape[0] # 平均値の配列を戻す return arr_mean
時系列GeoTIFFの作成
149年分のデータについて、8月上旬（1日～10日）の地上気温の平均マップを作成し、1つのGeoTIFFにします。
# 前のコードブロックのライブラリの読み込みと関数を設定しておく。 # "wgrib_docker”コンテナを起動する client = docker.from_env() wgrib_cont =client.containers.get("wgrib_dock_mntd") wgrib_cont.start() # 地理座標系のパラメータを設定 ul_x = -1930000.1287495417 ul_y = 5833346.794147097 x_res = 20000 y_res = -20000 x_size = 191 y_size = 155 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") # 処理するファイルの一覧を作成する os.chdir("/mnt/d/Users/kazuh/ONEFIFTY") griblist_df = pd.DataFrame({"grib_path":list(pathlib.Path("./").glob("**/*.grib"))}) ctllist_df = pd.DataFrame({"ctl_path":list(pathlib.Path("./").glob("**/*.ctl"))}) griblist_df["dir_path"] = griblist_df["grib_path"].map(lambda x: x.parent.absolute()) ctllist_df["dir_path"] = ctllist_df["ctl_path"].map(lambda x: x.parent.absolute()) griblist_df = griblist_df.merge(ctllist_df) 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) # 空の出力ファイルを用意する name = "NHRCM20_RCP85_tmp_early_Aug.gtiff" output = gdal.GetDriverByName('GTiff').Create(name,x_size,y_size,149,gdal.GDT_Float64) output.SetProjection(outRasterSRS.ExportToWkt()) output.SetGeoTransform((ul_x, x_res, 0, ul_y, 0, y_res)) # ファイル一覧を1行ずつ処理 for h_index,h_row in griblist_df.iterrows(): #records_dfを作成 parameter_df = res_parameter_df(h_row["ctl_path"].absolute()) records_df = res_records_df(h_row["grib_path"].absolute(),h_row["year"],parameter_df) # gribファイルをGDALデータセットとして読み込む grib_ds = gdal.Open(str(h_row["grib_path"].absolute()),gdalconst.GA_ReadOnly) # 8/１ 1:00～8/11 0:00の地表気温の平均値マップを作成 tmp_arr = res_mean_arr(grib_ds,records_df,"tmp", datetime.datetime(h_row["year"],h_row["month"],1,0), datetime.datetime(h_row["year"],h_row["month"],11,0)) # GeiTIFFのバンドに保存 outband = output.GetRasterBand(h_index+1) outband.WriteArray(tmp_arr.astype("float64")) # バンド名を設定 outband.SetDescription("{}-{:02d}".format(h_row["year"],h_row["month"])) outband.FlushCache() # 出力ファイルを閉じる output = None
アニメーション化
matplotlibを用いてgifアニメーションを作成します。
# ライブラリをインポート import os import pathlib from osgeo import gdal, gdalconst, gdal_array import matplotlib.pyplot as plt from matplotlib import animation, rc import cartopy.crs as ccrs %matplotlib widget plt.rcParams["font.family"] = 'TakaoGothic' # geotiffを読み込む os.chdir("/mnt/d/Users/kazuh/ONEFIFTY") geotiff_ds = gdal.Open("NHRCM20_RCP85_tmp_early_Aug.gtiff",gdalconst.GA_ReadOnly) outband = geotiff_ds.GetRasterBand(1) # マップを更新する関数 def update(i): if i != 0: plt.cla() outband = geotiff_ds.GetRasterBand(i+1) ax.set_extent([110, 160, 19, 50]) ax.imshow(outband.ReadAsArray(),origin="upper",vmin=5.,vmax=40., transform=ccrs.LambertConformal(central_longitude=135.0,central_latitude=0, standard_parallels=(30.0,60.0)), extent=[-1930000, 1890000, 2733347, 5833347],cmap='afmhot') ax.coastlines(resolution="10m") ax.gridlines() ax.set_title('地上気温 {}'.format(outband.GetDescription())) return fig # アニメーションを作成 fig = plt.figure(figsize=(10, 5)) ax = plt.axes(projection=ccrs.Mercator()) ax.set_extent([110, 160, 19, 50]) im = ax.imshow(outband.ReadAsArray(),origin="upper",vmin=5.,vmax=40., transform=ccrs.LambertConformal(central_longitude=135.0,central_latitude=0, standard_parallels=(30.0,60.0)), extent=[-1930000, 1890000, 2733347, 5833347],cmap='afmhot') ax.coastlines(resolution="10m") ax.gridlines() cb = plt.colorbar(im) cb.ax.set_title('[℃]') ax.set_title('地上気温 {}'.format(outband.GetDescription())) ani = animation.FuncAnimation(fig, update, 148, interval=600) ani.save("NHRCM20_RCP85_tmp_early_Aug.gif", writer = 'imagemagick')
アニメーションができました!
他の月、気象要素についても同様に作成することができます。
「全球及び日本域150年連続実験データ」の20km150年連続実験データのアニメーションを作る方法を紹介しました。試してみてくださいね～。
※本記事では文部科学省「統合的気候モデル高度化研究プログラム」において、地球シミュレータを用いて作成されたデータを使用しました。
「全球及び日本域150年連続実験データ」を可視化する その１―統合プログラム領域20km150年連続実験データの入手とマップ化
はじめに
「全球及び日本域150年連続実験データ」は1950 年から 2099 年までの連続したシミュレーション実験の出力データで、60km格子の全球気候情報と20km格子の日本域気候情報が公開されています。
このデータを用いることにより。
’20 世紀中頃から 21 世紀末までの連続した実験を行うことで、極端な気象現象が過去から 21 世紀末にかけてどのような推移で変化していくか、温暖化による変化が自然変動の変動幅を超えて明瞭になる時期が、どういう指標でいつ頃にはっきりするのか等を調べることが可能となる。（気候予測データセット解説書第1章．P1-10より引用）
とのこと。20km格子の日本域気候情報について、まずはマップ化してみましょう。
データ入手
ここでは、RCP8.5シナリオの8月のデータを入手します。
気候予測データセット（DS2022）のホームページにアクセス
diasjp.netデータセット紹介→全球及び日本域150年連続実験データに進む
データダウンロードをクリック
"データをダウンロード"をクリック
統合プログラム全球60km領域20km150年連続実験データセットのダウンロード一覧ページにログイン
1951年、8月の地上大気データを検索し、全ファイルにチェックを入れる
- ディレクトリ指定：/GCM60_NHRCM20_150yr_TOUGOU/NHRCM20/RCP8.5/1950
- キーワード指定：surf 195108(部分一致）
翌年8月の地上大気データを検索し、全ファイルにチェックを入れる
- ディレクトリ指定：/GCM60_NHRCM20_150yr_TOUGOU/NHRCM20/RCP8.5/2051
- キーワード指定：surf 195208(部分一致)
これを繰り返し、５～10年ごとに一括ダウンロードする
データを解凍する
# ユーザホームディレクトリに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\NHRCM20\RCP8.5
GrADSで表示
GrADSでデータを読み込む際に必要な CTL ファイル, IDXファイルも一緒にダウンロードしたので、まずはGrADSで表示してみます。
#データフォルダに移動 cd /mnt/c/Users/hoge/ONEFIFTY/dias/data/GCM60_NHRCM20_150yr_TOUGOU/NHRCM20/RCP8.5/1950/ #GrDASを立ち上げ grads
#データの読み込み open surf_HPD_195108.ctl #コントロールファイルの内容を表示 q file 1 #timestepを指定 set t 3 #グリッドを塗り潰して色別で表示。 set gxout grfill #メルカトル図法 set mproj latlon #地表気温マップを表示 d tmp
表示されます。
geotiffに変換する
ためしに時間値をGeoTIFFに変換します。
コントロールファイルから情報を抽出
surf_HPD_195108.ctlの内容は
***output grads ctl file for surf*** dset surf_HPD_195108.grib index surf_HPD_195108.idx pdef 191 155 LCCR 35 135 97 77 30 60 135 20000 20000 dtype grib 255 undef -9.9999e+20 xdef 321 linear 105 0.2 ydef 205 linear 15 0.2 zdef 1 linear 1 1 tdef 744 linear 01:00z01AUG1951 60mn vars 17 smqr 0 59, 1, 0 accumulated rain [mm/h] smqi 0 58, 1, 0 accumulated ice [mm/h] smqs 0 64, 1, 0 accumulated snow [mm/h] smqg 0 78, 1, 0 accumulated graupel [mm/h] smqh 0 79, 1, 0 accumulated hail [mm/h] rain 0 61, 1, 0 precipitation [mm/h] psea 0 2, 1, 0 sea level pressure [hPa] psurf 0 1, 1, 0 surface pressure [hPa] u 0 33, 1, 0 u-component of wind [m/s] v 0 34, 1, 0 v-component of wind [m/s] tmp 0 11, 1, 0 temperature [K] ttd 0 18, 1, 0 dew point depression [K] cll 0 73, 1, 0 low cloud cover [0-1] clm 0 74, 1, 0 medium cloud cover [0-1] clh 0 75, 1, 0 high cloud cover [0-1] cla 0 71, 1, 0 total cloud cover [0-1] tpw 0 54, 1, 0 precipitable water [kg/m2] endvars ***end of ctl file***
ここから、
- x方向のグリッド数：155
- y方向のグリッド数：191
- タイムステップ：1時間間隔
- 変数数：17
でデータが収録されていることが読み取れます。
この情報から、変数の内容をまとめたデータフレームを作成しておきます。
# ライブラリをインポート import os import pandas as pd # ctlファイルの変数の情報が記載されている行を特定 os.chdir("/mnt/c/Users/hoge/ONEFIFTY") substr_df = pd.read_fwf("./dias/data/GCM60_NHRCM20_150yr_TOUGOU/NHRCM20/RCP8.5/1950/surf_HPD_195108.ctl",widths=[10],header=None,skip_blank_lines=False).fillna("#") n_skip = substr_df.loc[substr_df[0].str.startswith("var")].index.values[0] + 1 # 変数の情報が記載されている行を読み込み parm_df = pd.read_fwf("./dias/data/GCM60_NHRCM20_150yr_TOUGOU/NHRCM20/RCP8.5/1950/surf_HPD_195108.ctl",header=None,skiprows=n_skip,nrows=17,colspecs=[(0,8),(11,13),(21,100)],names=["varname","parm_no","description"],dtype={"varname":str, "parm_no":int, "description":str})
収録データのレコード情報を整理する
docker-pyを使って"wgrib_docker”コンテナに接続し、gribファイルのヘッダーから収録データのレコード情報を抽出し、整理しておきます。
"wgrib_docker”コンテナの作成方法は以下になります。
# ライブラリをインポート import docker import pandas as pd # "wgrib_docker”コンテナを起動する client = docker.from_env() wgrib_cont =client.containers.get("wgrib_dock") wgrib_cont.start() # gribファイルのヘッダ情報を入手 res_code,res_out = wgrib_cont.exec_run('wgrib -min ./1950/surf_HPD_195108.grib', workdir='/mnt/c/Users/kazuh/ONEFIFTY/dias/data/GCM60_NHRCM20_150yr_TOUGOU/NHRCM20/RCP8.5/',stderr=False) # レコード番号（ラスタのバンド番号）、日時、パラメータ番号を整理 records_df = pd.DataFrame({"v_info":res_out.decode().split("\n")}) records_df = records_df.loc[records_df["v_info"].str.match("^\d")] records_df["band_no"] = records_df["v_info"].str.extract("(^\d*):.*",expand=False).astype(int) records_df["date_time"] = pd.to_datetime(("19" + records_df["v_info"].str.extract(".*:d=(\d*-\d*):.*",expand=False)),format="%Y%m%d%H-%M") records_df["parm_str"] = records_df["v_info"].str.extract(".*:d=\d*-\d*:([A-Z]*):.*",expand=False) records_df["parm_no"] = records_df["v_info"].str.extract(".*kpds5=(\d*)",expand=False).astype(int) records_df = records_df.reset_index(drop=True) # ctlファイルから作成したparm_dfを結合しておく records_df = records_df.merge(parm_df)
レコードの読み込みとGeoTIFFへの変換
PyGDALを用いて"1951-08-01 06:00:00"の地上気温データを読み込み、GeoTIFF化します。
# ライブラリをインポート from osgeo import gdal, gdalconst, gdal_array from osgeo import osr # gribファイルをGDALDatasetオブジェクトとしてオープン grib_ds = gdal.Open("./dias/data/GCM60_NHRCM20_150yr_TOUGOU/NHRCM20/RCP8.5/1950/surf_HPD_195108.grib", gdalconst.GA_ReadOnly) # レコード情報から"1951-08-01 06:00:00"の地上気温データのラスタバンド番号を取り出す h_band = records_df.loc[(records_df["varname"] == "tmp") & (records_df["date_time"] == "1951-08-01 06:00:00"),"band_no"].values[0].astype(int) # GDALDatasetオブジェクトから対応するバンド番号のレイヤデータを切り出す tmp_arr = grib_ds.GetRasterBand(int(h_band)).ReadAsArray() # 地理座標系のパラメータを設定 ul_x = -1930000.1287495417 ul_y = 5833346.794147097 x_res = 20000 y_res = -20000 y_size,x_size = tmp_arr.shape 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") # GeoTIFF化し、保存 name = "surf_HPD_tmp_195108010600.gtiff" 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(tmp_arr.astype("float64")) outband.FlushCache() output = None
QGISに読み込む
作成したGeoTIFFをQGISで表示してみます。
投影法の登録
このデータで用いられている投影法は、これまでに使っていたカスタム投影法と中央子午線が異なるため、新たに登録します。
QGISを立ち上げ、「設定」→「カスタム投影法」で”オプション－ユーザ定義CRS”を開きます。
＋ボタン（CRSを追加）をクリックし、次のように設定します
- 名前：NHRCM20_LCCR
- 形式：Proj文字列（非推奨）
- パラメータ：+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
「検証」をクリックして確認してから、「OK」で登録します。
GeoTIFFの読み込み
「レイヤ」→「レイヤを追加」→「ラスタレイヤを追加」
ソース：ラスタデータセットに"surf_HPD_tmp_195108010600.gtiff"を指定して「追加」
地表気温のマップが表示されます
表示投影法をWGS84/Pseudo-Mercatorに変更
他の日時、気象要素についても同様に作成することができます。
「全球及び日本域150年連続実験データ」の20km150年連続実験データを入手し、マップ化するまでを紹介しました。試してみてくださいね～。
※本記事では文部科学省「統合的気候モデル高度化研究プログラム」において、地球シミュレータを用いて作成されたデータを使用しました。
気候予測データの解析環境を構築する。その５―Docker-wgrib環境を構築する。
気候予測データの解析環境を構築する。その４―Docker-wgrib2環境を構築する。では、GRIB2形式を読み込むための環境を構築しました。
気候予測データセット（DS2022）の「全球及び日本域150年連続実験データ」はGRIB1形式で公開されており、wgrib2からうまく読み込めないことが判明*1したため、GRIB1を扱うために"wgrib"のコンテナを作成します。
wgrib用コンテナを作成
前回利用したochonjo/wgrib2-with-gdalから、新たにwgrib用コンテナを作成します。
docker run --name "wgrib_dock" -v /mnt/c:/mnt/c -it ochonjo/wgrib2-with-gdal:3.4.2
wgrib_dockターミナルでwgribのソースコードを入手し、コンパイルしインストールします。
# ソースコードを入手 cd ~ wget https://ftp.cpc.ncep.noaa.gov/wd51we/wgrib/wgrib.tar # 解凍 mkdir wgrib tar -xvf ./wgrib.tar -C ./wgrib # コンパイル cd ./wgrib ./src2all alpha make # 実行ファイルをコピー cp wgrib /usr/local/bin/wgrib # 終了 exit
テスト
以下のページを参考にうまく動くかテストします。
https://ftp.cpc.ncep.noaa.gov/wd51we/wgrib/porting.txt
コンテナに接続する
docekrのサイドバーの"wgrib_dock"を右クリック→「start」
"wgrib_dock"を右クリック→「Attach shell」
"wgrib_dock"コンテナのターミナルが起動する
# 呼び出してみる
wgirb
# テストデータを入手 cd ~ wget http://ftp.cpc.ncep.noaa.gov/wd51we/wgrib/land.grb # ファイルを読み込んでみる wgrib land.grb # ヘッダ情報を表示 wgrib land.grb -V
動きましたか？
これでGRIB1形式のファイルを処理する環境が整いました。 さあ、気象予測データの解析をはじめてみましょう！
*1:私がwgrib2でGRIB1を読み込む方法を知らないだけ・・かもしれません。
「全球及び日本域気候予測データ」 創生・統合プログラム2km格子NHRCM日本域気候予測データを可視化する
はじめに
２㎞格子の日本域気候情報について、５㎞格子と同様にマップ化してみましょう。
データの入手
ここでは、シミュレーション上における1993年6月29日のデータを入手します。
※あくまでも気候モデルシミュレーション上の日付なので、実際の1993年6月29日の天候とは一致しません。
気候予測データセット（DS2022）のホームページにアクセス
diasjp.netデータセット紹介→"全球及び日本域気候予測データ"に進む
データダウンロード＿領域2㎞をクリック
"データをダウンロード"をクリック
創生・統合プログラム2km格子NHRCM日本域気候予測データセットのダウンロード一覧ページにログイン
※初めての場合、利用申請を行なってください。
ディレクトリに”/NHRCM02_SOUSEI/SPA/1992”フォルダを指定し、キーワードに”19930629”を指定、”部分一致”を選択してファイル検索→”SPA_2km_1992_19930629_surf.grb2”にチェックを入れる
キーワードに”ctl”を指定、”部分一致”を選択してファイル検索→”SPA_2km_1992_surf.ctl”にチェックを入れる
キーワードに”idx”を指定、”部分一致”を選択してファイル検索→"SPA_2km_1992_surf.idx"にチェックを入れ、一括ダウンロード
利用規約に同意してダウンロードする
データを解凍する
# ユーザホームディレクトリにMHRCMフォルダを作成し、 # ダウンロードした”data**.tgz”を移動しておく。 cd /mnt/c/Users/hoge/MHRCM/
find ./ -type f -name ".tgz" | xargs -n 1 tar zxf
次のディレクトリにデータが解凍されます
C:\Users\hoge\MHRCM\dias\data\NHRCM02_SOUSEI\SPA\1992
GrADSで表示
GrADSでデータを読み込む際に必要な CTL ファイル, IDXファイルも一緒にダウンロードしたので、まずはGrADSで表示してみます。
#データフォルダに移動 cd /mnt/c/Users/hoge/MHRCM/dias/data/NHRCM02_SOUSEI/SPA/1992#GrDASを立ち上げ grads
#データの読み込み open SPA_2km_1992_surf.ctl #コントロールファイルの内容を表示 q file 1 #timestepを指定 1992年7月20日0時から344日と3時間後 set t 8260 #グリッドを塗り潰して色別で表示。 set gxout grfill #メルカトル図法 set mproj latlon #海岸線を描かない set mpdraw off #地表気温マップを表示 d tsfc
表示されます。
※軸ラベルがおかしな値になってしまっていますが、これは投影法がうまく設定できていないためです。
GeoTIFFに変換する
Pythonを用いて、GISで使用されるラスターデータの代表的なファイル形式であるGeoTIFFに変換します。
コントロールファイルから情報を抽出
SPA_2km_1992_surf.ctlの内容は
dset ^../1992/SPA_2km_1992%y4%m2%d2_surf.grb2 index ^SPA_2km_1992_surf.idx undef 9.999E+20 title GRIB2 data of SPA_2km_1992_surf * produced by g2ctl v0.0.8.8 * command line options: SPA_2km_1992%y4%m2%d2_surf.grb2 * griddef=1:0:(525 x 1721):grid_template=0:winds(N/S): lat-lon grid:(525 x 1721) units 1e-06 input WE:SN output WE:SN res 48 lat 1.000000 to 1721.000000 by 1.000000 lon 1.000000 to 165.000000 by 1.000000 #points=903525:winds(N/S)dtype grib2 ydef 1721 linear 1.000000 1 xdef 525 linear 1.000000 1 tdef 9793 linear 00Z20JUL1992 1hr zdef 1 linear 1 1 options template vars 17 RR1 0,1 0,194,1 surface 1 HOUR Precipitation [mm] QR1 0,1 0,194,2 surface 1 HOUR Precipitation of Liquid Water [mm] QI1 0,1 0,194,3 surface 1 HOUR Precipitation of Ice [mm] QS1 0,1 0,194,4 surface 1 HOUR Precipitation of Snow [mm] QG1 0,1 0,194,5 surface 1 HOUR Precipitation of Graupel [mm] QH1 0,1 0,194,6 surface 1 HOUR Precipitation of Hail [mm] PSFC 0,1 0,194,7 surface Surface Pressure [hPa] PSEA 0,1 0,194,8 surface Sea Level Pressure [hPa] USFC 0,1 0,194,9 surface Velocity U at 10m [m/s] VSFC 0,1 0,194,10 surface Velocity V at 10m [m/s] TSFC 0,1 0,194,11 surface Temperature at 1.5m [K] TTDSFC 0,1 0,194,12 surface Dew Point Depression at 1.5m [K] CLL 0,1 0,194,13 surface Low Cloud Cover [%] CLM 0,1 0,194,14 surface Middle Cloud Cover [%] CLH 0,1 0,194,15 surface High Cloud Cover [%] CLA 0,1 0,194,16 surface Total Cloud Cover [%] TPW 0,1 0,194,17 ** surface Precipitable Water [kg/m2] ENDVARS
ここから、
- x方向のグリッド数：1721
- y方向のグリッド数：525
- タイムステップ：1時間間隔
- 変数数：17
であることが読み取れます。
この情報から、変数の内容をまとめたデータフレームを作成しておきます。
# ライブラリをインポート import os import pandas as pd # ctlファイルの変数の情報が記載されている行を特定 os.chdir("/mnt/c/Users/hoge/MHRCM") substr_df = pd.read_fwf("./dias/data/NHRCM02_SOUSEI/SPA/1992/SPA_2km_1992_surf.ctl",widths=[10],header=None,skip_blank_lines=False).fillna("#") n_skip = substr_df.loc[substr_df[0].str.startswith("var")].index.values[0] + 1 # 変数の情報が記載されている行を読み込み parm_df = pd.read_fwf("./dias/data/NHRCM02_SOUSEI/SPA/1992/SPA_2km_1992_surf.ctl",header=None,skiprows=n_skip,nrows=17,colspecs=[(0,8),(27,31),(34,100)],names=["varname","parm_no","description"])
収録データのレコード情報を整理する
docker-pyを使って"wgrib2_docker”コンテナに接続し、grb2ファイルのヘッダーから収録データのレコード情報を抽出し、整理しておきます。
# ライブラリをインポート import docker import pandas as pd # "wgrib2_docker”コンテナを起動する client = docker.from_env() wgrib2_cont =client.containers.get("wgrib2_dock") wgrib2_cont.start() # grb2ファイルのヘッダ情報を入手 res_code,res_out = wgrib2_cont.exec_run('wgrib2 SPA_2km_1992_19930629_surf.grb2 -VT -var', workdir='/mnt/c/Users/hoge/MHRCM/dias/data/NHRCM02_SOUSEI/SPA/1992/',stderr=False) # レコード番号（ラスタのバンド番号）、日時、パラメータ番号を整理 records_df = pd.DataFrame({"v_info":res_out.decode().split("\n")}) records_df = records_df.loc[records_df["v_info"].str.match("^\d")] records_df["band_no"] = records_df["v_info"].str.extract("(^\d*):.*",expand=False).astype(int) records_df["date_time"] = pd.to_datetime(records_df["v_info"].str.extract(".*:vt=(\d*):.*",expand=False),format="%Y%m%d%H%M%S") records_df["parm_no"] = records_df["v_info"].str.extract(".*parm=(\d*)",expand=False).astype(int) records_df = records_df.reset_index(drop=True) # ctlファイルから作成したparm_dfを結合しておく records_df = records_df.merge(parm_df)
レコードの読み込みとGeoTIFFへの変換
PyGDALを用いて"1993-06-29 03:00:00"の地上気温データを読み込み、GeoTIFF化します。
# ライブラリをインポート from osgeo import gdal, gdalconst, gdal_array from osgeo import osr # grb2ファイルをGDALDatasetオブジェクトとしてオープン grb2_ds = gdal.Open("./dias/data/NHRCM02_SOUSEI/SPA/1992/SPA_2km_1992_19930629_surf.grb2", gdalconst.GA_ReadOnly) # レコード情報から1993-06-29 03:00:00 の地上気温データのラスタバンド番号を取り出す h_band = records_df.loc[(records_df["varname"] == "TSFC") & (records_df["date_time"] == "1993-06-29 03:00:00"),"band_no"].values[0].astype(int) # GDALDatasetオブジェクトから対応するバンド番号のレイヤデータを切り出す TSFC_arr = grb2_ds.GetRasterBand(int(h_band)).ReadAsArray() # 地理座標系のパラメータを設定 ul_x = 4034196.596039401 ul_y = 7570619.639166898 x_res = 2000 y_res = -2000 y_size,x_size = TSFC_arr.shape outRasterSRS = osr.SpatialReference() outRasterSRS.ImportFromProj4( "+proj=lcc +lat_1=30.0 +lat_2=60.0 +lon_0=80. +lat_0=0 +x_0=0 +y_0=0 +a=6371000 +b=6371000 +units=m +no_defs") # GeoTIFF化し、保存 name = "SPA_2km_TSFC_1993060290300.gtiff" 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(TSFC_arr.astype("float64")) outband.FlushCache() output =None
QGISに読み込む
「レイヤ」→「レイヤを追加」→「ラスタレイヤを追加」
ソース：ラスタデータセットに""SPA_2km_TSFC_1993060290300.gtiff"を指定して「追加」
地表気温のマップが表示されます
表示投影法をWGS84/Pseudo-Mercatorに変更。
他の日時、気象要素についても同様に作成することができます。
「全球及び日本域気候予測データ」の2km格子NHRCM日本域気候予測データをを入手し、マップ化するまでを紹介しました。 試してみてくださいね～。
※本記事では文部科学省「リスク情報創生プログラム」及び「統合的気候モデル高度化プログラム」において、地球シミュレータを用いて作成されたデータを使用しました。
