🌏

PyGMTと気象庁のJSONを使って地震の分布を描いてみる

2023/11/14に公開

Python

pygmt

tech

イントロダクション

以前の記事では、USGSの提供する地震情報を使って地震の分布を描きましたが、今回は気象庁の地震情報を使って描いてみます。

無料で利用できる気象庁の地震情報には下記のようなものがあります。

情報	対象期間	特徴
震度データベース	1919年〜	リアルタイム性はないが、詳細な条件で情報を絞り込むことができる。CSV形式で出力可能
Atomフィード	過去数日間	リアルタイム性があり、詳細な仕様が公開されているので、自動処理に適している。XML形式
気象庁Webページ表示用JSON	過去30日間	リアルタイム性はあるが、APIとして提供されているわけではないので、告知なしに提供中止や仕様変更の可能性あり。JSON形式

その中でも今回は、気象庁Webページ表示用のJSONを使用します。公式に提供されているわけではないので安定性は低いですが、過去30日間のデータが手軽に取得できるのは魅力的です。

コードはGitHubにもアップロードしています。

動作環境

前回同様、Google Colaboratoryを使用しました。
下記のコードを実行してPyGMTをインストールします。

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

コードと実行結果

import re
import requests
import pygmt

def fetch_earthquake_data(url):
    response = requests.get(url)
    data = response.json()
    return [quake for quake in data if quake['ttl'] == '震源・震度情報']

def parse_coordinates(coord_string):
    match = re.match(r'([+|-]\d{2}.\d+)([+|-]\d{3}.\d+)(.*)/', coord_string)
    return tuple(map(float, match.groups()))

url = 'https://www.jma.go.jp/bosai/quake/data/list.json'
quakes = fetch_earthquake_data(url)

coordinates = [parse_coordinates(quake['cod']) for quake in quakes]

latitudes, longitudes, depths = zip(*coordinates)
depths = [-0.001 * depth for depth in depths]

region = [
    min(longitudes) - 1,
    max(longitudes) + 1,
    min(latitudes) - 1,
    max(latitudes) + 1,
]

fig = pygmt.Figure()
fig.coast(
    projection='M15c',
    region=region,
    resolution='l',
    land='gray',
    shorelines='1/thinnest',
    frame='af'
)
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()

コードの説明

地震情報を取得している部分です。
各情報の種別がttlキー（titleの略？）の値として格納されていることを利用して、今回必要な「震源・震度情報」に該当する情報のリストを作成しています。

def fetch_earthquake_data(url):
    response = requests.get(url)
    data = response.json()
    return [quake for quake in data if quake['ttl'] == '震源・震度情報']
    
url = 'https://www.jma.go.jp/bosai/quake/data/list.json'
quakes = fetch_earthquake_data(url)

震源の座標を取り出している部分です。
各情報の震源座標はcodキー（coordinateの略？）の値として格納されているのですが、+32.9+130.9-10000/のような表記になっているのでこれを正規表現を使って取り出し、float型に変換しています。
深さはm単位からkm単位に変換し、正負を反転させます。

def parse_coordinates(coord_string):
    match = re.match(r'([+|-]\d{2}.\d+)([+|-]\d{3}.\d+)(.*)/', coord_string)
    return tuple(map(float, match.groups()))
    
coordinates = [parse_coordinates(quake['cod']) for quake in quakes]

latitudes, longitudes, depths = zip(*coordinates)
depths = [-0.001 * depth for depth in depths]

地震の分布を描画している部分です。
前回の記事から変更している箇所は下記のとおりです。

描画範囲（region）をデータの存在する領域のみに変更
地図の投影法（projection）をメルカトル図法に変更
地図のフレームの目盛り間隔（frame）をPyGMTに自動で選択してもらうよう変更（データに応じて描画範囲が変動するため）

region = [
    min(longitudes) - 1,
    max(longitudes) + 1,
    min(latitudes) - 1,
    max(latitudes) + 1,
]

fig = pygmt.Figure()
fig.coast(
    projection='M15c',
    region=region,
    resolution='l',
    land='gray',
    shorelines='1/thinnest',
    frame='af'
)

参考資料

PyGMTでの描画について、下記の説明を参考にしました。

イントロダクション

動作環境

コードと実行結果

コードの説明

参考資料

Discussion