はじめに

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

前回に続き、「日本域気候予測データ―格子点データ」2㎞版の投影パラメータを同定します。

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

モデル格子情報の入手

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

ファイルは次の3つ。

cnst_2km.dat
cnst_2km.ctl
cnst_2km.csv

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

ctlファイルのPDEF情報

cnst_2km.ctl

DSET ^cnst_2km.dat
OPTIONS big_endian
title constant data
undef -999.9
pdef 485 1681 LCCR 44 144 178 1474 30 60 80 2000 2000
xdef  2300 linear 114  0.02
ydef  2000 linear 16   0.02
zdef    1 linear 1 1
tdef 1 linear 00z01SEP1980 1hr
vars   4
FLAT  1 99 deg
FLON  1 99 deg
ZS    1 99 m
SL    1 99 sl
endvars

PDEF行を読み解くと、

pdef情報.

pdef 485 1681 LCCR 44 144 178 1474 30 60 80 2000 2000
↓
isize:485	The size of the native grid in the x direction
jsize:1681	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:178	i of ref point
jref:1474	j of ref point
Struelat:30	S true lat
Ntruelat:60	N true lat
slon:80	standard longitude
dx:2000	grid X increment in meters
dy:2000	grid Y increment in meters

となります。

Projectionの同定

座標系は5㎞版と同じで、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パラメータの同定

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

Python

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

になります。

GT(0)とGT(3)については、cnst_2km.csvからupper-left pixelの座標を特定できます。

Python

import os
import pandas as pd

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

# cnst_2km.csvを読み込む
cnst_df = pd.read_csv("./dias/data/jmagwp/gwp9/meta_2km/cnst_2km.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"]==1681),["flat","flon"]].values[0]

print(ul_lat,ul_lon)
# →49.128177642822266 144.71006774902344

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

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)
#　→4074196.5514666005 7530619.71549598

結果、GeoTransFormのパラメータは

Python

[4074196.5514666005,2000,0,7530619.71549598,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/JWP9")

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

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

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

outband.FlushCache()
output =None

QGISに読み込んでみる

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

結果

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

[4074196.5514666005,2000,0,7530619.71549598,0,-2000]

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

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