🌏

PyGMTとUSGSのAPIを使って地震の分布を描いてみる

2023/10/14に公開

イントロダクション

教科書やテレビなどで、地震の分布図を一度は目にしたことのある方も多いのではないでしょうか。この記事では、USGS（アメリカ地質調査所）の公開している地震情報と、PythonのライブラリであるPyGMTを使って作図する方法を紹介します。

2018〜2022年のマグニチュード5より大きい地震（データ出典：USGS）
地震情報をダウンロードするには、USGSや気象庁のサービスを利用するなどの方法がありますが、USGSの提供しているAPIを利用すると、比較的に簡単に情報を取得することができます。

また、Pythonで地図上にデータをプロットするのにはいくつか選択肢がありますが、今回は地球科学などの分野でよく使われてきたCUIのソフトウェアであるGMT（Generic Mapping Tools）のPythonラッパー、PyGMTを使用しました。

USGSのAPIと本家GMTを使用して、シェルスクリプトで可視化する方法を紹介してくださっている記事はありましたが、USGS API+PyGMTでの例は見当たらなかったので、この記事がどなたかの参考になれば幸いです。コードはGitHubにもアップロードしています。

動作環境

PyGMTを実行できる環境が必要になります。今回はこちらの情報を参考に、Google Colaboratoryを使用しました。下記のコードを実行してPyGMTをインストールします。

!pip install -q condacolab
import condacolab
condacolab.install()
!mamba install pygmt

なお、他の環境でPyGMTを導入するには、公式のガイドやこちらの情報が参考になりそうです。

USGSのAPI

USGSが公開している仕様を参考に、今回はこのようにURLを指定して取得してみました。

https://earthquake.usgs.gov/fdsnws/event/1/query?format=geojson&starttime=2018-01-01&endtime=2022-12-31&minmagnitude=5&eventtype=earthquake

"format="以降がパラメータを指定した部分になります。

format：geojson形式を指定
starttime, endtime：期間を2018/1/1〜2022/12/31に指定
minmagnitude：マグニチュード5よりも大きい地震を指定
eventtype：デフォルトだと火山現象などでの地震記録も含まれるので、狭い意味での地震のみを取得するのに"earthquake"を指定

ちなみに、これとは別にリアルタイムのフィードもあり、過去1時間、1日、7日、30日の地震情報を取得することもできます。例えば下記のURLで、過去30日間のマグニチュードが4.5より大きい地震情報が取得できます。

https://earthquake.usgs.gov/earthquakes/feed/v1.0/summary/4.5_month.geojson

PyGMT

前述のとおり、PyGMTはGMTのPythonラッパーです。Pythonには可視化ライブラリが多数ありますが、PyGMTはこちらの説明にあるように、地図投影やグリッドデータの処理を扱う機能が充実しています。この記事の内容の範囲では、他の可視化ライブラリと比較したときの強みを十分発揮できていないと思いますが、綺麗な図を作成できると感じていただけたら嬉しいです。

なお、PythonにはBokehやPlotly、Foliumなどの、インタラクティブな図を描画できるライブラリもあります。そうしたケースでは、PyGMTではなく別のライブラリを選択した方が良いでしょう。

コードと実行結果

同じ大きさの円でプロット

APIで取得したデータから、個々の地震の緯度・経度・深さの情報を取り出してプロットしてみます。描画する部分の詳細については公式ドキュメントを参照していただきたいのですが、流れとしては下記のようになります。

fig.coast：地図の投影法や枠線のレイアウトを設定し、海岸線を描画
pygmt.makecpt：カラーパレットの作成
fig.plot：作成したカラーパレットを使って地震データをプロット
fig.colorbar：カラーバーを表示

import pygmt
import requests

url = 'https://earthquake.usgs.gov/fdsnws/event/1/query?format=geojson&starttime=2018-01-01&endtime=2022-12-31&minmagnitude=5&eventtype=earthquake'
json = requests.get(url).json()
earthquake_data = json['features']

coordinates = [quake['geometry']['coordinates'] for quake in earthquake_data]
longitudes = [coord[0] for coord in coordinates]
latitudes = [coord[1] for coord in coordinates]
depths = [coord[2] for coord in coordinates]

