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

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

#実行環境
OS win10 64bit
Windows Subsystem for Linux (Ubuntu)
GMT5
gfortran

#海底地形データの入手

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

![image.png](https://qiita-image-store.s3.amazonaws.com/0/180426/90bd376e-f0f6-7cc2-7251-159eeaf402d0.png)

地形データを選択して、
![image.png](https://qiita-image-store.s3.amazonaws.com/0/180426/3727fb2d-437a-062c-a070-444ba35cc31b.png)

全部のファイルを詳細→ダウンロードでダウンロードしてください（くっそ使いにくい！）
![image.png](https://qiita-image-store.s3.amazonaws.com/0/180426/4a3a7231-81d0-fc53-6b10-9a7946c189f0.png)


#データフォーマット
とりあえず、なんでも良いのでダウンロードしてきてzipを解凍したデータをメモ帳等で表示してみましょう。
![image.png](https://qiita-image-store.s3.amazonaws.com/0/180426/06ad6ae5-62bb-435a-b9b7-aa48b39c5f6c.png)
(ﾟДﾟ)ﾊｧ?
位置情報が無いんですけどー！

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

![image.png](https://qiita-image-store.s3.amazonaws.com/0/180426/cf4e1a66-db97-7bbf-d66a-1b922ddbb1a2.png)

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

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

```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でしか見たことありません。

```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しか使ったことがないので、不安でいっぱいです。

```bash
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](https://qiita-image-store.s3.amazonaws.com/0/180426/05b9d310-9eec-8b84-4ab6-6a0ddffd13d6.png)
無事に東日本とその左側の海底が表示されています。たったこれだけの絵を出すために、まるまる一晩かかりました。原因は、いつの間にかGMT5がEPS出力を辞めてたことと、一部しか表示されてないことに気づかずにあーでもないこーでもないしてました。
あと、[GSViewer](http://www.ghostgum.com.au/software/gsview.htm)が見つからなかったこと。

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