はじめに

気候予測データセット2022の「日本域気候予測データ5㎞版」格子点データをラスタ化する際の投影パラメータを同定しました。

「気候予測データセット2022」はこちら

「日本域気候予測データ－格子点データ」5㎞版の入手とマップ化の記事はこちら。

格子点データから作ったGeoTiffをQGISに正しく読み込むためには、GDALの”Projection”と、”GeoTransform"のパラメータを正しく設定せねばならぬ。
ということで、投影パラメータを同定しました。

モデル格子情報の入手

気候予測データセット（DS2022）から、日本域気候予測データのモデル格子情報を入手します。

ファイルは次の3つ。

cnst.dat
cnst.ctl
cnst.csv

データの入手方法は以下を参考にしてください。

ctlファイルのPDEF情報

日本域気候予測データのモデル格子点データには、GrADSで解析・表示するためのコントロールファイル（cnst.ctl）がついています。

cnst.ctl

***output grads ctl file for imt_sf***
dset ^cnst.dat
title surface result
OPTIONS big_endian 
undef -999.9
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   1 linear 01z01SEP1980    1hr
vars   4
FLAT	1	99	deg
FLON	1	99	deg
ZS	1	99	m
SL	1	99	sl
endvars
***end of ctl file***

このうち、このうち"pdef"行にバイナリデータのグリッド設定が記載されています。
その文法は以下を参照ください。

PDEF行を読み解くと、

pdef情報.

pdef 467 744 LCCR 44 144 214 616 30 60 80 5000 5000
↓
isize:467	The size of the native grid in the x direction
jsize:744	The size of the native grid in the y direction
projection:LCCR the Lambert Conformal projection
latref:44	reference latitude
lonref:144	reference longitude (in degrees, E is positive, W is negative)
iref:214	i of ref point
jref:616	j of ref point
Struelat:30	S true lat
Ntruelat:60	N true lat
slon:80	standard longitude
dx:5000	grid X increment in meters
dy:5000	grid Y increment in meters

となります。

Projectionの同定

pdef情報から、緯度60°N、緯度 30°N の両緯線を2標準緯線とし、中央子午線を80°Eとしたランベルト正角円錐図法を使っていることがわかります。

気候モデルは球座標系で定式化されるので、楕円体は真球で、極半径＝赤道半径＝地球の平均半径（6371㎞）です。

これを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

になります。

PROJ.4文字列についてはこれを参考にさせていただきました。

GeoTransFormパラメータの同定

GDALのラスターデータには、ピクセル/ライン座標を地理参照空間に正しく投影するためのGeoTransFormが必要です。

計算式は以下

Xgeo = GT(0) + Xpixel*GT(1) + Yline*GT(2)
Ygeo = GT(3) + Xpixel*GT(4) + Yline*GT(5)

GT(0) x-coordinate of the upper-left corner of the upper-left pixel.
GT(1) w-e pixel resolution / pixel width.
GT(2) row rotation (typically zero).
GT(3) y-coordinate of the upper-left corner of the upper-left pixel.
GT(4) column rotation (typically zero).
GT(5) n-s pixel resolution / pixel height (negative value for a north-up image).

ctlファイルのPDEF情報から、

Python

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

になります。

あとはGT(0)とGT(3)ですが、cnst.csvからupper-left pixelの座標が特定できます。

Python

import os
import pandas as pd

os.chdir("/mnt/c/Users/hoge/JWP9")

# cnst.csvを読み込む
cnst_df = pd.read_csv("./dias/data/jmagwp/gwp9/meta/cnst.csv")
cnst_df.columns = ["x","y","flat","flon","zs","sl"]

# upper-left pixelの座標を特定
ul_lat,ul_lon = cnst_df.loc[(cnst_df["x"]==1)&(cnst_df["y"]==744),["flat","flon"]].values[0]

print(ul_lat,ul_lon)
# →55.145286560058594 139.4412841796875

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

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 - 2500
ul_y = ul_y + 2500

print(ul_x,ul_y)
#→3367936.9987040577 7752086.241150136

結果、GeoTransFormのパラメータは

Python

[3367936.9987040577,5000,0,7752086.241150136,0,-5000]

となります。

海陸比データを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/JWP9")

#モデル格子点情報を読み込む
cnst_dat = np.fromfile("./dias/data/jmagwp/gwp9/meta/cnst.dat",dtype='>f4')
var_arr =  cnst_dat.reshape(4,347448)
sl_arr = var_arr[3].reshape(744,467)[::-1,:]

#GeoTiFFを作成
name = "cnst_5km_sl.gtiff"
output = gdal.GetDriverByName('GTiff').Create(name, 467, 744,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([3367936.9987040577,5000,0,7752086.241150136,0,-5000])

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

outband.FlushCache()
output =None

QGISに読み込んでみる

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

結果

気候予測データセット2022の「日本域気候予測データ5㎞版」格子点データを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パラメータ.

[3367936.9987040577,5000,0,7752086.241150136,0,-5000]

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

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