はじめに

「全球及び日本域確率的気候予測データ（d4PDF シリーズ）」―日本域ダウンスケーリングデータをラスタ化する際の投影パラメータを同定しました。

「全球及び日本域150年連続実験データ」の20km150年連続実験データも同じモデルを使用しており、このパラメータを使用できます。

2.2. モデル諸元
本データセットの作成のために使用された NHRCM の諸元は VI.2.2 日本域確率的気候予測データ（d4PDF シリーズ、20km 格子 NHRCM）と同様である。　(気候予測データセット解説書第２章 P2-162より引用)

「全球及び日本域150年連続実験データ」ー統合プログラム領域20km150年連続実験データの入手とマップ化の記事はこちら。

モデル格子点情報の入手

気候予測データセット（DS2022）から、全球及び日本域確率的気候予測データ（d4PDF シリーズ）日本域ダウンスケーリングの”実験共通：fixed"の"標高分布・海陸分布データ : topo_essp20"データを入手します。

ファイルは次の2つ

topo_essp20.dat
topo_essp20.ctl

ctlファイルの情報

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

topo_essp20.ctl

dset ^topo_essp20.dat
pdef 191 155 LCCR 35 135 97 77 30 60 135 20000 20000
options  big_endian sequential
title  LAT LON ZS SL
undef -999.999
xdef 321 linear  105.  0.2
ydef 205 linear   15.  0.2
zdef 1 linear 1 1
tdef 1 linear 01:00z20JUL1950 60mn
vars 4
flat 0 99 latitude
flon 0 99 longitude
zs   0 99 height
sl   0 99 land fraction over grid
endvars

PDEF行を読み解くと、

pdef情報.

pdef 191 155 LCCR 35 135 97 77 30 60 135 20000 20000
↓
isize:191	The size of the native grid in the x direction
jsize:155	The size of the native grid in the y direction
projection:LCCR the Lambert Conformal projection
latref:35	reference latitude
lonref:135	reference longitude (in degrees, E is positive, W is negative)
iref:97	i of ref point
jref:77	j of ref point
Struelat:30	S true lat
Ntruelat:60	N true lat
slon:135	standard longitude
dx:20000	grid X increment in meters
dy:20000	grid Y increment in meters

となります。

Projectionの同定

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

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

これをPROJ.4文字列にすると、

PROJ.4 String

"+proj=lcc +lat_1=30.0 +lat_2=60.0 +lon_0=135.0 +lat_0=0 +x_0=0 +y_0=0 +a=6371000 +b=6371000 +units=m +no_defs"

となります。

GeoTransFormパラメータの同定

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

Python

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

になります。

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

Python

# ライブラリの呼び出し
import os
import numpy as np
from osgeo import gdal, gdalconst
from osgeo import osr

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

# 4バイト実数として読み込み、変数毎に分割
topo_vec = np.fromfile("./d4PDF_RCM/fixed/topo_essp20.dat",dtype='>f4')
topo_arr = topo_vec.reshape(4,29607)

# 緯度と経度情報を配列に変換
FLAT_arr = topo_arr[0,1:-1].reshape(155,191)
FLON_arr = topo_arr[1,1:-1].reshape(155,191)

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

print(ul_lat,ul_lon)
#→46.48844528198242 108.55301666259766

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

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=135. +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 - 10000
ul_y = ul_y + 10000

print(ul_x,ul_y)
#→-1930000.1287495417 5833346.794147097

結果、GeoTransFormのパラメータは

Python

[-1930000.1287495417,20000,0,5833346.794147097,0,-20000]

となります。

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

#海陸比情報を読み込む
topo_vec = np.fromfile("./d4PDF_RCM/fixed/topo_essp20.dat",dtype='>f4')
topo_arr = topo_vec.reshape(4,29607)
SL_arr = topo_arr[3,1:-1].reshape(155,191)

#GeoTiFFを作成
name = "topo_essp20_sl.gtiff"
output = gdal.GetDriverByName('GTiff').Create(name, 191, 155 ,1,gdal.GDT_Float64)

#Projectionを設定
outRasterSRS = osr.SpatialReference()
outRasterSRS.ImportFromProj4(
"+proj=lcc +lat_1=30.0 +lat_2=60.0 +lon_0=135. +lat_0=0 +x_0=0 +y_0=0 +a=6371000 +b=6371000 +units=m +no_defs")
output.SetProjection(outRasterSRS.ExportToWkt())

#GeoTransformを設定
output.SetGeoTransform([-1930000.1287495417,20000,0,5833346.794147097,0,-20000])

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

outband.FlushCache()
output =None

QGISに読み込んでみる

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

結果

気候予測データセット2022の「全球及び日本域確率的気候予測データ（d4PDF シリーズ）」―日本域ダウンスケーリングデータをラスタ化する際の投影パラメータ

proj4文字列.

+proj=lcc +lat_1=30.0 +lat_2=60.0 +lon_0=135. +lat_0=0 +x_0=0 +y_0=0 +a=6371000 +b=6371000 +units=m +no_defs”

GeoTransFormパラメータ.

[-1930000.1287495417,20000,0,5833346.794147097,0,-20000]

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

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