2018/01/09

地形モデル

遠地津波計算をするためにGEBCOのデータを変換する手順は以前書きましたので,今回は近地津波の計算をするための直交座標系での地形データの作成手順を書いておきます.

元データは内閣府の南海トラフの巨大地震モデル検討会で作成された地形データがG空間情報センターからダウンロード出来ます.無償にはなっていますが,ダウンロードするためにはアカウントを作成する必要があります.地形データ以外の粗度や堤防,初期水位データ等もダウンロード出来ます..

例として,JAGURS で使える形式に変換してみます.
手順は大きく分けで3つ.(1)まずはデータをダウンロードして,(2)まともなGISフォーマットに変換して,(3)最後に必要なフォーマットへの変換をします.
地形データは,ここからダウンロードできます.利用規約はもちろん確認が必要ですが,「計算入力データ(地形,粗度,堤防データ)の説明」を確認して下さい.ここに何系がどの領域に対応しているか書かれています.例えば四国を計算する場合は第04系になるのでその地形データをダウンロードします.ダウンロードしたデータにはメタ情報が何も書かれていないので,別途「計算範囲設定」を確認する必要があります.下記の ncols, nrows, xllcorner, yllcorner の値は,エクセルのメッシュ個数(X方向,Y方向),南西端の位置(X座標,Y座標)に対応します.

データは 10f8.2 で記述してあるということなので,Fortranで読み込むのは簡単ですが, f8.2 の場合,標高(陸域)が1000mを超える場合に,符号も含めると -1000.00 となり8桁になってしまうため,前の数字とスペースなしで連結されてしまいます.Fortranなら読めるでしょって言ってしまえばそれまでなのですが,一般的なアスキーフォーマットと考えた場合,スペースで区切られていないと使えません.というわけで, 10f8.2read して, 10f9.2write します.さらに,アスキーフォーマットとしてGISデータになるようにヘッダを付与します.プログラム(reformat.f)も掲載しておきます.
     ncols 720
     nrows 540
     xllcorner -730200
     yllcorner -857700
     cellsize 2430
     NODATA_value -9999

ヘッダの中身は,さほど難しくなく,ncolsはX方向のメッシュ数,nrowsはY方向のメッシュ数,xllcorner は南西端のX座標,yllcorner は南西端のY座標,cellsize は解像度(上記例は 2430m),NODATA_value は不要ですがお決まりの -9999 を設定しています.

reformat.f:
      PROGRAM reformat
      real*8 data(10)
      character(len=128) :: fname
      read(*,'(a)') fname
      open(10,file=fname)
      do
         read(10,'(10f8.2)',end=100) data
         write(*,'(10f9.2)') data
      end do
 100  END

作成したアスキーファイル(bathy.SD01.asc)をJAGURS用のcf形式に変換します.
% grdreformat bathy.SD01.asc=gd bathy.SD01.grd=cf
これで,ネスティングしない場合は計算できるようになりました.ネスティングする場合は,メッシュの細かい調整が必要なので後日掲載します(一応,手元では動くデータは作れましたがこれで良いか確認します).

日本海講側は,中央防災会議のものが公開されています(古いです).