「日本域気候予測データ―格子点データ」5km版の投影パラメータを同定する

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

「気候予測データセット2022」はこちら

https://diasjp.net/ds2022/dataset/ds02.html

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

https://cci-labo.hateblo.jp/entry/2023/07/11/083451

格子点データから作ったGeoTiffをQGISに正しく読み込むためには、GDALの”Projection”と、”GeoTransform"のパラメータを正しく設定せねばならぬ。
ということで、投影パラメータを同定しました。

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

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

ファイルは次の3つ。

+ cnst.dat
+ cnst.ctl
+ cnst.csv

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

https://cci-labo.hateblo.jp/entry/2023/06/06/210638

## ctlファイルのPDEF情報

日本域気候予測データのモデル格子点データには、GrADSで解析・表示するためのコントロールファイル（cnst.ctl）がついています。

```ctl:cnst.ctl
***output grads ctl file for imt_sf***
dset ^cnst.dat
title surface result
OPTIONS big_endian 
undef -999.9
pdef 467 744 LCCR 44 144 214 616 30 60 80 5000 5000
xdef  920 linear 114  0.05
ydef  800 linear 16   0.05
zdef   1 linear 1 1
tdef   1 linear 01z01SEP1980    1hr
vars   4
FLAT	1	99	deg
FLON	1	99	deg
ZS	1	99	m
SL	1	99	sl
endvars
***end of ctl file***
```

このうち、このうち"pdef"行にバイナリデータのグリッド設定が記載されています。
その文法は以下を参照ください。

http://cola.gmu.edu/grads/gadoc/pdef.htm

PDEF行を読み解くと、
```pdef情報. 
pdef 467 744 LCCR 44 144 214 616 30 60 80 5000 5000
↓
isize:467	The size of the native grid in the x direction
jsize:744	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:214	i of ref point
jref:616	j of ref point
Struelat:30	S true lat
Ntruelat:60	N true lat
slon:80	standard longitude
dx:5000	grid X increment in meters
dy:5000	grid Y increment in meters
```
となります。

## Projectionの同定

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

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

これを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
```

になります。

PROJ.4文字列についてはこれを参考にさせていただきました。

https://pygis.io/docs/d_understand_crs_codes.html

## GeoTransFormパラメータの同定

GDALのラスターデータには、ピクセル/ライン座標を地理参照空間に正しく投影するためのGeoTransFormが必要です。

https://gdal.org/tutorials/geotransforms_tut.html#geotransform-tutorial

計算式は以下

```GeoTranForm式．
Xgeo = GT(0) + Xpixel*GT(1) + Yline*GT(2)
Ygeo = GT(3) + Xpixel*GT(4) + Yline*GT(5)

GT(0) x-coordinate of the upper-left corner of the upper-left pixel.
GT(1) w-e pixel resolution / pixel width.
GT(2) row rotation (typically zero).
GT(3) y-coordinate of the upper-left corner of the upper-left pixel.
GT(4) column rotation (typically zero).
GT(5) n-s pixel resolution / pixel height (negative value for a north-up image).
```
ctlファイルのPDEF情報から、
```python:Python
GT(2) = GT(4) = 0
GT(1) = 5000
GT(5) = -5000
```
になります。

あとはGT(0)とGT(3)ですが、cnst.csvからupper-left pixelの座標が特定できます。

```python:Python
import os
import pandas as pd

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

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

print(ul_lat,ul_lon)
# →55.145286560058594 139.4412841796875
```

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

```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=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 - 2500
ul_y = ul_y + 2500

print(ul_x,ul_y)
#→3367936.9987040577 7752086.241150136
```

結果、GeoTransFormのパラメータは
```python:Python
[3367936.9987040577,5000,0,7752086.241150136,0,-5000]
```
となります。

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

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

#GeoTiFFを作成
name = "cnst_5km_sl.gtiff"
output = gdal.GetDriverByName('GTiff').Create(name, 467, 744,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([3367936.9987040577,5000,0,7752086.241150136,0,-5000])

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

outband.FlushCache()
output =None
```

## QGISに読み込んでみる

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

![QGIS_sl_jpn.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3511025/378ec1cb-005a-1dea-5db3-fc8679a437e1.jpeg)
![QGIS_sl_hok.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3511025/97bdecb2-6437-e35a-846a-fa43a1db5958.jpeg)
![QGIS_sl_tk.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3511025/cebea0c6-2b1b-ca50-0c95-abc61b345a53.jpeg)
![QGIS_sl_osk.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3511025/7d9a8c72-1acc-28d4-e8ab-3b7fdbd3a684.jpeg)
![QGIS_sl_kyuu.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3511025/38dd7322-48ce-8d5b-75b1-b290ee94d24f.jpeg)

## 結果

気候予測データセット2022の「日本域気候予測データ5㎞版」格子点データを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パラメータ.
[3367936.9987040577,5000,0,7752086.241150136,0,-5000]
```

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

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

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