GMT projectの練習 伊豆から富山の震源分布 (再)

震源の深さ方向の分布を見るのに便利なprojectというコマンドがGMTでは用意されている、ということを最近知りました。以前の記事ではいったんxyなりutmなりに変換してたので、少し面倒だったんですが、このコマンドを使うと作業が一つ減ります。

　←クリックで拡大

これを稽古するのにこんな断面を気象庁地震カタログから作成して見ました。経度(x)、緯度(y)、深さ(z)、マグニチュードの並びのファイルをある2点を通る大円に沿う方向(p)と、直行する方向の角度(q)で表した座標に変換してくれます。断面はp×地球半径で変換した距離と深さでpsxyをすれば、直接描くことができます。上の図ではqの値で色分けをして奥行きを表現しました。
上図右が四角で囲った部分を東から見た断面での震源分布、左は西から見た断面です。
また、pの値、qの値で出力を制限できるので図の四角の範囲に含まれる震源だけを抽出するということができます。

マニュアルではマグニチュードに相当する説明がありませんが、x,y,[z]の[z]はx,yに続く部分をまとめてzと扱うようで、もとのデータが、x,y,z,mとなっていた場合、出力にzを指定すれば、もれなくmが付いてきます。
それから、データの取り出しは、大円についての角度などで指定できるのですが、それが地図上でどこに来るのか、上図の四角を描くための座標ですが、これは赤道上の2点を投影しておいて、pが0～L,qが-Wから+Wという投影後考えている四角形の座標を今度は投影された赤道上の2点で投影しなおす、ということをやれば元の緯度、経度に戻すことができます。

プロットしたのは一度断面を見てみたい、と思っていた個所なんですが、伊豆半島のところで、地震発生そうの下端が盛り上がってたり、何か色々妄想できそうなものができました。
火山の場所は日本の火山位置データと活火山データベースを参考にしました。


コメントを頂いたんですが、何か手違いがあって、非表示になっていたのを復活しました。