内閣府 中央防災会議が作った海底地形データを読めるようにしてGMTで可視化する

19247 ワード
Fortran Tsunami GMT5 Fortran テキストリンク

目的

JAGURSを使って津波伝播シミュレーションんをするための海底地形データの入手と変換と可視化する。

実行環境

OS win10 64bit
Windows Subsystem for Linux (Ubuntu)
GMT5
gfortran

海底地形データの入手

G空間情報センター(公式)に行って、アカウントを作ってください。
「データセット」(公式)行って「津波」で検索して、左側の、「内閣府　中央防災会議　日本海溝・千島海溝周辺海溝型地震に関する専門調査会」をクリックして、

地形データを選択して、
image.png

全部のファイルを詳細→ダウンロードでダウンロードしてください（くっそ使いにくい！）

データフォーマット

とりあえず、なんでも良いのでダウンロードしてきてzipを解凍したデータをメモ帳等で表示してみましょう。
image.png
(ﾟДﾟ)ﾊｧ?
位置情報が無いんですけどー！

DOCファイルを見ると、Fortran書式の(10f8.2)だそうです。Fortranではよくあるファイル書き出し形式です。（使いにくいんだよ馬鹿野郎）

変換するにも、可視化するにも、この形式はくっそ不便なので、とりあえず(i,j,depth)形式と(東西,南北,depth)に変換します。
変換するための基本データはdocに画像として保存されています（常識的に考えて、csvファイルつけるよね？）
ちなみに、南海トラフのデータにはエクセルファイルがついてた気がする。
ありがとう！online ocr (公式)

丸1年ぶりにFortranというよりコンパイルして動かすでプログラムを書いた。。。。

convert_to_xyz.f
      implicit none
      integer i,j,k
      integer px,py
      integer ix,iy,dxdy
      character(len=200) fname,iformat
      integer iosnum      
      open(11,file='bathylist.csv',status='old')
      do
        read(11,*,iostat=iosnum)fname,dxdy,px,py,ix,iy
        if (iosnum < 0)then
          exit
        end if
        write(*,*) trim(fname)
        iformat = '(10f8.5)'
        call convertfile(fname,dxdy,px,py,ix,iy,iformat)
      end do

      close(11)
      stop
      end

      subroutine convertfile(infname,dxdy,px,py,ix,iy,iformat)
        implicit none
        character(len=200) infname,iformat
        character(len=200) fname,fname2,fname3
        integer dxdy,ix,iy,px,py
        integer i,j
        real z(ix,iy)

        iformat = '(10f8.3)'
        fname = 'depth_'//trim(infname)//'.dat'
        write(*,*)fname
        open(21,file=trim(fname),status='old')
        read(21,fmt=trim(iformat))((z(i,j),i=1,ix),j=1,iy)
        close(21)
        fname2 = 'depth_'//trim(infname)//'_ijz.dat'
        open(22,file=trim(fname2),status='unknown')
        do i=1,ix
          do j=1,iy
            write(22,'(i6,i6,f15.6)')i,j,z(i,iy-j+1)
          end do
        end do
        close(22)

        fname3 = 'depth_'//trim(infname)//'_xyz.dat'
        open(23,file=trim(fname3),status='unknown')
        do i=1,ix
          do j=1,iy
            write(23,'(f12.2,f12.2,f15.4)')real((i-1)*dxdy)+px,
     c                                     real(-(j-1)*dxdy)+py,
     c                                     z(i,j)
          end do
        end do
        close(23)
      end subroutine convertfile

二重のdo文が重なったread文なんてFortranでしか見たことありません。

doubleLoopReadMethod.f
read(unit=11,fmt='(10F8.2)')((z(i,j),i=1,ix),j=1,iy)

ファイルには一行で１０個の数字しかないのですが、こうやることで水深データがz(1000,1350)の配列に収まります。
ただし事前に配列の大きさ(ix,iy)がわかってる必要があります。

参考
http://www.nag-j.co.jp/fortran/FI_17.html#ImplDos

とりあえず海底地形の可視化

使ったのがGMT5
ちなみに僕は、GMT4しか使ったことがないので、不安でいっぱいです。

gmt xyz2grd epth_1350-01_ijz.dat -Gepth_1350-01_ijz.grd
gmt makecpt -Crainbow -T-3000/9000/500  > colors.cpt
gmt grdimage depth_1350-01_ijz.grd -Jx0.01 -Ccolors.cpt -V -P -Ba > aaa.eps
gmt psconvert aaa.eps  -A+s4c -P -V -Tg

image.png
無事に東日本とその左側の海底が表示されています。たったこれだけの絵を出すために、まるまる一晩かかりました。原因は、いつの間にかGMT5がEPS出力を辞めてたことと、一部しか表示されてないことに気づかずにあーでもないこーでもないしてました。
あと、GSViewerが見つからなかったこと。

まとめ

僕はデータフォーマットやGMTにさんざん文句を言っていますが、約１５年前にもまったく同じ経験をして全く同じ文句を言っています。
経験的にGMTの調整には一晩かかることも知ってます。
僕って進歩ないなぁ