LoginSignup
0
5

More than 3 years have passed since last update.

@k_reiji

【Python, ObsPy】Matplotlib+ObsPyで震源球(beachball)を書いてみた

Last updated at Posted at 2021-01-02

地震の発震機構を表す震源球 ( beachball ) を描画するツールとしては、GMT が一般的かと思われますが、Python 及びその中の可視化ライブラリ Matplotlib を用いてもできます。Python を用いて震源球を描画する主な方法としては、調べた限りでは次の通りです。

  1. 震源球の投影方法を勉強して自力でスクリプトを描く。
  2. 自前環境に GMT (ただし>=6.1.11) がインストールされている状態で、Pythonのライブラリの一つである PyGMT をインポートし、そのAPIの一つである pygmt.Figure.meca を用いて GMT を呼び出す。
  3. QGISを用いる
  4. Python のライブラリの一つである ObsPy を用いる。

GMT については現役の研究者等が書いた指南書が腐るほどあり、QGISについては解説記事が既にありますので、ここでは 4. の ObsPy について、単体の震源球を描画する方法を紹介します。

※表示している地震のデータは防災科学技術研究所が公開している地震のメカニズム ( 防災科研：メカニズム解の検索 ) からランダムに選びダウンロードしたものです。本ページ中の画像は個人がデモンストレーションとして描画したもので、その結果に学術的な意味を求めていません。

※本ページのスクリプトの利用は自己責任でお願いします。

描画例
some.png

1. 目次

1. 目次
2. このページの目的
3. 使用環境
4. ObsPyのインストール方法
5. とりあえずなにかを描画
6. デモンストレーション
7. おわりに
8. 参考文献

2. このページの目的

Matplotlib の subplot として一つの震源球を書くこと。

3. 使用環境

  • Windows 10
  • conda 4.9.2
  • python 3.8.5
  • matplotlib 3.3.2
  • obspy 1.2.2

Matplotlib は Anaconda をインストールした際にデフォルトで入ったものを使用しました。ObsPy の required として future, scipy, sqlalchemy, numpy, lxml, setuptools, decorator, matplotlib, requests などが表示されていたので、インストールでエラーが出たらこれらを疑いましょう。通常手順でインストールすればこれらも自動でインストールされると思います。

4. ObsPyのインストール方法

公式は Anacondaを使って仮想環境にインストールすることを推奨しています2。具体的には、

  1. コマンドでconda config --add channels conda-forgeと打って Anaconda の configuration に conda-forge チャンネルを加え、
  2. conda create -n obspy python=3.8 等で指定したpythonのバージョンが動く仮想環境( 環境名：obspy )を作り、
  3. conda activate obspy 等で仮想環境を起動し
  4. conda install obspy でインストールする。

という流れになります。詳しくは「conda 仮想環境 作り方」や「conda パッケージ インストール方法」で調べれば出てきます。上記の流れは下記のページをそのまま訳したものなので、やる前に確認しておくことを強く推奨します。

Installation via Anaconda · obspy/obspy Wiki · GitHub

なお自前環境では pip を用いて一発で入りました ( 壊れてもいいパソコンなので )。コマンドプロンプトでpip install obspyと打ち込めば依存性もろとも処理してくれました。ただし Windows で環境構築する場合には、仮想でやった方が一般によいという内容がネット上に散見されますので、おすすめはしません。またより上のバージョンが欲しい場合は下記リンクにgithubのレポジトリがありますので、そこからインストールできるかもしれません ( 筆者は未経験 )。

GitHub - obspy/obspy: ObsPy: A Python Toolbox for seismology/seismological observatories.

5. とりあえずなにかを描画

ObsPy の公式チュートリアル (リンク 1、リンク 2) を参考にしながら、以下のスクリプトを実行してみます。自前環境では matplotlib.pyplot.show() がうまく動かなかったので、とりあえず PNG の画像 (test.png) を作成するように組みました。

test.py
import matplotlib.pyplot as plt
from obspy.imaging.beachball import beach

fig = plt.Figure(figsize=(6,6)) # matplotlib.pyplot.Figure 600x600(px)
np1 = [249,15,127] # 描きたい発震機構その1
mt = [-1.2676, 1.8667, -0.5991, 0.7903, 1.6198, -0.9791] # 描きたい発震機構その2
beach1 = beach(np1, xy=(25, 25), width=30) 
beach2 = beach(mt, xy=(75, 75), width=50)
ax = fig.add_subplot(111, facecolor="white")
ax.plot([25, 75], [75, 25], "rv", ms=20)
ax.add_collection(beach1)
ax.add_collection(beach2)
ax.set_xlim((0, 100))
ax.set_ylim((0, 100))
fig.savefig("test.png")
plt.close()

結果
test.png

環境が整っていれば難なく実行できるはずです。上の画像からもわかる通り、次のような特徴があります。

  1. 設定した座標値に震源球の中心が来るように配置される。
  2. デフォルトでは青色になるので、よく使われる白黒表示にしたい場合は色をいちいち黒に変更しないといけない(環境によっては異なるかも)。
  3. 断層面の構成 ( 走向、傾斜、すべり角 ) を与える場合とモーメントテンソル ( mrr(mzz), mtt(mxx), mff(myy), mrt(mxz), mrf(-myz), mtf(-mxy) ) を与える場合の２種類の方法が用意されている。
    ※要素の順番、方位については公式が上記の順番を採用しているらしい ( obspy.imaging.beachball.beach — ObsPy Documentation (1.2.0) ) が、再確認を推奨する。特にモーメントテンソルの書式については複数通りあるので要注意。
  4. デフォルトではおそらく下半球投影である(要確認)。
  5. 震源球の大きさは座標値に依存している ( つまりx軸y軸のスケールが違うと震源球が歪む )。

