震源断面図を描く その2

前回の続きです。

今回は南北断面図を描いてみましょう。

前回のスクリプト

震源データは気象庁一元化震源データを使用しました。
記して感謝申し上げます。

#!/bin/bash
gmt begin foo png
    gmt basemap -JM12 -R132/135/32/34 -Bafg -BWSNE -Y60
    gmt makecpt -Cgeo -T-8000/8000/200 -Z
    gmt grdcut @earth_relief_01s -R132/135/32/34 -Gdem.nc
    gmt grdgradient dem.nc -Ggrad.grd -A45 -Ne0.8
    gmt grdimage dem.nc -Igrad.grd -C
    gmt coast -Df -W0.25 -LJBR+jTR+o0/1+c31+w100+f
    gmt makecpt -Cseis -T0/100/1 -Z
    gmt plot hypo.txt -Sc0.1 -C -W0.25 -:
    gmt colorbar -DJBR+jBL+o2/0+w-5/0.2 -Baf
    gmt basemap -JX12/-12 -R132/135/0/100 -Bxafg+l"longitude" -Byafg+l"depth" -BWSne -Y-15
    awk '{print $1,$3,$3}' hypo.txt | gmt plot -Sc0.1 -W0.01 -C
    gmt colorbar -DJBR+jBL+o2/0+w-5/0.2 -Baf
gmt end

地図のタテの長さを求める

gmt mapprojectの使い方

GMTのmapprojectモジュールについて説明しています。

dandango.pw

2020.09.03

東西断面図と異なり、南北断面図を描くには地図のタテの長さを求めなければなりません。

#!/bin/bash
gmt mapproject -JM12 -R132/135/32/34 << END
135 34
END

この結果が、12 9.49479714703となるので、地図のタテの長さは9.49479714703cmとなります。

南北断面図の横軸は深さ、縦軸は緯度であることに注意して、

#!/bin/bash
gmt begin foo png
    gmt basemap -JM12 -R132/135/32/34 -Bafg -BWSNE
    gmt makecpt -Cgeo -T-8000/8000/200 -Z
    gmt grdcut @earth_relief_01s -R132/135/32/34 -Gdem.nc
    gmt grdgradient dem.nc -Ggrad.grd -A45 -Ne0.8
    gmt grdimage dem.nc -Igrad.grd -C
    gmt coast -Df -W0.25 -LJBR+jTR+o0/1+c31+w100+f
    gmt makecpt -Cseis -T0/100/1 -Z
    gmt plot hypo.txt -Sc0.1 -C -W0.25
    gmt basemap -JX9.49479714703/-9.49479714703 -R0/100/32/34 -Bxafg+l"depth" -Byafg+l"latitude" -BwsNE -X15
    awk '{print $3,$2,$3}' hypo.txt | gmt plot -Sc0.1 -W0.01 -C
    gmt colorbar -DJBR+jBL+o2.5/0+w-5/0.2 -Baf+l"depth"
gmt end

とすれば、この結果が、

となります。

まとめ

折角なので1枚の図に東西断面図と南北断面図を描いてみよう。

#!/bin/bash
gmt begin foo png
    # 震源分布図
    gmt basemap -JM12 -R132/135/32/34 -Bafg -BWSNE -Y60
    gmt makecpt -Cgeo -T-8000/8000/200 -Z
    gmt grdcut @earth_relief_01s -R132/135/32/34 -Gdem.nc
    gmt grdgradient dem.nc -Ggrad.grd -A45 -Ne0.8
    gmt grdimage dem.nc -Igrad.grd -C
    gmt coast -Df -W0.25 -LJBR+jTR+o0/1+c31+w100+f
    gmt makecpt -Cseis -T0/100/1 -Z
    gmt plot hypo.txt -Sc0.1 -C -W0.25
    # 東西断面図
    gmt basemap -JX12/-12 -R132/135/0/100 -Bxafg+l"longitude" -Byafg+l"depth" -BWSne -Y-15
    awk '{print $1,$3,$3}' hypo.txt | gmt plot -Sc0.1 -W0.01 -C
    gmt colorbar -DJBR+jBL+o5/0+w-10/0.2 -Baf+l"depth"
    # 南北断面図
    gmt basemap -JX9.49479714703/-9.49479714703 -R0/100/32/34 -Bxafg+l"depth" -Byafg+l"latitude" -BwsNE -X15 -Y15
    awk '{print $3,$2,$3}' hypo.txt | gmt plot -Sc0.1 -W0.01 -C
gmt end

この結果が、

となる。それっぽい！

今回はここまで。