「全球及び日本域確率的気候予測データ（d4PDF シリーズ）」―全球データをラスタ化する際の投影パラメータを同定しました。
「全球及び日本域150年連続実験データ」の60km150年連続実験データも同じモデルを使用しており、このパラメータを使用できます。
モデル格子点情報の入手
気候予測データセット（DS2022）から、全球及び日本域確率的気候予測データ（d4PDF シリーズ）全球モデルの”実験共通：fixed"の"標高分布・海陸分布データ : TopogRatiol_gsmuv_TL319"データを入手します。
ファイルは次の2つ
- TopogRatiol_gsmuv_TL319.ctl
- TopogRatiol_gsmuv_TL319.gd
ctlファイルの情報
GrADSで解析・表示するためのコントロールファイルは次のとおり。
TopogRatiol_gsmuv_TL319.ctl
DSET ^TopogRatiol_gsmuv_TL319.gd
TITLE gsmuv TL959 topography data
OPTIONS BIG_ENDIAN
UNDEF -9.99E33
XDEF 640 LINEAR 0 0.5625
YDEF 320 LEVELS
-89.570 -89.013 -88.453 -87.892 -87.331
-86.769 -86.208 -85.647 -85.085 -84.523
-83.962 -83.400 -82.839 -82.277 -81.716
-81.154 -80.592 -80.031 -79.469 -78.908
-78.346 -77.784 -77.223 -76.661 -76.100
-75.538 -74.976 -74.415 -73.853 -73.291
-72.730 -72.168 -71.607 -71.045 -70.483
-69.922 -69.360 -68.799 -68.237 -67.675
-67.114 -66.552 -65.990 -65.429 -64.867
-64.306 -63.744 -63.182 -62.621 -62.059
-61.498 -60.936 -60.374 -59.813 -59.251
-58.689 -58.128 -57.566 -57.005 -56.443
-55.881 -55.320 -54.758 -54.196 -53.635
-53.073 -52.512 -51.950 -51.388 -50.827
-50.265 -49.704 -49.142 -48.580 -48.019
-47.457 -46.895 -46.334 -45.772 -45.211
-44.649 -44.087 -43.526 -42.964 -42.402
-41.841 -41.279 -40.718 -40.156 -39.594
-39.033 -38.471 -37.909 -37.348 -36.786
-36.225 -35.663 -35.101 -34.540 -33.978
-33.416 -32.855 -32.293 -31.732 -31.170
-30.608 -30.047 -29.485 -28.924 -28.362
-27.800 -27.239 -26.677 -26.115 -25.554
-24.992 -24.431 -23.869 -23.307 -22.746
-22.184 -21.622 -21.061 -20.499 -19.938
-19.376 -18.814 -18.253 -17.691 -17.129
-16.568 -16.006 -15.445 -14.883 -14.321
-13.760 -13.198 -12.636 -12.075 -11.513
-10.952 -10.390 -9.828 -9.267 -8.705
-8.144 -7.582 -7.020 -6.459 -5.897
-5.335 -4.774 -4.212 -3.651 -3.089
-2.527 -1.966 -1.404 -0.842 -0.281
0.281 0.842 1.404 1.966 2.527
3.089 3.651 4.212 4.774 5.335
5.897 6.459 7.020 7.582 8.144
8.705 9.267 9.828 10.390 10.952
11.513 12.075 12.636 13.198 13.760
14.321 14.883 15.445 16.006 16.568
17.129 17.691 18.253 18.814 19.376
19.938 20.499 21.061 21.622 22.184
22.746 23.307 23.869 24.431 24.992
25.554 26.115 26.677 27.239 27.800
28.362 28.924 29.485 30.047 30.608
31.170 31.732 32.293 32.855 33.416
33.978 34.540 35.101 35.663 36.225
36.786 37.348 37.909 38.471 39.033
39.594 40.156 40.718 41.279 41.841
42.402 42.964 43.526 44.087 44.649
45.211 45.772 46.334 46.895 47.457
48.019 48.580 49.142 49.704 50.265
50.827 51.388 51.950 52.512 53.073
53.635 54.196 54.758 55.320 55.881
56.443 57.005 57.566 58.128 58.689
59.251 59.813 60.374 60.936 61.498
62.059 62.621 63.182 63.744 64.306
64.867 65.429 65.990 66.552 67.114
67.675 68.237 68.799 69.360 69.922
70.483 71.045 71.607 72.168 72.730
73.291 73.853 74.415 74.976 75.538
76.100 76.661 77.223 77.784 78.346
78.908 79.469 80.031 80.592 81.154
81.716 82.277 82.839 83.400 83.962
84.523 85.085 85.647 86.208 86.769
87.331 87.892 88.453 89.013 89.570
ZDEF 1 LEVELS 0
TDEF 1 LINEAR jan2001 1mo
VARS 2
height 0 0 Topography [m]
ratiol 0 0 Ratio of land area [0-1]
ENDVARS
Projectionの同定
ctlファイルにpdef行はありません。XDEFとYDEFの値から、単純な緯度経度座標系（地理座標系）であると考えられます。また、気候モデルは球座標系で定式化されるので、楕円体は真球で、極半径＝赤道半径＝地球の平均半径（6371㎞）です。
PROJ.4文字列では、
PROJ.4 String
"+proj=longlat +a=6371000 +b=6371000 +over +no_defs"
となります。
GeoTransFormパラメータの同定
経度方向の格子配置は、XDEF行の記載から
XDEF 640 LINEAR 0 0.5625
子午線から東方向に0.5625°間隔で640格子（0～360°）が設定されています。
一般的な緯度経度座標系は経度が-180~180°の範囲を対象にしているため、これに合わせておかないと変換がうまくいかない場合があります。そこで、データ自体を-180~180°の範囲に合わせてGeoTransFormパラメータを設定することにします。
右端の座標は、320格子分マイナス（西）側に設定し、経度方向のパラメータは
Python
GT(0) = 0 - 0.5625*320 - 0.5625/2 = -180.28125
GT(1) = 0.5625
GT(2) = 0
となります。
一方で緯度方向については、YDEF行に各格子点の緯度が指定されており、ここから推定する必要があります。
Python
import os
import pandas as pd
os.chdir("/mnt/c/Users/kazuh/d4PDF_GCM")
substr_df = pd.read_fwf("./d4PDF_GCM/fixed/TopogRatiol_gsmuv_TL319.ctl",widths=[10],header=None,skip_blank_lines=False).fillna("#")
YDEF_ind = substr_df.loc[substr_df[0].str.startswith("YDEF")].index.values[0]
YDEF_nrow = substr_df.loc[substr_df[0].str.startswith("ZDEF")].index.values[0] - YDEF_ind - 1
FLAT_vec = pd.read_fwf("./d4PDF_GCM/fixed/TopogRatiol_gsmuv_TL319.ctl",header=None,skiprows=YDEF_ind+1,nrows=YDEF_nrow,colspecs=[(6,15),(15,23),(24,32),(33,41),(42,50)]).stack().values
格子間隔を計算し、グラフ化してみると
Python
import pandas as pd
import plotly.express as px
lat_df = pd.DataFrame({"lat":FLAT_vec})
lat_df["interval"] = lat_df["lat"].diff()
px.scatter(lat_df,x="lat",y="interval")
高緯度域で間隔が微妙に狭くなっているようです。ただ、その差は最大0.005度で、距離にすると約0.5km、低緯度域の格子間隔の1％程度。
GISのラスタデータとして扱うには、セルサイズを同じにする必要があり、低緯度域の格子間隔で統一します。
低緯度域の格子間隔は
Python
ind_from = lat_df.loc[lat_df["interval"] >= 0.561,"lat"].index[0]
ind_to = lat_df.loc[lat_df["interval"] >= 0.561,"lat"].index[-1]
lat_res = (lat_df.loc[ind_to,"lat"] - lat_df.loc[ind_from-1,"lat"]) / (ind_to - ind_from + 1)
print(lat_res)
#→0.5616063492063492
Upper側のy座標は半球分の格子数160を掛けることで求められます。
Python
ul_lat = lat_res * 160
print(89.85701587301588)
緯度方向のパラメータは
Python
GT(3) = 89.85701587301588
GT(4) = 0
GT(5) = -0.5616063492063492
結果、GeoTransFormのパラメータは
Python
[-180.28125,0.5625,0,89.85701587301588,0,-0.5616063492063492]
となります。
海陸比データを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_GCM")
#海陸比情報を読み込み、左右を入れ替える
topo_vec = np.fromfile("./d4PDF_GCM/fixed/TopogRatiol_gsmuv_TL319.gd",dtype='>f4')
topo_arr = topo_vec.reshape(2,320,640)
ratiol_arr = np.concatenate([topo_arr[1,::-1,320:],topo_arr[1,::-1,:320]],axis=1)
#GeoTiFFを作成
name = "Ratiol_gsmuv_TL319.gtiff"
output = gdal.GetDriverByName('GTiff').Create(name, 640, 320, 1, gdal.GDT_Float64)
#Projectionを設定
outRasterSRS = osr.SpatialReference()
outRasterSRS.ImportFromProj4(
"+proj=longlat +a=6371000 +b=6371000 +no_defs")
output.SetProjection(outRasterSRS.ExportToWkt())
#GeoTransformを設定
output.SetGeoTransform([-180.28125,0.5625,0,89.85701587301588,0,-0.5616063492063492])
#海陸比データを設定
outband = output.GetRasterBand(1)
outband.WriteArray(ratiol_arr.astype("float64"))
outband.FlushCache()
output =None
QGISに読み込んでみる
うまく投影されていることが確認できます。
結果
気候予測データセット2022の「全球及び日本域確率的気候予測データ（d4PDF シリーズ）」―日本域ダウンスケーリングデータをラスタ化する際の投影パラメータ
proj4文字列.
"+proj=longlat +a=6371000 +b=6371000 +over +no_defs"
GeoTransFormパラメータ.
[-180.28125,0.5625,0,89.85701587301588,0,-0.5616063492063492]
であると同定できました！
（個人）気候変化情報研究室では、できるだけ簡単に、将来、地域の気候がどのように変化するかを示すグラフ・地図を作成する方法を開発し、紹介しています。
興味のある方はぜひ読んでみてくださいね～。