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

はじめに

先日、「全国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.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

GT(2) = GT(4) = 0
GT(1) = 5000
GT(5) = -5000

になります。

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

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

# ライブラリの呼び出し
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

[3153257.170699811,5000,0,7844357.000364543,0,-5000]

になります。

海陸比データをGeoTIFF化

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に読み込んでみる

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

結果

「全国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]

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

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