気象庁地震カタログを用いて経緯度と震央地名の対応表を作る

この投稿は、[地震情報アプリ界隈 Advent Calendar 2021](https://adventar.org/calendars/6729) の25日目の記事です。ご参加頂いた方々、ありがとうございました。
昨日は、はるさんの [テロップソフトのすゝめ](https://qiita.com/kanto_haru_dev/items/6da5a6d85912df1c3e70) でした。

# はじめに
気象庁が発表する地震情報では、地震の発生した場所を地名で表現します。例えば「茨城県南部で地震」「宮城県沖で地震」などという言葉をよく聞きますよね。この「茨城県南部」「宮城県沖」といった地名は**震央地名**と呼ばれ、気象情報などとは別に地震情報ならではの地域区分がなされています[^1]。

どの地域がどの震央地名で表されるかは、気象庁HPで公開されています。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/565515/313df382-e72b-a20b-3309-fa4f7dbc8c33.png)
出典： [気象庁 | 地震情報で用いる震央地名（日本全体図）](https://www.data.jma.go.jp/svd/eqev/data/joho/region/index.html)

例えば東北や関東の周辺を見ると、震央地名の区分を示す赤い線が0.1度刻みで折れ曲がっていることが分かります。これは、震央を経緯度で表す際に、速報値の最小単位が0.1度であることと対応していると思われます。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/565515/99b3ac7e-6551-3eaa-10ab-5daf8d31bb0d.png)
出典： [気象庁 | 地震情報で用いる震央地名（東北・関東・北陸・東海地方）](https://www.data.jma.go.jp/svd/eqev/data/joho/region/index.html)

しかし、仙台管区気象台のページでは、区分を示す線が折れ曲がっていません。速報値では0.1度刻みですが、暫定値ではもっと細かい（例えば0.1分）刻み幅ですから、実際はこのように連続した線をもって震央地名が設定されていそうです。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/565515/ba7ab89b-956b-5476-7bb5-0505309205e0.png)
出典： [地震情報で発表する東北地方と周辺海域の震央地名（仙台管区気象台）](https://www.data.jma.go.jp/sendai/knowledge/earthquake/sinnou.html)

では、ある経緯度（震央）がどの震央地名に属しているのかを知りたいとき、どのようなデータがあれば良いでしょうか？
今回は**気象庁が実際に発表した地震情報をすべて解析**し、経緯度と震央地名の対応表を**ゴリ押しで**作成します。

# データの取得と対応表の作成
今回は、地震月報（カタログ編）の震度データを使用し、0.1度刻みの経緯度に対して震央地名の候補を返すようなテーブルを作ります。
なお、地震月報（カタログ編）の震度データには、震央地名が日本語で収録されています。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/565515/64d03730-ad52-9bb2-c708-149cf297c3fd.png)
（[気象庁｜地震月報(カタログ編)](https://www.data.jma.go.jp/svd/eqev/data/bulletin/index.html)より）

1919年から2019年まである全ての震度データをダウンロードし、Pythonで対応表を作成します。

```python:ll2epiname.py
from decimal import Decimal, ROUND_HALF_UP
import re

d = {}
for i in range(1919, 2020):
    with open(f"i{str(i)}.dat", encoding="shift_jis") as f:
        for line in[x.rstrip() for x in f.readlines()]:
            # 震源レコードかの判定
            if re.match(r'[A-Z]', line) is None:
                continue

            # 不明データかの判定
            if "不明" in line:
                continue

            # 緯度
            lat_d = line[21:24].strip()
            lat_m = line[24:28].strip()
            lat_i = float(lat_d)
            if lat_m:
                lat_i += float(lat_m) / 10 ** (len(lat_m) - 2) / 60
            lat = str(Decimal(lat_i).quantize(Decimal('0.1'), rounding=ROUND_HALF_UP))

            # 経度
            lon_d = line[32:36].strip()
            lon_m = line[36:40].strip()
            lon_i = float(lon_d)
            if lon_m:
                lon_i += float(lon_m) / 10 ** (len(lon_m) - 2) / 60
            lon = str(Decimal(lon_i).quantize(Decimal('0.1'), rounding=ROUND_HALF_UP))

            # 震央地名
            name = re.match(r"(\S+)", line[68:]).groups()[0]

            t = (lat, lon)
            if t in d:
                d[t].add(name)
            else:
                d[t] = {name}

# 緯度の昇順、経度の昇順でソート
l = sorted(d.items(), key=lambda x: (float(x[0][0]), float(x[0][1])))

# 整形
l = [f"{x[0][0]},{x[0][1]},{len(x[1])},{' '.join(x[1])}" for x in l]

# 出力
with open("ll2epiname.csv", "w") as f:
    f.write("\n".join(l))

```


有効な震度データ**合計119,933件**から、経緯度の組合せ**合計10,638件**について、経緯度と震央地名の対応表を作成できました。
対応表を一部抜粋すると以下のようになっています。


```plaintext:ll2epiname.csv
北緯,東経,震央地名の候補の個数,震央地名の候補（スペース区切り）
35.8,139.0,2,東京都多摩西部 山梨県東部・富士五湖
35.8,139.1,1,東京都多摩西部
35.8,139.2,1,東京都多摩西部
35.8,139.3,3,東京都多摩東部 東京都多摩西部 埼玉県南部
35.8,139.4,2,東京都多摩東部 埼玉県南部
35.8,139.5,2,東京都多摩東部 埼玉県南部
35.8,139.6,3,東京都多摩東部 埼玉県南部 東京都２３区
35.8,139.7,2,埼玉県南部 東京都２３区
35.8,139.8,2,東京都２３区 埼玉県南部
35.8,139.9,3,東京都２３区 千葉県北西部 埼玉県南部
35.8,140.0,1,千葉県北西部
35.8,140.1,1,千葉県北西部
35.8,140.2,1,千葉県北西部
35.8,140.3,1,千葉県北西部
35.8,140.4,2,千葉県北東部 千葉県北西部
35.8,140.5,1,千葉県北東部
35.8,140.6,1,千葉県北東部
35.8,140.7,2,茨城県南部 千葉県北東部
35.8,140.8,3,茨城県南部 千葉県北東部 千葉県東方沖
35.8,140.9,1,千葉県東方沖
35.8,141.0,1,千葉県東方沖
```

対応関係を調べられた地域を地図に表すと、このようになりました。ここで、日本周辺の矩形は、その矩形が属する経緯度（0.1度メッシュ）に、震央地名データが存在することを示します。矩形の色は震央地名の候補の個数に対応していて、白が1候補、黄が2候補、赤が3候補、茶色が4候補を表します。
この図を見ると、震央地名の境界にあたる黒線の周辺や、北海道〜東北の沖合で、震央地名の混在がある、すなわち震央地名の境界があることが分かります。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/565515/22c82499-b8f1-f1dc-f1f7-93d2443a00a6.png)

関東に注目すると、栃木県南部、群馬県南部、埼玉県北部、茨城県南部の4震央地名が混在する経緯度メッシュで、4候補を示す茶色となっていますから、解析はだいたい合っていそうです。
![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/565515/6c5a8683-8cf0-97a6-f87d-211a8ea99c7e.png)

なお、震央地名の境界線にあたるのに色が白い（候補が1つしかない）地域については、震度1以上の地震が地域をまたぐように発生していないことによるデータの不足を意味します。

# おわりに
今回は、実際に発生した地震のデータを用いて、経緯度と震央地名の対応表を作成しました。
重要なのは、特定の経緯度がどの震央地名に属するかよりも、震央地名の境界がどうなっているかですね。

もしかしたらお気づきの方もいるかもしれませんが、内陸については震央地名の境界がどこにあるかという正解が既に[気象庁GISデータ](https://www.data.jma.go.jp/developer/gis.html)という形で存在しますから、これを使えば空間解像度をより上げられることが出来るかもしれません。
海域については、序盤で紹介した震央地名の地図を見ると、境界はある程度の長さの直線で定義されています。そこで、線の傾きと切片の2変数をOpenCVの画像認識や何かで求めてあげれば、空間解像度の改善になるのかもしれません。難しいですが……。

今回解析した結果は、GitHubに公開しますので、使いたい方がいましたらぜひお使いください。
[https://github.com/compo031/ll2epiname](https://github.com/compo031/ll2epiname)

[^1]: 震央地名には地震・火山月報（防災編）によるものと、地震月報（カタログ編）によるものが存在します。ここでは、地震・火山月報（防災編）による震央地名を指します。