はじめに

「全球及び日本域気候予測データ」創生・統合プログラム2km格子NHRCM日本域気候予測データをラスタ化する際の投影パラメータを同定しました。

「全球及び日本域気候予測データ」創生・統合プログラム2km格子NHRCM日本域気候予測データの入手とマップ化の記事はこちら。

GRIB2形式のファイルをそのままQGISに読み込むと、

5km格子NHRCM日本域気候予測データと同様に、一応読み込むことはできますが、このままでは日本列島の位置が合いません。

そこで、5㎞格子の時と同様に、投影パラメータを同定しました。

モデル格子点情報の入手

気候予測データセット（DS2022）から、創生・統合プログラム2km格子NHRCM日本域気候予測データの定数データを入手します。

ファイルは次の2つ

cnst2km.dat
cnst2km.ctl

ctlファイルの情報

GrADSで解析・表示するためのコントロールファイルは次のとおり。

cnst2km.ctl

***output grads ctl file for imt_cf***
dset cnst2km.dat
title constatnt data
OPTIONS big_endian sequential 
undef -999.9
xdef  525 linear 1 1
ydef 1721 linear 1 1
zdef    1 linear 1 1
tdef    1 linear 00:00z20JUL2089 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方向のグリッド数：525
y方向のグリッド数：1721
z方向のグリッド数：1
タイムステップ数：1
変数数：11

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

データ形式は、気候予測データセット2022解説書のP2-56”モデル諸元”に記載があり、"32bit 実数バイナリ"とのこと。

Projectionの同定

座標系は「日本域気候予測データー格子点データ」と同じ、緯度60°N、緯度 30°N の両緯線を2標準緯線とし、中央子午線を80°Eとしたランベルト正角円錐図法を使います。

PROJ.4文字列にすると、

PROJ.4 String

+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パラメータの同定

2km格子なので、resolutionとrotationのパラメータは

Python

GT(2) = GT(4) = 0
GT(1) = 2000
GT(5) = -2000

になります。

upper-left pixelの座標はcnst2km.datから特定できます。

Python

# ライブラリの呼び出し
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/NHRCM02_SOUSEI/cnst2km.dat",dtype='>f4')
cnst_arr = cnst_vec.reshape(11,903527)

# 緯度と経度情報を配列に変換
FLAT_arr = cnst_arr[2,1:-1].reshape(1721,525)
FLON_arr = cnst_arr[3,1:-1].reshape(1721,525)

# 左上格子の緯度経度を取り出す
ul_lat = FLAT_arr[-1,0].astype(float)
ul_lon = FLON_arr[-1,0].astype(float)

print(ul_lat,ul_lon)
#→49.653743743896484 144.72854614257812

ランベルト正角円錐図法の座標系に変換し、右上角の座標を同定します。

Python

# ライブラリの呼び出し
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 - 1000
ul_y = ul_y + 1000

print(ul_x,ul_y)
#→4034196.596039401 7570619.639166898

結果、GeoTransFormのパラメータは

Python

[4034196.596039401,2000,0,7570619.639166898,0,-2000]

となります。

海陸比データをGeoTIFF化

同定したパラメータの検証のため、海陸比データをGeoTiff化します。

Python

# ライブラリの呼び出し
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/NHRCM02_SOUSEI/cnst2km.dat",dtype='>f4')
cnst_arr = cnst_vec.reshape(11,903527)
sl_arr = cnst_arr[1,1:-1].reshape(1721,525)

#GeoTiFFを作成
name = "cnst_2km_sl.gtiff"
output = gdal.GetDriverByName('GTiff').Create(name, 525, 1721 ,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([4034196.596039401,2000,0,7570619.639166898,0,-2000])


#海陸比データを設定
outband = output.GetRasterBand(1)
outband.WriteArray(sl_arr[::-1,:].astype("float64"))

outband.FlushCache()
output =None

QGISに読み込んでみる

うまく投影変換されていることが確認できます。

結果

気候予測データセット2022の「全球及び日本域気候予測データ」創生・統合プログラム2km格子NHRCM日本域気候予測データをGDALでラスタ化する際の投影パラメータは、

proj4文字列.

+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パラメータ.

[4034196.596039401,2000,0,7570619.639166898,0,-2000]

であると同定できました！

（個人）気候変化情報研究室では、できるだけ簡単に、将来、地域の気候がどのように変化するかを示すグラフ・地図を作成する方法を開発し、紹介しています。
興味のある方はぜひ読んでみてくださいね～。