はじめに
「全球及び日本域気候予測データ」創生・統合プログラム2km格子NHRCM日本域気候予測データをラスタ化する際の投影パラメータを同定しました。
「全球及び日本域気候予測データ」創生・統合プログラム2km格子NHRCM日本域気候予測データの入手とマップ化の記事はこちら。
GRIB2形式のファイルをそのままQGISに読み込むと、
5km格子NHRCM日本域気候予測データと同様に、一応読み込むことはできますが、このままでは日本列島の位置が合いません。
そこで、5㎞格子の時と同様に、投影パラメータを同定しました。
モデル格子点情報の入手
気候予測データセット(DS2022)から、創生・統合プログラム2km格子NHRCM日本域気候予測データの定数データを入手します。
ファイルは次の2つ
- cnst2km.dat
- cnst2km.ctl
ctlファイルの情報
GrADSで解析・表示するためのコントロールファイルは次のとおり。
***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=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のパラメータは
GT(2) = GT(4) = 0
GT(1) = 2000
GT(5) = -2000
になります。
upper-left pixelの座標はcnst2km.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/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
ランベルト正角円錐図法の座標系に変換し、右上角の座標を同定します。
# ライブラリの呼び出し
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のパラメータは
[4034196.596039401,2000,0,7570619.639166898,0,-2000]
となります。
海陸比データを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/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でラスタ化する際の投影パラメータは、
+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”
[4034196.596039401,2000,0,7570619.639166898,0,-2000]
であると同定できました!
(個人)気候変化情報研究室では、できるだけ簡単に、将来、地域の気候がどのように変化するかを示すグラフ・地図を作成する方法を開発し、紹介しています。
興味のある方はぜひ読んでみてくださいね~。