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


## はじめに

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

「全球及び日本域気候予測データ」創生・統合プログラム2km格子NHRCM日本域気候予測データの入手とマップ化の記事はこちら。

https://cci-labo.hateblo.jp/entry/2023/08/01/102338

GRIB2形式のファイルをそのままQGISに読み込むと、

![QGIS_TSFC_latlon.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3511025/a3169d55-742e-9cc5-d3a3-6042ec479c3a.jpeg)

5km格子NHRCM日本域気候予測データと同様に、一応読み込むことはできますが、このままでは日本列島の位置が合いません。

そこで、5㎞格子の時と同様に、投影パラメータを同定しました。

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

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

ファイルは次の2つ

+ cnst2km.dat
+ cnst2km.ctl

## ctlファイルの情報

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

```ctl:cnst2km.ctl
***output grads ctl file for imt_cf***
dset cnst2km.dat
title constatnt data
OPTIONS big_endian sequential 
undef -999.9
xdef  525 linear 1 1
ydef 1721 linear 1 1
zdef    1 linear 1 1
tdef    1 linear 00:00z20JUL2089 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方向のグリッド数：525
y方向のグリッド数：1721
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パラメータの同定

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

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

upper-left pixelの座標はcnst2km.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/NHRCM02_SOUSEI/cnst2km.dat",dtype='>f4')
cnst_arr = cnst_vec.reshape(11,903527)

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

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

print(ul_lat,ul_lon)
#→49.653743743896484 144.72854614257812
```

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

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

print(ul_x,ul_y)
#→4034196.596039401 7570619.639166898
```
結果、GeoTransFormのパラメータは

```Python:Python
[4034196.596039401,2000,0,7570619.639166898,0,-2000]
```

となります。

## 海陸比データを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/NHRCM02_SOUSEI/cnst2km.dat",dtype='>f4')
cnst_arr = cnst_vec.reshape(11,903527)
sl_arr = cnst_arr[1,1:-1].reshape(1721,525)

#GeoTiFFを作成
name = "cnst_2km_sl.gtiff"
output = gdal.GetDriverByName('GTiff').Create(name, 525, 1721 ,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([4034196.596039401,2000,0,7570619.639166898,0,-2000])


#海陸比データを設定
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/5426627a-bf7f-3539-4a45-b3c50a677305.jpeg)

![QGIS_SL_hok.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3511025/afbd63b9-5c3d-7f25-8175-a7703cd0c58a.jpeg)

![QGIS_SL_kanto.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3511025/873e0bca-18a6-de89-b24c-f13e70c99a1a.jpeg)

![QGIS_SL_kansai.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3511025/ba8602bd-b2e5-1731-75d3-d14328e98e9c.jpeg)

![QGIS_SL_kyusyu.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3511025/91c84302-76cc-cfb1-b81a-eac2f5e36e8d.jpeg)

## 結果

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

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

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

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