fig = pygmt.Figure()
fig.coast(
    projection='Cyl_stere/20c',
    region='g',
    resolution='l',
    land='gray',
    shorelines='1/thinnest',
    frame='a60f30'
)
pygmt.makecpt(
    cmap='no_green',
    series='0/700/50',
    reverse=True)
fig.plot(
    x=longitudes,
    y=latitudes,
    style='c0.3',
    pen='thinnest',
    cmap=True,
    fill=depths
)
fig.colorbar(
    position='JBR+jBL+o0.8/0+w-5/0.3+e',
    frame='a100f50+l"Depth (km)"'
)
fig.show()

上記の図では情報を取得したままの順番でプロットしたので、深い地震は数が少ないため後ろに隠れがちになってしまいました。そこで、pandasライブラリを使って深さでソートしてからプロットしてみます。PyGMTに関連するコードで変わっているのはfig.plotで引数にとっているx, y, fillの値をpandasのseriesに変えている部分のみです。

import pandas as pd

df = pd.json_normalize(json['features'])

df['longitude'] = df['geometry.coordinates'].apply(lambda coord: coord[0])
df['latitude'] = df['geometry.coordinates'].apply(lambda coord: coord[1])
df['depth'] = df['geometry.coordinates'].apply(lambda coord: coord[2])

df_sorted = df.sort_values(by='depth')

fig = pygmt.Figure()
fig.coast(
    projection='Cyl_stere/20c',
    region='g',
    resolution='l',
    land='gray',
    shorelines='1/thinnest',
    frame='a60f30'
)
pygmt.makecpt(
    cmap='no_green',
    series='0/700/50',
    reverse=True)
fig.plot(
    x=df_sorted['longitude'],
    y=df_sorted['latitude'],
    style='c0.3',
    pen='thinnest',
    cmap=True,
    fill=df_sorted['depth']
)
fig.colorbar(
    position='JBR+jBL+o0.8/0+w-5/0.3+e',
    frame='a100f50+lDepth (km)'
)
fig.show()

このような図にすると、プレートが沈み込んでいるエリアでは（日本列島付近など）浅い地震から深い地震まで連続的に分布していることが分かりやすいですね。

マグニチュードに応じて大きさを変えた円でプロット

今度は、APIで取得したデータから個々の地震のマグニチュードの値も取り出して、マグニチュードが大きいほど大きな円になるようにプロットしてみます。

大きな円が背面側になった方が見やすいので、マグニチュードの降順でソートしてプロットします。先ほどのようにpandasでソートしても良いのですが、APIで取得するときに"&orderby=magnitude"をURLに含めると、マグニチュードの降順で情報を得ることができます（ちなみに"&orderby=magnitude-asc"で昇順にできます）。

また、円のサイズを可変にするため、PyGMTでfig.plotのstyleを'c0.3'から'c'に変更し（シンボルの形のみ指定）、新たにsizeでマグニチュードに応じた適当な値のリストを与えています。

url = 'https://earthquake.usgs.gov/fdsnws/event/1/query?format=geojson&starttime=2018-01-01&endtime=2022-12-31&minmagnitude=5&eventtype=earthquake&orderby=magnitude'

json = requests.get(url).json()
earthquake_data = json['features']

magnitudes = [quake['properties']['mag'] for quake in earthquake_data]

coordinates = [quake['geometry']['coordinates'] for quake in earthquake_data]
longitudes = [coord[0] for coord in coordinates]
latitudes = [coord[1] for coord in coordinates]
depths = [coord[2] for coord in coordinates]

fig = pygmt.Figure()
fig.coast(
    projection='Cyl_stere/20c',
    region='g',
    resolution='l',
    land='gray',
    shorelines='1/thinnest',
    frame='a60f30'
)
pygmt.makecpt(
    cmap='no_green',
    series='0/700/50',
    reverse=True)
fig.plot(
    x=longitudes,
    y=latitudes,
    style='c',
    pen='thinnest',
    size=[0.01 * ( 2 ** magnitude ) for magnitude in magnitudes],
    cmap=True,
    fill=depths
)
fig.colorbar(
    position='JBR+jBL+o0.8/0+w-5/0.3+e',
    frame='a100f50+l"Depth (km)"'
)
fig.show()

特に大きな地震はプレートの境界付近で起こっているように見えることが分かります。

謝辞

下記に掲載されているコードを参考にさせていただきました。

https://dandango.pw/archives/606