気象庁の地震データを見様見真似でPythonでスクレイピングしてみた。

はじめに

　地震の発生頻度と規模の関係を表す法則であるグーテンベルグ・リヒター則はご存知でしょうか。

　此の図は, 「気象庁の地震データを見様見真似でPythonでスクレイピングしてみた」結果得られたデータから LuaLaTeX 上の MetaPost に描いて貰ったものです。線形な関係に見えない？そりゃまぁ当然のことで, あとのプログラムや気象庁のデータを実際に眺めて貰えば判るとおり, 近年はマグニチュードが0以下の地震も補足してますが昔はそうでもないし, 実際古いデータに関しては補足している地震が多分少ない（観測地点や状況がそもそも貧弱？）からかデータが少ないわけです。

　流石に規模の大きい地震だとほぼ正確だろうし, 各所に書かれている通り, 「1つの地震」か「幾つかの地震の総体」かも判然としないところはあるので, データを扱う際には注意したいところです。それでもマグニチュード5.5辺りから先を見れば, 何となく線形にデータが並んでいる様子が見えませんか？

---

　さて, この話は, 対数を学んで貰うときに, その実用として紹介していて, 実際ウケも良いしなるほど味も大きいのですが... で, そのデータとしてずっと理科年表にあったデータを利用していたのですが, これが当然だけど古くなっていました。

　もちろん神戸の震災や東日本大震災があったことも関係しているのかもしれませんが, あまりに古いデータだと説得力が欠けそうです。そこで, 最新とは言わないまでも最近までの日本の地震のマグニチュードのデータはないものかと思いましたがこれがなかなか見当たらない。何年か前に手作業的なデータのダウンロードと小手先のテキスト処理で気象庁のデータから数えたこともあったのですがそれも古くなって... しかもこんな事毎年やってられへんと。
　最近流行のLLMな人工知能擬きにお願いするという手も考えたけど... そこまで使った事ないしあてになるかどうか... というわけでスクレイピングだと。

参考にしたもの

気象庁の地震データをPythonでスクレイピングする方法
urllib.request や requests でダウンロードしたzipファイルを一時ファイルを介さずに解凍する
　データ量としては, 2つ目のzipファイルを保存解凍してっていう方で、何十年分のデータを処理できて, とても有難かったのですが, スクレイピングの技術的な点では最初のものがとても有難かった。それは記事にもあるとおり気象庁の直近何ヶ月かのデータは妙なページ構成になっているからです。まるでデータ利用を邪魔しているかのような... 少なくともグーテンベルグ=リヒター則が気になった小中学生や高校生にはデータの利用は難しそうです。

早速ですがPythonのコードを

　普段のプログラミング的なことはRubyやJuliaが主で, 文書などを作成することも多いので図や文書はLuaTeXやMetaPostで済ませています。Pythonも使うのですが... 暫く触ってない間にMacOSXのアップデートのせいらしいですが, VScodeでの利用もターミナルでの反応も変になっていて, この遊び?の為に, 四苦八苦してしましました。結局のところVScode周りは後回しにして, AnacondaからのJupyterのアップデートというか再導入から初めて, Jupyter Notebook上で処理をしました。
　あとは参考にしたページはどちらもデータというかファイルの取得まででその処理の部分がなかったので, そちらは適当に追加してます。

見様見真似コード

import io
import zipfile
import urllib.request
import math
# ここまでがzip展開して処理する為の分
import numpy as np
import pandas as pd
import requests
from bs4 import BeautifulSoup
from tqdm import tqdm
# こちらは謎のHTMLの処理に関する為の分

extract_dir = "/Volumes/Macintosh HD - Data/Users/satie/Jupyter"
# ここは私のMacのディレクトリ？なのでご利用の際は各自で適宜...
# 以下同様です...（隠すような情報でもないのでそのままにしてます。）
eqMcl=[0]*100

# ------------------------------------------------------------------
# 震源リストに掲載されている分の年月日をshingen_listへ格納
URL_index = "https://www.data.jma.go.jp/svd/eqev/data/daily_map/index.html"
res = requests.get(URL_index)
soup = BeautifulSoup(res.content, "html.parser")
shingen_list = [i.get("href") for i in soup.find_all("a") if len(i.get("href"))==13]

# レコード形式
format_dict_1 = {"年":4, "月":2, "日":2, "時分":5, "秒":4, 
                 "緯度":10, "経度":10, 
                 "深さ(km)":7, "M":5, "震央地名":34}

## 入力された1行のレコードを分割する関数を定義
def split_record_1(record):
    record_list = []
    cnt = 0
    for p,c in enumerate(record[59:]):
        if c == " ":
            break

    for name,i in format_dict_1.items():
        record_list.append(record[cnt:cnt+i])
        cnt += i+1
    return record_list

# スクレイピング
URL1 = "https://www.data.jma.go.jp/svd/eqev/data/daily_map/"
df_1 = []

def is_num(s):
    try:
        float(s)
    except ValueError:
        return False
    else:
        return True
    
## ページでforループ
for day in tqdm(shingen_list):
#    print(day)
    url = URL1 + day

    res = requests.get(url)
    soup = BeautifulSoup(res.content, "html.parser")
    r = soup.pre.text  # レコード部分を抽出

    # ループで1行ずつ取り出す
    eqMcl=[0]*100
    for record in [i for i in r.split("\n") if len(i)>1][2:]:
        splr=split_record_1(record)
        df_1.append(splr)
        if is_num(splr[8]):
            mg=int(float(splr[8])*10)
            if mg>0 and mg<100:
                eqMcl[mg]=eqMcl[mg]+1
            else:
                if mg>0:
                    print(splr)