6. デモンストレーション

チュートリアルのスクリプトを踏まえて、デモ用のスクリプトを作成しました。

beachball_demo.py
import matplotlib.pyplot as plt
from obspy.imaging.beachball import beach

fig = plt.Figure(figsize=(18,6))

# 走向、傾斜、すべり角による描画
ax1 = fig.add_subplot(1,3,1, facecolor="white") # 背景色は白く
np1 = [249,15,127] # strike, dip, rake
explanation = "Earthquake\n2003-09-26 04:50" # 震源球の上に記載するテキスト
explanation2 = "strike : " + str(int(np1[0])) + '\n' + "dip : " + str(int(np1[1])) + '\n' + "rake : " + str(int(np1[2])) #同様 
beach1 = beach(np1, xy=(50, 50), width=50, facecolor="black") # obspy.imaging.beachballによる描画アイテム 自前環境だとデフォルトでは青
ax1.set_aspect("equal")# xy軸のスケールをそろえる
ax1.add_collection(beach1) # beachballを挿入
ax1.text(50,85,explanation,fontsize=20,horizontalalignment='center',verticalalignment='center') # テキストの配置。ボックスの中心が(50,85)になるように
ax1.text(50,12,explanation2,fontsize=20,horizontalalignment='center',verticalalignment='center') # 同様
ax1.set_xlim((0, 100)) # 便宜上、このような座標値を与えている
ax1.set_ylim((0, 100)) # 同様
ax1.axes.xaxis.set_visible(False) # 目盛の非表示
ax1.axes.yaxis.set_visible(False) # 同様

# モーメントテンソルによる描画
ax2 = fig.add_subplot(1,3,2, facecolor="white")
mt = [-1.2676, 1.8667, -0.5991, 0.7903, 1.6198, -0.9791] # mrr(mzz) mtt(mxx) mff(myy) mrt(mxz) mrf(-myz) mtf(-mxy)
explanation = "Earthquake\n2020-12-31 23:45"
explanation2 = "mrr : " + str(mt[0]) + '\n' + "mtt : " + str(mt[1]) + '\n' + "mff : " + str(mt[2])
explanation3 = "mrt : " + str(mt[3]) + '\n' + "mrf : " + str(mt[4]) + '\n' + "mtf : " + str(mt[5])
beach2 = beach(mt, xy=(50, 50), width=50, facecolor="red")# 赤色にする
ax2.set_aspect("equal")
ax2.add_collection(beach2)
ax2.text(50,85,explanation,fontsize=20,horizontalalignment='center',verticalalignment='center')
ax2.text(2,12,explanation2,fontsize=20,horizontalalignment='left',verticalalignment='center') # 左揃えに 座標値も少し調整
ax2.text(52,12,explanation3,fontsize=20,horizontalalignment='left',verticalalignment='center')
ax2.set_xlim((0, 100))
ax2.set_ylim((0, 100))
ax2.axes.xaxis.set_visible(False)
ax2.axes.yaxis.set_visible(False)

#適当な散布図 地震とは関係がない
ax3 = fig.add_subplot(1,3,3, facecolor="white")
ax3.scatter([10, 49, 46, 90, 24], [45, 57, 32, 67, 2],color="black")
ax3.set_aspect("equal")
ax3.set_xlim((0, 100))
ax3.set_ylim((0, 100))
ax3.set_xlabel('x', fontsize=20)
ax3.set_ylabel('y', fontsize=20)
ax3.set_title("Some plot", fontsize=20)

plt.tight_layout()# ラベルとかが重ならないように
fig.savefig("beachball_demo.png") # png画像として保存
#fig.show() # 自前環境だと発動しなかった
plt.close()

結果
beachball_demo.png

うまく動くみたいですv(-∀-)v

7. おわりに

このスクリプトを応用することにより、ほかのデータをプロットしている画像内で震源球を描画できるかと思います。ただ震源球の需要の大部分は地図にプロットすることなので、Python上でそれらを完結させるとしたら地図ライブラリである Cartopy などと連携させる必要があります ( 筆者は現状ではそこまで追えていません ) 。
地図へのプロットについては ObsPy の公式ドキュメントでは mpl_toolkits.basemap を使用したスクリプトが紹介されています ( リンク ) が、使用している mpl_toolkits.basemap が 2020年 をもって新規リリースが停止される3ことになっており、後継として Cartopy の利用が推奨されています。機会があれば筆者が試してみます。

地図上の描画については現状 GMT やQGISのほうが参考文献が多いと思いますので、それらを参照することをお勧めします。

=====
2021/01/15追記：Cartopyとの連携が可能であることが分かりました。
【Python, ObsPy】Cartopy+ObsPy で地図上に震源球 (beachball) を描いてみた

8. 参考文献

  1. Installing — PyGMT ↩

  2. Installation via Anaconda · obspy/obspy Wiki · GitHub ↩

  3. Introduction — Basemap Matplotlib Toolkit 1.2.1 documentation ↩

0
5
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
Sign upLogin
0
5