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

## はじめに

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

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

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

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

https://cci-labo.hateblo.jp/entry/2023/08/14/084211

## モデル格子点情報の入手

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

ファイルは次の2つ

+ topo_essp20.dat
+ topo_essp20.ctl

## ctlファイルの情報

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

```ctl: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行を読み解くと、

```ctl: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文字列にすると、

```ctl: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:Python
GT(2) = GT(4) = 0
GT(1) = 20000
GT(5) = -20000
```
になります。

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

```python: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: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:Python
[-1930000.1287495417,20000,0,5833346.794147097,0,-20000]
```

となります。

## 海陸比データをGeoTIFF化

同定したパラメータの検証のため、海陸比データをGeoTiff化します。

```python: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に読み込んでみる

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

![QGIS_LCC_jpn.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3511025/98983af1-241d-a7e6-632a-af8c9cd6d815.jpeg)

![QGIS_LCC_hok.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3511025/022c9bda-2b1e-18b0-bc27-13229a0e65a5.jpeg)

![QGIS_LCC_tohoku.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3511025/41f1acd7-b190-6db3-f271-f7fb35f34df8.jpeg)

![QGIS_LCC_cyubu.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3511025/019decaf-0c71-e37b-6ede-3cce51357fd7.jpeg)

![QGIS_LCC_shikokukyusyu.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3511025/e1764dfa-e4c3-9b1d-e0fc-84e2ada8f073.jpeg)

## 結果

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

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

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

https://cci-labo.hateblo.jp/