「全球及び日本域気候予測データ」 - 創生・統合プログラム5km格子NHRCM日本域気候予測データを可視化する

気候予測データについて

はじめに

「全球及び日本域気候予測データ」は、創生・統合プログラムによる気候変動予測モデルの実験出力データで、 20㎞格子の全球気候情報と５㎞格子と２㎞格子の日本域気候情報が公開されています。

このうち日本域気候情報は３０分～１時間の間隔の出力値が公開されており、 これまでに紹介した「日本域気候予測データ」は、この値を日、月、季節、年単位に集計したものです。

「日本域気候予測データー格子点データ」と同様にマップ化してみましょう。

データの入手

ここでは、シミュレーション上における1993年6月29日のデータを入手します。

※日付は気候モデルシミュレーション上の仮想的な日付で、1993年6月29日の実際の天候と一致しません。

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

  2. データセット紹介→"全球及び日本域気候予測データ"に進む

  3. データダウンロード＿領域5㎞をクリック

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

  5. 創生・統合プログラム5km格子NHRCM日本域気候予測データセットのダウンロード一覧ページにログイン
    ※初めての場合、利用申請を行なってください。

  6. ディレクトリに”/NHRCM05_SOUSEI/SPA/1992”フォルダを指定し、キーワードに”19930629”を指定、”部分一致”を選択してファイル検索→”SPA_5km_1992_19930629_surf.grb2”にチェックを入れる

  7. キーワードに”ctl”を指定、”部分一致”を選択してファイル検索→”SPA_5km_1992_surf.ctl”にチェックを入れる

  8. キーワードに”idx”を指定、”部分一致”を選択してファイル検索→"SPA_5km_1992_surf.idx"にチェックを入れ、一括ダウンロード

  9. 利用規約に同意してダウンロードする

データを解凍する

# ユーザホームディレクトリに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\NHRCM05_SOUSEI\SPA\1992

GrADSで表示

GrADSでデータを読み込む際に必要な CTL ファイル、IDXファイルも一緒にダウンロードしたので、まずはGrADSで表示してみます。

#データフォルダに移動
cd /mnt/c/Users/hoge/MHRCM/dias/data/NHRCM05_SOUSEI/SPA/1992

#GrDASを立ち上げ
grads

#データの読み込み
open SPA_5km_1992_surf.ctl
#コントロールファイルの内容を表示
q file 1
#timestepを指定　（1992年7月20日0時から344日と3時間後）
set t 16520
#グリッドを塗り潰して色別で表示。
set gxout grfill
#メルカトル図法
set mproj latlon
#海岸線を描かない
set mpdraw off
#地表気温マップを表示
d tsfc

表示されます。
※軸ラベルがおかしな値になってしまっていますが、投影法がうまく設定できていないためです。

geotiffに変換する

Pythonを用いて、GISで使用されるラスターデータの代表的なファイル形式であるGeoTIFFに変換します。

コントロールファイルから情報を抽出

「日本域気候予測データ―格子点データ」と同様に、5km格子NHRCM日本域気候予測データにも、GrADSでデータを読み込むためのコントロールファイルが公開されています。

SPA_5km_1992_surf.ctlの内容は

dset ^../1992/SPA_5km_1992%y4%m2%d2_surf.grb2
index ^SPA_5km_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_5km_1992%y4%m2%d2_surf.grb2
* griddef=1:0:(527 x 807):grid_template=0:winds(N/S): lat-lon grid:(527 x 807) 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 804 linear 1.000000 1
xdef 527 linear  1.000000 1
tdef 19585 linear 00Z20JUL1992 30mn
*tdef 48 linear 00Z20JUL1992 30mn
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方向のグリッド数：804
  • y方向のグリッド数：744
  • タイムステップ：30分間隔
  • 変数数：17

でデータが収録されていることが読み取れます。

この情報から、変数の内容をまとめたデータフレームを作成しておきます。

# ライブラリをインポート
import os
import pandas as pd

# ctlファイルの変数の情報が記載されている行を特定
os.chdir("/mnt/c/Users/hoge/MHRCM")
substr_df = pd.read_fwf("./dias/data/NHRCM05_SOUSEI/SPA/1992/SPA_5km_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/NHRCM05_SOUSEI/SPA/1992/SPA_5km_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ファイルのヘッダーから収録データのレコード情報を抽出し、整理しておきます。

docker-pyの使い方は以下を参考にさせていただきました。

qiita.com

"wgrib2_docker”コンテナの作成方法は以下です。

cci-labo.hateblo.jp

# ライブラリをインポート
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_5km_1992_19930629_surf.grb2 -VT -var', workdir='/mnt/c/Users/hoge/MHRCM/dias/data/NHRCM05_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-02 06:00:00"の地上気温データを読み込み、GeoTIFF化します。

# ライブラリをインポート
from osgeo import gdal, gdalconst, gdal_array 
from osgeo import osr

# grb2ファイルをGDALDatasetオブジェクトとしてオープン
grb2_ds = gdal.Open("./dias/data/NHRCM05_SOUSEI/SPA/1992/SPA_5km_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 = 3217950.058222905
ul_y = 7902104.780171389
x_res = 5000
y_res = -5000
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_5km_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_5km_TSFC_1993060290300.gtiff"を指定して「追加」

地表気温のマップが表示されます。

表示投影法をWGS84/Pseudo-Mercatorに変更。

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

「全球及び日本域気候予測データ」の5km格子NHRCM日本域気候予測データをを入手し、マップ化するまでを紹介しました。 試してみてくださいね～。

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