##はじめに
このページに書いてあるスクリプトは未完成です。
GMT5について少しでも噛んでいる方に読んでいただいて、おかしい所を指摘していただけると大変助かります。
#本題
地図全体の震源断面図を作るのは簡単でawkでごにょごにょしてpsbasemap pscoast psxyでどうにかすればいい。
自分はポンコツなのでエレガントには書けないが。
⇒じゃあ、任意の直線を引いてその線の震源断面図を書くにはどうすればよいか。
##Projectコマンド
公式マニュアル曰く、projectを使いなさいとのこと。
http://gmt.soest.hawaii.edu/doc/5.4.5/project.html
project - Project table data onto lines or great circles, generate tracks, or translate coordinates
GMT4版ではあるが日本語に訳していただいたサイトを見てみると
project は任意の (x, y[, z]) データを標準入力[又は infile ]から読み込ん で (x, y, z, p, q, r, s) の任意の組合せを標準出力に書き出す。ここで (p, q) は投影された座標、 (r, s) は (x, y) 座標系における位置で、プロファイ ル (q = 0 の経路) は (x, y) に最も近い位置になり、 z は入力された残り全 ての列である(必要な x と y の列以外)。代わりに、 project はプロファイル に沿った (r, s, p) を等しい間隔 dist で作成することもできる。
http://gmt.shin-gen.jp/GMT4.1.4/project.html
らしい。ざっくりと言うと以下の通り。
- 投影法は次の 3 つのうちの任意の( ただし単独の) 1 つで定義される。
- 定義1. 中心
−Cと北から時計回りの方位角-A。 - 定義2. 投影経路の始点
−Cと終点−E。 - 定義3. 中心
−Cと 回転極の位置−T。
- 定義1. 中心
どうやらこの3つのうちの1つを使いなさいよってことらしい。今回は直線断面図なので定義2を使う。
(他の定義を誰かやってくださると助かります)
データは
−L及び−Wオプションを用いて選択的にある範囲だけを取り出すことができる。−Wを使うと、投影される幅はw_min<q<w_maxの範囲内だけになる。−Lを設定すると、長さがl_min<p<l_maxの範囲内の点だけを使うように設定される。−Eオプションが投影法の定義に使われているときは、−Lwを選ぶと投影範囲はちょうどOからBまでの範囲になる。
つまり、本来は-Lを使って長さを指定しなければいけないが、-Lwにすれば始点から終点まで読み込むということっぽい…？
##結果と課題
というわけで何となく作ってみた。-Wの使い方がよく分からない…。
もっとエレガントな書き方があると思います。ColorPalletなどはご自由に。
#!/bin/bash
output='out.ps'
year='2016'
makecpt -Cseis.cpt -T0/160/10 -Z > depth.cpt
gmt psbasemap -R134/138/33/36 -Jm3.5 -Ba1f0.5g1WSen -X3 -Y16 -Lf137.5/33.2/25/50 -P -K -V > ${output}
gmt pscoast -R -J -Df -W1 -G165/165/165 -K -O -V >> ${output}
#直線を引く始点と終点の設定
lon1='134'
lon2='138'
lat1='33'
lat2='35.5'
#ここら辺最高に無駄な気がする
awk '{printf $1" "$2"\n"}' << END > lonlat.dat
${lon1} ${lat1}
${lon2} ${lat2}
END
#図に直線を引く
gmt psxy lonlat.dat -R -J -W1,red -K -O -V >> ${output}
#気象庁の震源レコードから経度・緯度・深さを抽出
#awkじゃなくても「経度 緯度 深さ」の順番になっているデータを置けば大丈夫なはず
awk '{print $11" "$9" "$13}' h${year}.hypo | \
gmt psxy -Cdepth.cpt -R -J -h1 -Sc0.045 -K -O -V >> ${output}
#地図の下に書く断面図作成
gmt psbasemap -R0/500/0/140 -JX14c/-5c -Bxa100f50g100+l"distance(km)" -Bya20f10g20+l"depth(km)" -BWSen -X0 -Y-7 -P -K -O -V >> ${output}
#気象庁の震源レコードから再度抽出し直線付近の震源のみを絞り込む
#-Wオプションあたりで地震数が決まったりしているので、今度詳しく調べないといけない
awk '{print $11" "$9" "$13}' h${year}.hypo > temp0.dat
gmt project temp0.dat -C${lon1}/${lat1} -E${lon2}/${lat2} -W-30/30 -Lw -Q -V > temp1.dat
#断面図に書き込み
awk '{print $4" "$3" "$3}' temp1.dat | \
gmt psxy -Cdepth.cpt -R -J -Sc0.06 -K -O -V >> ${output}
#ここらへん自分でも怪しい
gmt psscale -Dx15/10/-5/0.2 -Jm -R134/138/33/36 -Cdepth.cpt -By20+l"Depth[km]" -O -V >> ${output}
#いらないやつを削除
rm temp* lonlat.dat depth.cpt
この結果が、
になった。線が赤くなってしまって震源と被ってるが一応書けたは書けたので満足。
##(既知の問題点)
- 横軸の距離がおかしい
##今後の課題
-
-Wとは - distanceの長さの一致
- てかこれ本当に震源断面図があってるのか
⇒詳しい方がいたらご教授願えると幸いです。
GMT（Generic Mapping Tool; Wessel, P. and W. H. F. Smith, New, improved version of Generic Mapping Tools released, EOS Trans. Amer. Geophys. U., Vol. 79 (47), pp. 579, 1998）