「日本域気候予測データ―格子点データ」を可視化する　その２－2km版データの入手とマップ化

はじめに

「日本域気候予測データ」では、それまでの『「地球温暖化予測情報第 9 巻」データセット（「第 9巻」データセット）』に水平解像度 2km への力学的ダウンスケーリングを行った予測が追加されました。

「2km 版」データセットの計算で用いた NHRCM02 は、「5km 版」で用いられた NHRCM05 では考慮されていなかった都市化による影響を含めることができる都市モデルが新たに導入されるなど、様々な改良が加えられている（表 2.2）。（気候予測データセット解説書第２章．P2-74）

とのこと。

2km版データについて、５㎞版と同様にマップ化してみましょう。

データの入手

格子点データをダウンロードします。2km版のデータは”/jmagwp/gwp9/data_2km/binary”フォルダに格納されています。

5㎞版と同様に、SFA_rcp85_ens,SFA_rcp26_ens,SPAのannualデータを選択し、入手します。

データを解凍する

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

次のディレクトリにデータが解凍されます C:\Users\hoge\JWP9\dias\data\jmagwp\gwp9\data_2km\binary\SFA_rcp85_ens\annually
C:\Users\hoge\JWP9\dias\data\jmagwp\gwp9\data_2km\binary\SFA_rcp26_ens\annually
C:\Users\hoge\JWP9\dias\data\jmagwp\gwp9\data_2km\binary\SPA\annually

GrADSで表示

5km版と同様に、まずはGrADSで表示してみます。

#SPAのデータフォルダに移動
cd /mnt/c/Users/hoge/JWP9/dias/data/jmagwp/gwp9/data_2km/binary/SPA/annually

#GrDASを立ち上げ
grads

#データの読み込み
open annual_2km-SPA.ctl
#コントロールファイルの内容を表示
q file 1
#timestepを指定
set t 1
#グリッドを塗り潰して色別で表示。
set gxout grfill
#メルカトル図法
set mproj latlon
#マップを表示
d prec

表示されます。

geotiffに変換する

５㎞版と同様に20年の平均値を計算し、GeoTIFFに変換します。

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

annual_2km-SPA.ctlの内容

DSET ^annual_2km-SPA.dat

OPTIONS little_endian 
title surface result
undef 65535
pdef 485 1681 LCCR 44 144 178 1474 30 60 80 2000 2000 
xdef 2300 linear 114  0.02 
ydef 2000 linear 16   0.02 
zdef    1 linear 1 1
tdef 20 linear 01SEP1980 1yr
vars   6
TS         1 -1,40,2 (Average temperature at 2m) / 0.01 K
TShMAXd    1 -1,40,2 (Maximum temperature at 2m) / 0.01 K
TShMINd    1 -1,40,2 (Minimum temperature at 2m) / 0.01 K
PREC       1 -1,40,2 (Precipitation            )        mm
SNDEPhMAXy 1 -1,40,2 (Maximum snow depth       )        cm
SNFALL     1 -1,40,2 (Snow precipitation       )        cm
endvars
***end of ctl file***

ctlの記載から、このファイルは

x方向のグリッド数：485
y方向のグリッド数：1681
z方向のグリッド数：1
タイムステップ数：20
変数数：6
変数のデータ形式：2-byte unsigned integers

であることが読み取れます。

この情報を基に、バイナリデータを読み込み、整形していきます。

バイナリデータの読み込み

データを読み込み、変数別のnumpy配列に変換します。

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

# バイナリファイルを読み込み
os.chdir("/mnt/c/Users/hoge/JWP9")
dvec = np.fromfile("./dias/data/jmagwp/gwp9/data_2km/binary/SPA/annually/annual-SPA.dat",dtype='<u2')

# 4次元の配列に変換
darr = dvec.reshape(20,6,1681,485)

# 気象要素に切り分け、実数に変換
#TS_arr = darr[:,0,:,:] *0.01 -273.15 
#TShMAXd_arr = darr[:,1,:,:] *0.01 -273.15
#TShMINd_arr = darr[:,2,:,:] *0.01 -273.15
PREC_arr = darr[:,3,:,:]*1.0
#SNDEPhMAXy_arr = darr[:,4,:,:]*1.0
#SNFALL_arr = darr[:,5,:,:]*1.0

GeoTIFFへ変換

20年の平均値を計算し、GeoTIFF化します。

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

#地理座標系のパラメータを設定
ul_x = 4074196.551466600
ul_y = 7530619.71549598
x_res = 2000
y_res = -2000

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")

# 格子点の平均値を計算
out_arr = np.mean(PREC_arr,axis=0)[::-1,:]

# GeoTIFF化し、保存
name = "SPA_annually_PREC_mean_2km.gtiff"
output = gdal.GetDriverByName('GTiff').Create(name, 485,1681,1,gdal.GDT_Float64)
output.SetGeoTransform((ul_x, x_res, 0, ul_y, 0, y_res))
outband = output.GetRasterBand(1)
outband.WriteArray(out_arr.astype("float64"))
output.SetProjection(outRasterSRS.ExportToWkt())
outband.FlushCache()
output =None