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

## はじめに

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

前回に続き、「日本域気候予測データ―格子点データ」2㎞版の投影パラメータを同定します。

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

https://cci-labo.hateblo.jp/entry/2023/07/13/092644

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

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

ファイルは次の3つ。

+ cnst_2km.dat
+ cnst_2km.ctl
+ cnst_2km.csv

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

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

## ctlファイルのPDEF情報

```ctl:cnst_2km.ctl
DSET ^cnst_2km.dat
OPTIONS big_endian
title constant data
undef -999.9
pdef 485 1681 LCCR 44 144 178 1474 30 60 80 2000 2000
xdef  2300 linear 114  0.02
ydef  2000 linear 16   0.02
zdef    1 linear 1 1
tdef 1 linear 00z01SEP1980 1hr
vars   4
FLAT  1 99 deg
FLON  1 99 deg
ZS    1 99 m
SL    1 99 sl
endvars
```

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

## Projectionの同定
座標系は5㎞版と同じで、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パラメータの同定

ctlファイルのPDEF情報から、
```python:Python
GT(2) = GT(4) = 0
GT(1) = 2000
GT(5) = -2000
```
になります。

GT(0)とGT(3)については、cnst_2km.csvからupper-left pixelの座標を特定できます。

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

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

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

print(ul_lat,ul_lon)
# →49.128177642822266 144.71006774902344
```

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

```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)
#　→4074196.5514666005 7530619.71549598
```

結果、GeoTransFormのパラメータは
```python:Python
[4074196.5514666005,2000,0,7530619.71549598,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/JWP9")

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

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

#海陸比データを設定
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/d5ac591f-e77f-f89b-dd99-d521f901f52b.jpeg)
![QGIS_sl_hok.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3511025/4ffe4aff-1634-678e-6edd-79e21d5b7804.jpeg)
![QGIS_sl_tk.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3511025/c93142bc-4ed6-f6da-b762-6b83e49994bd.jpeg)
![QGIS_sl_os.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3511025/8f34a32c-1dc9-7a11-ce21-373fa39a064c.jpeg)
![QGIS_sl_fuk.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3511025/1d20f0b5-a8ed-e109-8938-486eba7311fa.jpeg)

## 結果

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


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

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

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