Edited at

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


目的

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


実行環境

OS win10 64bit

Windows Subsystem for Linux (Ubuntu)

GMT5

gfortran


海底地形データの入手

G空間情報センター(公式)に行って、アカウントを作ってください。

「データセット」(公式)行って「津波」で検索して、左側の、「内閣府 中央防災会議 日本海溝・千島海溝周辺海溝型地震に関する専門調査会」をクリックして、

image.png

地形データを選択して、

image.png

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

image.png


データフォーマット

とりあえず、なんでも良いのでダウンロードしてきてzipを解凍したデータをメモ帳等で表示してみましょう。

image.png

(゚Д゚)ハァ?

位置情報が無いんですけどー!

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

image.png

変換するにも、可視化するにも、この形式はくっそ不便なので、とりあえず(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)


ファイルには一行で10個の数字しかないのですが、こうやることで水深データが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にさんざん文句を言っていますが、約15年前にもまったく同じ経験をして全く同じ文句を言っています。

経験的にGMTの調整には一晩かかることも知ってます。

僕って進歩ないなぁ