「全球及び日本域気候予測データ」創生・統合プログラム5km格子NHRCM日本域気候予測データの投影パラメータを同定する

## はじめに

「全球及び日本域気候予測データ」創生・統合プログラム5km格子NHRCM日本域気候予測データをラスタ化する際の投影パラメータを同定しました。

データの入手とマップ化の記事はこちら。

https://cci-labo.hateblo.jp/entry/2023/07/26/084301

データ自体はGRIB2形式で提供されており、うまく設定すればそのままQGISで表示できるのですが、このデータには投影パラメータが設定されていないようで、そのままQGISに読み込むと、

![QGIS_raw_data.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3511025/1dad41e0-4e36-1982-ee1a-51af4d56a8b6.jpeg)

一応読み込むことはできますが、このままでは日本列島の位置が合いません。

そこで、「日本域気候予測データ－格子点データ」5㎞版の時と同様に、投影パラメータを同定しました。

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

気候予測データセット（DS2022）から、創生・統合プログラム5km格子NHRCM日本域気候予測データの定数データを入手します。

ファイルは次の2つ

+ cnst.dat
+ cnst.ctl

## ctlファイルの情報

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


```ctl:cnst.ctl
***output grads ctl file for imt_cf***
dset cnst.dat
title constatnt data
OPTIONS big_endian sequential 
undef -999.9
xdef  527 linear 1 1
ydef  804 linear 1 1
zdef    1 linear 1 1
tdef    1 linear 00:00z20JUL1980 1mn
vars 11
ZS     1 99 ZS
SL     1 99 SL
FLAT   1 99 FLAT
FLON   1 99 FLON
KIND   1 99 KIND
ROUGH  1 99 ROUGH
WET    1 99 WET
ALBED  1 99 ALBED
FKTG   1 99 FKTG
ROCTG  1 99 ROCTG
URBAN  1 99 URBAN
endvars
***end of ctl file***
```

ctlの記載から、このファイルは
+ x方向のグリッド数：527
+ y方向のグリッド数：804
+ z方向のグリッド数：1
+ タイムステップ数：1
+ 変数数：11

であることが読み取れます。

データ形式は、気候予測データセット2022解説書のP2-56”モデル諸元”に記載があり、"32bit 実数バイナリ"とのこと。

## Projectionの同定

座標系は「日本域気候予測データー格子点データ」と同じ、緯度60°N、緯度 30°N の両緯線を2標準緯線とし、中央子午線を80°Eとしたランベルト正角円錐図法を使います。

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
```

になります。

## GeoTransFormパラメータの同定

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

```python:Python
GT(2) = GT(4) = 0
GT(1) = 5000
GT(5) = -5000
```
になります。

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

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

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

# 4バイト実数として読み込み、変数毎に分割
cnst_vec = np.fromfile("./dias/data/NHRCM05_SOUSEI/cnst.dat",dtype='>f4')
cnst_arr = cnst_vec.reshape(11,423710)

# 緯度と経度情報を配列に変換
FLAT_arr = cnst_arr[2,1:-1].reshape(804,527)
FLON_arr = cnst_arr[3,1:-1].reshape(804,527)

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

print(ul_lat,ul_lon)
#→57.08396 139.28857
```

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

```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)
```

結果、GeoTransFormのパラメータは

```Python:Python
[3217950.058222905,5000,0,7902104.780171389,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/MHRCM")

#海陸比情報を読み込む
cnst_vec = np.fromfile("./dias/data/NHRCM05_SOUSEI/cnst.dat",dtype='>f4')
cnst_arr = cnst_vec.reshape(11,423710)
sl_arr = cnst_arr[1,1:-1].reshape(804,527)

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

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

outband.FlushCache()
output =None
```

## QGISに読み込んでみる

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

![QGIS_SL_jpn.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3511025/40feb4bc-7113-0a96-d1bd-70acf344488d.jpeg)

![QGIS_SL_hok.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3511025/a63b9024-a43c-2c76-f132-9d8119cbb4de.jpeg)

![QGIS_SLkanto.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3511025/8de55cdd-89d4-aae8-7253-265c279df9dc.jpeg)

![QGIS_SL_kansai.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3511025/dccaa67b-efab-4eaf-f25c-f55247f7e872.jpeg)

![QGIS_SL_kyusyu.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3511025/e8b3ff56-e626-ffad-6852-9f27c292006d.jpeg)

## 結果

気候予測データセット2022の「全球及び日本域気候予測データ」創生・統合プログラム5km格子NHRCM日本域気候予測データを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パラメータ.
[3217950.058222905,5000,0,7902104.780171389,0,-5000]
```

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

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

興味のある方はぜひ読んでみてくださいね～。

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