「全球及び日本域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”コンテナの作成方法は以下になります。

cci-labo.hateblo.jp

# ライブラリをインポート
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