gmt project の使い方 / 震源断面図を描く その3

project — GMT 6.2.0 documentation

gmt project は…まあ、色々なことができるモジュールです。

簡単に言うとデータを大円に沿って投影することができます。私自身も全てを理解していません。
詳しい説明は公式ドキュメントを読んでください。

任意の線を用いて色々する

ここでは、project のうち、任意の線を用いて色々する機能を説明します。
自分が唯一使いこなせる機能です。

  1. 始点から終点まで、○○km間隔の経度・緯度を求める。

gmt project -C[始点経度]/[始点緯度] -E[終点経度]/[終点緯度] -G[データを取得する間隔] -Q
-Qを付けると単位が [km] で統一される。

  1. データを任意の範囲内で絞り込む。

gmt project [データファイル名] -C[始点X座標]/[始点Y座標] -E[終点X座標]/[終点Y座標] -W[許容する範囲]/[許容する範囲] -Lw -Q
-Qを付けると単位が [km] で統一される。

○○km間隔の経度・緯度を求める

たとえば、

#!/bin/bash
gmt begin project png
    gmt basemap -JM12 -R139/143/38/41 -Bagf -BWSNE
    gmt makecpt -Cgeo -T-8000/8000/200 -Z
    gmt grdcut @earth_relief_03s -R139/143/38/41 -Gdem.grd
    gmt grdgradient dem.grd -Gmapgrad.grd -A45 -Ne0.8
    gmt grdimage dem.grd -C -Imapgrad.grd
    gmt coast -Df -W0.05 -Lx11/-1.0+c20+w100+f
    gmt plot -W2 -Wred << END
139.5 40
142.5 40
END
gmt end

として、

のとき、赤線の左側の端点から5kmごとの座標を求めるには、

#!/bin/bash
gmt project -C139.5/40 -E142.5/40 -G5 -Q > result.txt

とすれば、result.txt が作成され、「経度・緯度・距離」 の順でデータが記されているはずです。

データを任意の範囲内で絞り込む

たとえば、ここに2019年に起きた地震を描いてみます。
震源データには気象庁一元化震源データを使用しました。

となります。
(震源データはこちらから)

ここで、赤線周辺で発生した地震を絞り込んでみましょう。
厳密に赤線上で絞り込むと震源が少なくなってしまうので、今回は赤線から各々15kmまでは離れていても許容します。
赤線を中心に15kmずつ幅を持たせるイメージです。

#!/bin/bash
 gmt project hypo.txt -C139.5/40 -E142.5/40 -W-15/15 -Lw -Q > danmen.txt

とすれば、 danmen.txt というファイルに絞り込んだ震源データが記載されているはずです。
一応、「震源経度・震源緯度・震源深さ・震源マグニチュード・赤線始点からの距離・赤線からのズレ・赤線上での経度・赤線上での緯度」の順に並んでいるみたいです。

あとは、これを横軸に距離を取って描けばよいので、

#!/bin/bash
distance=`gmt mapproject -G139.5/40+uk << END | awk '{print $3}'
142.5 40
END
`
gmt begin danmen jpg
    gmt basemap -JX12/-12 -R0/${distance}/0/200 -Bxafg+l"distance[km]" -Byafg+l"depth[km]" -BWSne
    gmt makecpt -Cseis -T0/200/1 -Z
    awk '{print $5,$3,$3}' result.txt | gmt plot -Sc0.1 -W0.01 -C
gmt end

この結果が、

となります。

これを地図と組み合わせると、

#!/bin/bash
distance=`gmt mapproject -G139.5/40+uk << END | awk '{print $3}'
142.5 40
END
`
gmt begin danmen jpg
    gmt basemap -JM12 -R139/143/38/41 -Bagf -BWSNE -Y60
    gmt makecpt -Cgeo -T-8000/8000/200 -Z
    gmt grdcut @earth_relief_03s -R139/143/38/41 -Gdem.grd
    gmt grdgradient dem.grd -Gmapgrad.grd -A45 -Ne0.8
    gmt grdimage dem.grd -C -Imapgrad.grd
    gmt coast -Df -W0.05 -Lx11/-1.0+c20+w100+f
    gmt plot -W2 -Wred << END
139.5 40
142.5 40
END
    gmt basemap -JX12/-12 -R0/${distance}/0/200 -Bxafg+l"distance[km]" -Byafg+l"depth[km]" -BWSne -Y-16
    gmt makecpt -Cseis -T0/200/1 -Z
    awk '{print $5,$3,$3}' result.txt | gmt plot -Sc0.1 -W0.01 -C
gmt end

となります。

謝辞

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

コメント

タイトルとURLをコピーしました