はじめに
「全球及び日本域気候予測データ」創生・統合プログラム5km格子NHRCM日本域気候予測データをラスタ化する際の投影パラメータを同定しました。
データの入手とマップ化の記事はこちら。
データ自体はGRIB2形式で提供されており、うまく設定すればそのままQGISで表示できるのですが、このデータには投影パラメータが設定されていないようで、そのままQGISに読み込むと、
一応読み込むことはできますが、このままでは日本列島の位置が合いません。
そこで、「日本域気候予測データ－格子点データ」5㎞版の時と同様に、投影パラメータを同定しました。
モデル格子点情報の入手
気候予測データセット（DS2022）から、創生・統合プログラム5km格子NHRCM日本域気候予測データの定数データを入手します。
ファイルは次の2つ
- cnst.dat
- cnst.ctl
ctlファイルの情報
GrADSで解析・表示するためのコントロールファイルは次のとおり。
***output grads ctl file for imt_cf***
dset cnst.dat
title constatnt data
OPTIONS big_endian sequential
undef -999.9
xdef 527 linear 1 1
ydef 804 linear 1 1
zdef 1 linear 1 1
tdef 1 linear 00:00z20JUL1980 1mn
vars 11
ZS 1 99 ZS
SL 1 99 SL
FLAT 1 99 FLAT
FLON 1 99 FLON
KIND 1 99 KIND
ROUGH 1 99 ROUGH
WET 1 99 WET
ALBED 1 99 ALBED
FKTG 1 99 FKTG
ROCTG 1 99 ROCTG
URBAN 1 99 URBAN
endvars
***end of ctl file***
ctlの記載から、このファイルは
- x方向のグリッド数：527
- y方向のグリッド数：804
- z方向のグリッド数：1
- タイムステップ数：1
- 変数数：11
であることが読み取れます。
データ形式は、気候予測データセット2022解説書のP2-56”モデル諸元”に記載があり、"32bit 実数バイナリ"とのこと。
Projectionの同定
座標系は「日本域気候予測データー格子点データ」と同じ、緯度60°N、緯度 30°N の両緯線を2標準緯線とし、中央子午線を80°Eとしたランベルト正角円錐図法を使います。
PROJ.4文字列にすると、
+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
になります。
GeoTransFormパラメータの同定
5km格子なので、resolutionとrotationのパラメータは
GT(2) = GT(4) = 0
GT(1) = 5000
GT(5) = -5000
になります。
upper-left pixelの座標はcnst.datから特定できます。
# ライブラリの呼び出し
import os
import numpy as np
from osgeo import gdal, gdalconst
from osgeo import osr
os.chdir("/mnt/c/Users/hoge/MHRCM")
# 4バイト実数として読み込み、変数毎に分割
cnst_vec = np.fromfile("./dias/data/NHRCM05_SOUSEI/cnst.dat",dtype='>f4')
cnst_arr = cnst_vec.reshape(11,423710)
# 緯度と経度情報を配列に変換
FLAT_arr = cnst_arr[2,1:-1].reshape(804,527)
FLON_arr = cnst_arr[3,1:-1].reshape(804,527)
# 左上格子の緯度経度を取り出す
ul_lat = FLAT_arr[-1,0].astype(float)
ul_lon = FLON_arr[-1,0].astype(float)
print(ul_lat,ul_lon)
#→57.08396 139.28857
ランベルト正角円錐図法の座標系に変換し、右上角の座標を同定します。
# ライブラリの呼び出し
from osgeo import gdal, gdalconst, gdal_array
from osgeo import osr
#緯度経度座標系のインスタンスを作成
src_srs = osr.SpatialReference()
src_srs.ImportFromEPSG(4326)
#ランベルト正角円錐図法座標系のインスタンスを作成
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")
#緯度経度をランベルト正角円錐図法の座標系に変換
trans = osr.CoordinateTransformation(src_srs,outRasterSRS)
ul_x,ul_y,Z = trans.TransformPoint(ul_lat,ul_lon)
#左上の角の位置を計算
ul_x = ul_x - 2500
ul_y = ul_y + 2500
print(ul_x,ul_y)
結果、GeoTransFormのパラメータは
[3217950.058222905,5000,0,7902104.780171389,0,-5000]
となります。
海陸比データをGeoTIFF化
同定したパラメータの検証のため、海陸比データをGeoTiff化します。
# ライブラリの呼び出し
import os
import numpy as np
from osgeo import gdal, gdalconst, gdal_array
from osgeo import osr
os.chdir("/mnt/c/Users/hoge/MHRCM")
#海陸比情報を読み込む
cnst_vec = np.fromfile("./dias/data/NHRCM05_SOUSEI/cnst.dat",dtype='>f4')
cnst_arr = cnst_vec.reshape(11,423710)
sl_arr = cnst_arr[1,1:-1].reshape(804,527)
#GeoTiFFを作成
name = "cnst_5km_sl.gtiff"
output = gdal.GetDriverByName('GTiff').Create(name, 527, 804 ,1,gdal.GDT_Float64)
#Projectionを設定
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")
output.SetProjection(outRasterSRS.ExportToWkt())
#GeoTransformを設定
output.SetGeoTransform([3217950.058222905,5000,0,7902104.780171389,0,-5000])
#海陸比データを設定
outband = output.GetRasterBand(1)
outband.WriteArray(sl_arr[::-1,:].astype("float64"))
outband.FlushCache()
output =None
QGISに読み込んでみる
うまく投影変換されていることが確認できます。
結果
気候予測データセット2022の「全球及び日本域気候予測データ」創生・統合プログラム5km格子NHRCM日本域気候予測データをGDALでラスタ化する際の投影パラメータは、
+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”
[3217950.058222905,5000,0,7902104.780171389,0,-5000]
であると同定できました！
（個人）気候変化情報研究室では、できるだけ簡単に、将来、地域の気候がどのように変化するかを示すグラフ・地図を作成する方法を開発し、紹介しています。
興味のある方はぜひ読んでみてくださいね～。