モデル格子点情報の入手

気候予測データセット（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]

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

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