はじめに
先日、「全国5kmメッシュアンサンブル気候予測データ」が公開されました。
これは、地球温暖化に資するアンサンブル気候予測データベース(d4PDF)の現在気候（1950－2010 年）、産業革命時から世界平均で２度上昇（RCP8.5 シナリオで近未来 2040 年頃）及び４度上昇時（RCP8.5 シナリオで 21 世紀末 2090 年頃）の気候予測データのそれぞれ732年分について、地域気候モデルを用いて5㎞にダウンスケーリングしたもので、大雨や大雪、猛暑などの極端気象が温暖化が進むに伴いどのように変わるのか？その頻度や強度の変化を評価することができます。
データ形式は、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=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のパラメータは
GT(2) = GT(4) = 0
GT(1) = 5000
GT(5) = -5000
になります。
upper-left pixelの座標はcnst.ncから特定できます。
# ライブラリのインポート
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
ランベルト正角円錐図法の座標系に変換し、右上角の座標を同定します。
# ライブラリの呼び出し
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のパラメータは
[3153257.170699811,5000,0,7844357.000364543,0,-5000]
になります。
海陸比データをGeoTIFF化
# ライブラリの呼び出し
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に読み込んでみる
うまく投影変換されていることが確認できます。
結果
「全国5kmメッシュアンサンブル気候予測データ」をGDALでラスタ化する際の投影パラメータは、
+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”
[3153257.170699811,5000,0,7844357.000364543,0,-5000]
であると同定できました！
（個人）気候変化情報研究室では、できるだけ簡単に、将来、地域の気候がどのように変化するかを示すグラフ・地図を作成する方法を開発し、紹介しています。
興味のある方はぜひ読んでみてくださいね～