# DataFrameに格納
df_1 = pd.DataFrame(df_1, columns=format_dict_1.keys())
# ------------------------------------------------------------------
# ------------------------------------------------------------------
for mth in range(1,4):
    url = "https://www.data.jma.go.jp/eqev/data/bulletin/data/hypo/h20220"+str(mth)+".zip"
    with (
        urllib.request.urlopen(url) as res,
        io.BytesIO(res.read()) as bytes_io,
        zipfile.ZipFile(bytes_io) as zip,
    ):
        zip.extractall(extract_dir)
        f = open("/Volumes/Macintosh HD - Data/Users/satie/Jupyter/h20220"+str(mth),'r')
        dcnt=0
        for data in f:
            if data[52].isdigit():
                mgn=int(data[52:54])
#                print(mgn)
                if mgn>0 and mgn<100:
                    eqMcl[mgn]=eqMcl[mgn]+1
        f.close()
        print("h20220"+str(mth)+" : ")
        for mgn in range(10):
            print(eqMcl[mgn*10:mgn*10+9])
for yer in range(1998,2022):
    url = "https://www.data.jma.go.jp/eqev/data/bulletin/data/hypo/h"+str(yer)+".zip"

    with (
        urllib.request.urlopen(url) as res,
        io.BytesIO(res.read()) as bytes_io,
        zipfile.ZipFile(bytes_io) as zip,
    ):
        zip.extractall(extract_dir)
        f = open("/Volumes/Macintosh HD - Data/Users/satie/Jupyter/h"+str(yer),'r')
        dcnt=0
        for data in f:
            if data[52].isdigit():
                mgn=int(data[52:54])
#                print(mgn)
                if mgn>0 and mgn<100:
                    eqMcl[mgn]=eqMcl[mgn]+1
        f.close()
        print("h"+str(yer)+" : ")
#        for mgn in range(10):
#            print(eqMcl[mgn*10:mgn*10+9])
for mth in ["01","10"]:
    url = "https://www.data.jma.go.jp/eqev/data/bulletin/data/hypo/h1997"+mth+".zip"
    with (
        urllib.request.urlopen(url) as res,
        io.BytesIO(res.read()) as bytes_io,
        zipfile.ZipFile(bytes_io) as zip,
    ):
        zip.extractall(extract_dir)
        f = open("/Volumes/Macintosh HD - Data/Users/satie/Jupyter/h1997"+mth,'r')
        dcnt=0
        for data in f:
            if data[52].isdigit():
                mgn=int(data[52:54])
#                print(mgn)
                if mgn>0 and mgn<100:
                    eqMcl[mgn]=eqMcl[mgn]+1
        f.close()
        print("h1997"+mth+" : ")
        for mgn in range(10):
            print(eqMcl[mgn*10:mgn*10+9])
for yer in range(1983,1997):
    url = "https://www.data.jma.go.jp/eqev/data/bulletin/data/hypo/h"+str(yer)+".zip"

    with (
        urllib.request.urlopen(url) as res,
        io.BytesIO(res.read()) as bytes_io,
        zipfile.ZipFile(bytes_io) as zip,
    ):
        zip.extractall(extract_dir)
        f = open("/Volumes/Macintosh HD - Data/Users/satie/Jupyter/h"+str(yer),'r')
        dcnt=0
        for data in f:
            if data[52].isdigit():
                mgn=int(data[52:54])
#                print(mgn)
                if mgn>0 and mgn<100:
                    eqMcl[mgn]=eqMcl[mgn]+1
        f.close()
        print("h"+str(yer)+" : ")
for yer in [1967,1961,1951,1919]:
    url = "https://www.data.jma.go.jp/eqev/data/bulletin/data/hypo/h"+str(yer)+".zip"

    with (
        urllib.request.urlopen(url) as res,
        io.BytesIO(res.read()) as bytes_io,
        zipfile.ZipFile(bytes_io) as zip,
    ):
        zip.extractall(extract_dir)
        f = open("/Volumes/Macintosh HD - Data/Users/satie/Jupyter/h"+str(yer),'r')
        dcnt=0
        for data in f:
            if data[52].isdigit():
                mgn=int(data[52:54])
#                print(mgn)
                if mgn>0 and mgn<100:
                    eqMcl[mgn]=eqMcl[mgn]+1
        f.close()
        print("h"+str(yer)+" : ")
#for mgn in range(10):
#    print(eqMcl[mgn*10:mgn*10+9])
print(eqMcl)
print("fin")
prnl=""
for mgn in range(1,100):
    if eqMcl[mgn]>0:
        prnl=prnl+"mgn["+str(mgn)+"]:="+str(round(math.log10(eqMcl[mgn]),5))+";"
#        print(math.log10(eqMcl[mgn]))
    if mgn%10==9:
        print(prnl)
        prnl=""

　ご覧になれば分かりますが繰り返し同じパターンの部分が多いです。それは例えば, 気象庁の掲載データの形式の問題で, 何故か月毎だったり, 何年か纏めてだったり, ある年だけ前半後半？みたくその年のデータが2つに分けて載せてあったりで, その都度その都度ファイル名の形式が変わるからで, その処理は面倒だから, もう具体的にほぼ同じコードを繰り返し書いているからです。

　データの出力は, MetaPostでそのまま利用できるように, 数えた上で常用対数にしたものを, MetaPostの数値(numeric)な配列形式で吐かせてます。

データ

　折角なので, 今回得られたデータの素なものを置いておこうと思いますが... 気象庁さんが文句言ったりするのだろうか？