gmt project は…まあ、色々なことができるモジュールです。
簡単に言うとデータを大円に沿って投影することができます。私自身も全てを理解していません。
詳しい説明は公式ドキュメントを読んでください。
任意の線を用いて色々する
ここでは、project のうち、任意の線を用いて色々する機能を説明します。
自分が唯一使いこなせる機能です。
- 始点から終点まで、○○km間隔の経度・緯度を求める。
gmt project -C[始点経度]/[始点緯度] -E[終点経度]/[終点緯度] -G[データを取得する間隔] -Q
※ -Qを付けると単位が [km] で統一される。
- データを任意の範囲内で絞り込む。
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
となります。
謝辞
震源データには気象庁一元化震源データを改変したものを使用しました。
ここに記して感謝申し上げます。
コメント