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

気候予測データについて

はじめに

「日本域気候予測データ」には、バイアス補正をしていない「モデル格子点データ」と気象庁アメダス観測地点に対応するモデル格子点についてバイアス補正した「観測地点データ」の 2種類があります。

「観測地点データ」については、前の記事でマップ化、グラフ化する方法を紹介しました。

cci-labo.hateblo.jp

cci-labo.hateblo.jp

cci-labo.hateblo.jp

「モデル格子点データ」は、バイナリ形式（FORTRAN unformatted record）で公開されており、気候予測データセット2022解説書では、”格子点値解析描画ツール GrADS”と”GPV 読み込みマクロ”の２つのツールが紹介されています。（気候予測データセット解説書第２章．P2-139）

このデータをGeoTiffなどのラスタフォーマットに変換することにより、QGISなどの一般的なGISソフトウェアで扱うことが可能になります。

では、やってみましょう。

データの入手

観測地点データを入手した時と同様に、日本域気候予測データのダウンロード一覧ページにログインします。


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

  1. /jmagwp/gwp9/data/binary/SFA_rcp85_ens/annuallyフォルダを指定し、ファイル検索→全ファイルにチェックを入れる

  2. /jmagwp/gwp9/data/binary/SFA_rcp26_ens/annuallyフォルダを指定し、ファイル検索→全ファイルにチェックを入れる

  3. /jmagwp/gwp9/data/binary/SPA/annuallyフォルダを指定し、ファイル検索→全ファイルにチェックを入れる

  4. 一括ダウンロード

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

データを解凍する

# ユーザホームディレクトリに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\binary\SFA_rcp85_ens\annually
C:\Users\hoge\JWP9\dias\data\jmagwp\gwp9\data\binary\SFA_rcp26_ens\annually
C:\Users\hoge\JWP9\dias\data\jmagwp\gwp9\data\binary\SPA\annually

GrADSで表示

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

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

#GrDASを立ち上げ
grads

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

グリッドのデータ補完処理にしばらく時間がかかってから表示されます。

詳しい使い方は、このあたりを参考にさせていただいています。

seesaawiki.jp

wind.gp.tohoku.ac.jp

geotiffに変換する

気象予測の格子点データと人口の将来予測データや地形データを組み合わせることで、より効果的な解析が可能です。そのためには、気象予測データをGISに読み込める形式に変換する必要があります。

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

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

格子点データをGrADSに読み込むときに使った"CTLファイル”には、本体のバイナリデータに、どの範囲の、どのようなデータが、どのような形で収録されているかの情報が記載されています。

普通のテキストファイルですので、VSCodeで開くことができます。

DSET ^annual-SPA.dat

OPTIONS little_endian 
title surface result
undef 65535
pdef 467 744 LCCR 44 144 214 616 30 60 80 5000 5000
xdef  920 linear 114  0.05
ydef  800 linear 16   0.05
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***

このうち"pdef"行にバイナリデータのグリッド設定が記載されています。

その文法は以下を参照ください。

cola.gmu.edu

ctlの記載から、"annual-SPA.dat"には

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

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

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

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

VScodeの対話型ウインドウを起動します。

cci-labo.hateblo.jp

データを読み込み、変数別の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/binary/SPA/annually/annual-SPA.dat",dtype='<u2')

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

# 気象要素に切り分け、実数に変換
#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 = 3367936.9987040577
ul_y = 7752086.241150136
x_res = 5000
y_res = -5000

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_5km.gtiff"
output = gdal.GetDriverByName('GTiff').Create(name, 467, 744,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

QGISに読み込む

作成したGeoTIFFをQGISで表示してみます。

投影法の登録

格子点データのコントロールファイルのPDEF情報から、 緯度60°N、緯度 30°N の両緯線を2標準緯線とし、中央子午線を80°Eとしたランベルト正角円錐図法で投影されていることがわかります。

QGISでただしく表示するために、カスタム投影法として登録しておきます。

QGISを立ち上げ、「設定」→「カスタム投影法」で”オプション－ユーザ定義CRS”を開きます。

＋ボタン（CRSを追加）をクリックし、次のように設定します

  • 名前：NHRCM_LCCR
  • 形式：Proj文字列（非推奨）
  • パラメータ：+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

「検証」をクリックして確認してから、「OK」で登録します。

geotiffの読み込み

「レイヤ」→「レイヤを追加」→「ラスタレイヤを追加」

ソース：ラスタデータセットに"SPA_annually_PREC_mean_5km.gtiff"を指定して「追加」

年降水量のマップが表示されます。

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

平均気温など、他の気象要素についても同様に作成することができます。

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

※本記事では「気候予測データセット2022 ②日本域気候予測データ」を使用しました。同データセットは、気象庁気象研究所が開発した気候モデルを利用して、文部科学省気候変動リスク情報創生プログラム（RCP8.5シナリオ）及び統合的気候モデル高度化研究プログラム（RCP2.6シナリオ）において計算されたデータを元に作成されたものです。