GDALで「全国5kmメッシュアンサンブル気候予測データ」の投影パラメータを同定する

## はじめに

先日、「全国5kmメッシュアンサンブル気候予測データ」が公開されました。

https://diasjp.net/infomation/20231205/

これは、地球温暖化に資するアンサンブル気候予測データベース(d4PDF)の現在気候（1950－2010 年）、産業革命時から世界平均で２度上昇（RCP8.5 シナリオで近未来 2040 年頃）及び４度上昇時（RCP8.5 シナリオで 21 世紀末 2090 年頃）の気候予測データのそれぞれ732年分について、地域気候モデルを用いて5㎞にダウンスケーリングしたもので、大雨や大雪、猛暑などの極端気象が温暖化が進むに伴いどのように変わるのか？その頻度や強度の変化を評価することができます。

https://www.miroc-gcm.jp/d4PDF/

データ形式は、GIS界隈に優しいnetCDFで公開されています。が、そのままでは地図に重ねてデータを切り出すことは難しそう。
ということで、GDALを使って投影パラメータを同定します。

> Q. データの形式、出力変数を教えてください。
A. netCDF形式(***.nc)で、各変数ごとに時別値で提供しております。また、日本時間0時を日界とした日別値も提供しております。日別値は変数をまとめています。詳しくは利用の手引きをご覧ください。なお、netCDFの格子点情報はモデルの格子数となっており、緯度経度ではありません。対応する緯度経度は、cnst.ncからご取得ください。
>(d4PDF_5kmDDS_JP よくある質問（FAQ）より)

## モデル格子点情報の入手
DIASデータ俯瞰検索システムの全国5kmメッシュアンサンブル気候予測データページから、緯度経度、地形、海陸分布データ(cnst.nc)を入手します。

## Projectionの同定

利用の手引きの表1に、
地図投影法は

>ランベルト正角円錐図法
（標準緯度：北緯 30 度，60 度，標準経度：東経 80 度）

とありますので、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.ncから特定できます。

```python:Python
# ライブラリのインポート
import os
import numpy as np
import xarray as xr

# cnst.ncの読み込み
os.chdir("/mnt/c/Users/hoge/d4PDF_5km")
ds = xr.open_dataset("cnst.nc")

#左上隅グリッドの緯度経度を取り出す
ul_ds = ds.isel(lat=-1,lon=0)
ul_lat = ul_ds["flat"].data[0][0].astype(float)
ul_lon = ul_ds["flon"].data[0][0].astype(float)

print(ul_lat,ul_lon)
# 57.08604431152344 137.833984375
```

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

```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)
# 3153257.170699811 7844357.000364543
```
結果、GeoTransFormのパラメータは

```Python:Python
[3153257.170699811,5000,0,7844357.000364543,0,-5000]
```

になります。

## 海陸比データをGeoTIFF化

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

os.chdir("/mnt/d/Users/hoge/d4PDF_5km")

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

#海陸比データを設定
outband = output.GetRasterBand(1)
outband.WriteArray(ds["sl"].data[0,0,::-1,:].astype(float))

outband.FlushCache()
output =None
```

## QGISに読み込んでみる

![スクリーンショット 2023-12-09 172850.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3511025/4e379cfe-58c6-7676-2cab-e42b68a460a7.jpeg)

![スクリーンショット 2023-12-09 173051.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3511025/2dfcb0ca-fa3a-96f5-ab9b-3f7f146d601a.jpeg)

![スクリーンショット 2023-12-09 173200.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3511025/6c699c08-fbcf-d056-2f77-41a84358ddf8.jpeg)

![スクリーンショット 2023-12-09 173250.jpg](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/3511025/bb0032f0-4e7f-c1c2-3df2-5c76deaa98bd.jpeg)

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

## 結果

「全国5kmメッシュアンサンブル気候予測データ」を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パラメータ.
[3153257.170699811,5000,0,7844357.000364543,0,-5000]
```

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

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

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