Go to Qiita Advent Calendar 2024 Top
LoginSignup
3
3

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?

More than 1 year has passed since last update.

記事投稿キャンペーン 「エンジニア×非エンジニアのコミュニケーション」
@after-a-rain

Google ColaboratoryでGMT(PyGMT)を使う #2

Last updated at Posted at 2023-10-29

<- Google ColaboratoryでGMT(PyGMT)を使う #1

はじめに

#1では，Google Colabratory上でPyGMTを動かしで日本地図の描画を行いました．#2では，気象庁の震源データを取得して，震源プロットを作成したいと思います．

やりかた

1.準備，2.気象庁から震源データのスクレイピング，3.気象庁の震源データの整理，4.PyGMTを用いた震源プロット の4つで構成されています．
2.と3.を行うことで，気象庁データをcsv化したものがGoogle Driveに保存され，4.ではそのcsvを読み込んでPyGMTでプロットします．そのため，1度本記事でcsvを作成し，かつcsvの内容を変更しなくて良い場合は，1.→4.の順で行って下さい．

1.準備

Google Colabratoryにアクセスし，#1同様にpygmtのインストールを行います．

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

#1では無かった

!mamba install ghostscript==9.54.0

が追加されていることに注意してください．後ほど解説します．

Google Driveをマウントします．
これにより，気象庁の震源データのダウンロードや作成した画像の保存をGoogle Driveに保存することができます．

# Google Driveに接続
from google.colab import drive
drive.mount('/content/drive')

スクレイピングとPyGMTの両方で使うライブラリのインポートと作業ディレクトリの指定を行います．

import os
DIR_PATH = "/content/drive/MyDrive/plotting_epicenter/" # 作業ディレクトリの指定
os.makedirs(DIR_PATH, exist_ok=True) # 上で指定した作業ディレクトリを作成 make dir(ectory)s

# データ処理関係
import pandas as pd
from datetime import datetime

# PyGMT
import pygmt

全体を通してファイル・ディレクトリの構造は以下の様になります．

content/
├ drive/
│ └ MyDrive/
│ 　 └ plotting_epicenter/(1.で作成)
│ 　 　 ├ data/(移動後)(2.で作成)
│ 　 　 │ ├ h1919.zip
│ 　 　 │ ├ h1951.zip
│ 　 　 │ └ ...
│ 　 　 ├ out2000-2022.csv (3.で作成)
│ 　 　 └ epi2000-2022.png (4.で作成)
└ data/ (保存ディレクトリ(移動前))(2.で作成)
　 ├ h1919.zip
　 ├ h1951.zip
　 └ ...

2.気象庁から震源データのスクレイピング

出典: 気象庁｜震源データ
上のURLにアクセスすると，各年(時々，月)ごとにダウンロードURLのリンクが一覧になっていることが分かります．
そのため，まず震源データHPの中からリンクを取得して，そのリンクにアクセスしてダウンロードし，最後にGoogle Driveに保存しています．
Google Driveに直接保存したところ10分くらいかかってしまったので，一旦Google Colab側で一時ディレクトリを作り，最後にDriveにコピーするようにしています．

# 震源データスクレイピング関係のライブラリ
from bs4 import BeautifulSoup # スクレイピング用
import urllib  # webページ操作用
import zipfile  # zip操作用
import shutil # move用


# 気象庁の震源データスクレイピング
# 5分くらいかかる

url = "https://www.data.jma.go.jp/eqev/data/bulletin/hypo.html" # 気象庁震源データHP

# (1)震源データHPのリンクのURL一覧を取得
response = urllib.request.urlopen(url)
html_content = response.read()
html = html_content.decode()

soup = BeautifulSoup(html, 'html.parser')
contents = soup.find_all('table')

url_list = [] # URL一覧を格納するリスト
for content in contents:
  get_a = content.find_all('a')

  for i in range(len(get_a)):
      try:
          link_ = get_a[i].get("href")
          link_ = f'https://www.data.jma.go.jp/eqev/data/bulletin/{link_[2:]}'
          url_list.append(link_)
      except:
        pass

df_title_url = pd.DataFrame({'URL':url_list})
df_notnull = df_title_url.dropna(how='any')
df_s = df_notnull.sort_values('URL', ascending=True) # 昇順にソート

temp_dir = '/content/data'
os.makedirs(temp_dir, exist_ok=True)

# (2)URL一覧の上からダウンロード
for url in df_s['URL']:
    try:
        res = urllib.request.urlopen(url)

    except urllib.request.HTTPError:
        print('Not found:', url)

    dst_path = os.path.join(temp_dir, os.path.basename(url))
    
    if not os.path.isfile(dst_path):
      data = res.read()
      with open(dst_path, mode="wb") as f:
        f.write(data)
      print('OK:', url)
      res.close()

    else:
      print('OK but already file exist:', url)

# (3)保存先を指定->作成． DIR_PATH/data に移動
save_dir = os.path.join(DIR_PATH, 'data')
os.makedirs(save_dir, exist_ok=True)

shutil.copytree(temp_dir, save_dir, dirs_exist_ok=True) # save_dir(DIR_PATH/data)にコピー
shutil.rmtree(temp_dir) # 一時ディレクトリ /content/data を削除

3.気象庁の震源データの整理

気象庁からの震源データは，こちらHPのフォーマットに従って記録されています．例えば，3.11の東北地方太平洋沖地震は

となっていますが，このままでは扱いにいくいので整理してCSVファイルに保存します．
ここでは，2000年から2022年3月までのデータをcsvにします．
output_filenameでファイル名を指定します．
なお，out2000-2022JU.csvのJUはJMAとUSGSのデータを使ったことを示しています．

# ダウンロードした震源レコードを指定の期間(e.g. 2000-2022)だけデータベース化
# Google Driveに csvで保存, pandasのDataframeも残る
# 2000-2022では，作業時間が4分半, csvの容量は300MBくらい

# 震源レコード用 (震度レコードと間違えないように注意)

# 関数: マグニチュードの変換
def convert_magnitude(input):
    input = input.strip()
    if input == "":
        return ""
    else:
        if input[0].isdigit():
            return float(input)/10
        elif input[0] == "-":
            return float(input)/10
        elif input[0] == "A":
            return -1 - float(input[1])/10
        elif input[0] == "B":
            return -2 - float(input[1])/10
        elif input[0] == "C":
            return -3 - float(input[1])/10

# 関数: 小数の変換
# 引数: 整数部分, 小数部分, 度分秒orNot
# 返り値: 変換後の小数(float) or データが空白の場合は""

def convert_float(intpart, fracpart, DMS=False):
  if not (intpart.replace(" ", "") == "" and fracpart.replace(" ", "") == ""):
    if not DMS:
      return round(float(f'{intpart.replace(" ", "") or "0"}.{fracpart.replace(" ", "0") or "0"}'), 6)
    else:
      fracpart_ = str(float(fracpart.replace(" ", "0") or "0")/6000)
      fracpart_ = fracpart_.split(".")[1]
      return round(float(f'{intpart.replace(" ", "") or "0"}.{fracpart_}'), 6)
  else:
    return ""

# データベース化する範囲を指定, en_yearも含まれる
st_year = 2000
en_year = 2022

# 出力ファイル名の指定, 保存先は DIR_PATH/(ファイル名).csv
output_filename = f"out{st_year}-{en_year}JU.csv"

# 震源レコードのデータベース化, ここから
f_list = [] # 保存用リスト
for year in range(st_year, en_year+1):
  print(year)
  for month in range(13):
      if not month:
        input_zippath = os.path.join(save_dir, f'h{year}.zip') # monthループが0の時はh(年).zipを指定
        if not os.path.isfile(input_zippath):
          continue
      else:
        input_zippath = os.path.join(save_dir, f'h{year}{month:02d}.zip') # monthループが1-12の時はh(年)(月).zipを指定
        if not os.path.isfile(input_zippath):
          continue

      input_filename = os.path.splitext(os.path.basename(input_zippath))[0] # h2000.zipの中にh2000というファイル名を指定
      with zipfile.ZipFile(input_zippath, 'r') as zf: # zip開く
        with zf.open(input_filename) as f:            # zipの中のファイルを開く
            for i, s_line in enumerate(f):            # 1行ずつ読み込み

              s_line = s_line.decode("utf-8")
              if s_line[0] in ['J', 'U']:
                  # レコード種別ヘッダ
                  if s_line[0] == 'J': # J: 気象庁による震源
                    etype = "JMA"
                  elif s_line[0] == 'U': # U: USGSが決定した震源
                    etype = "USGS"
                  elif s_line[0] == 'I': # I: その他の国際機関(ISC,IASPEIなど)による震源
                    etype = "Other"

                  edt = f'{s_line[1:15]}{s_line[15:17].replace(" ", ""):0<6}' # 地震発生の年月日時分秒ミリ秒
                  edt2 = datetime.strptime(edt, '%Y%m%d%H%M%S%f')             # edtをdatatime型に変換
                  etime_e = convert_float(s_line[17:19], s_line[19:21])       # 時刻の標準誤差

                  lat2 = convert_float(s_line[21:24], s_line[24:28], DMS=True) # 緯度
                  lon2 = convert_float(s_line[32:36], s_line[36:40], DMS=True) # 経度

                  lat_e = convert_float(s_line[28:30], s_line[30:32]) # 緯度の標準誤差
                  lon_e = convert_float(s_line[40:42], s_line[42:44]) # 経度の標準誤差

                  depth = convert_float(s_line[44:47], s_line[47:49]) #深さ
                  depth = round(depth, 2)

                  Mag1 = convert_magnitude(s_line[52:54]) # マグニチュード1
                  Mag1type = s_line[54]                   # マグニチュード1種別
                  Mag2 = convert_magnitude(s_line[55:57]) # マグニチュード2
                  Mag2type = s_line[57]                   # マグニチュード2種別

                  trt = s_line[58]                        # 使用走時表(Travel time)
                  ec_eva = s_line[59]                     # 震源評価
                  sp_info = s_line[60]                    # 震源補助情報
                  max_intensity= s_line[61]               # 最大震度

                  tmp = [etype, edt2.strftime('%Y-%m-%d'), edt2.strftime('%H:%M:%S.%f'), etime_e,
                          lat2, lat_e, lon2, lon_e, depth, Mag1, Mag1type, Mag2, Mag2type,
                          trt, ec_eva, sp_info, max_intensity]
                  f_list.append(tmp) # リストに追加

header = ["type", "date","time", "time_err", "lat","lat_err", "lon", "lon_err" , "depth","mag1","mag1t","mag2","mag2t", 'trt', 'ec_eva', 'info', 'max_intensity']
df = pd.DataFrame(f_list ,columns = header)
df.to_csv(os.path.join(DIR_PATH, output_filename))

ここで，マグニチュードは次のようなフォーマットになっているため，関数convert_magnitudeで数値(float型)に変換しています．

気象庁マグニチュード、気象庁ＣＭＴのモーメントマグニチュードまたはUSGS等が計算したマグニチュード
0未満の場合は以下のように表記する
-0.1: -1   -0.9: -9   -1.0: A0
-1.9: A9   -2.0: B0   -3.0: C0
マグニチュード1が求まらなかった場合は空白（半角スペース×２）

また，震源深さや緯度経度は小数点を省いた形式になっていますので，関数convert_floatで変換します．
例えば，震源深さは F5.2 なので，314.15km であれば 31415，3.14km であれば □□314
3.01km であれば □□3□1 と記載されます．(□は半角スペース)

緯度経度は分以下の扱いが少しややこしく，

25-28 F4.2 緯度(分) 震央の緯度(分) 震源固定の場合は小数点以下空白

とありますので
24.18分 であれば 2418と記載され，度に直す場合は，2418/6000 = 0.403 となります．

気象庁のデータからデータベース化する部分は，最初にzipのファイル名を検索します．
年ごとのforループの中に，月ごとのforループがネストする構造になっていて，月のループの
month == 0 のとき h2022.zip
month == 1 ~ 12 のとき h202201.zip
ようにファイル名を指定するようにしています．
(書きながら，先にglobで検索すればいいのではとか思っています...もっといい方法はありそうです．)

次に，zipを解答し中のファイルを開きます．なお，h2022.zipの中にはh2022というファイルが入っています．

その後，ファイルの中を1行ずつ読んでいき，何文字目から何文字目はこの情報といった具合に，欄ごとに変数に格納し，その変数をまとめてリストに追加しています．

それをst_yearとen_yearで指定した範囲だけ繰り返し，出来上がったリストをpandasのDataFrameにして，
df.to_csvでcsvに保存します．

(ここでも，書きながらnumpyかpandasで一気に読み込んで，文字数ごと(インデックスごと)にスライスなりしたほうが速いしPython的な書き方だな思っています...for回しても数分で終わるのでわざわざ書き換えるのもと思い今回はそのままにしています．)

CSVにすることで，色々使いまわせると思います．手元では400万レコード，300MB程度でしたので，Excelで読むことはしたくありませんが，何かと使えそうです．
Pythonで使う分には，これくらいのデータベースでは処理が重くなることはなさそうですが，もっと絞り込んでデータベース化，CSV化したいのであればst_yearとen_yearを変えることで可能です．

4.PyGMTを用いた震源プロット

データベースの準備ができましたので，実際にPyGMTでプロットしていきます．

まず，できたCSVを読み込みます(30秒程度)．

# csvからpandasに読み込む
df = pd.read_csv(os.path.join(DIR_PATH, 'out2000-2022JU.csv'))
df['date'] = pd.to_datetime(df['date']) # date列をdatetime型に変換

3.で作ったdfを読み込んでもいいですが，うまく型指定がされなかったこと，CSVを読み込むことで次回以降ここからでいいので汎用性が高まること，読み込みにそんなに時間かからないことから，都度読み込むようにしています．
date(年月日)はDateFrame型にならなかったのでここで変換しておきます．

続けて描画します．

# PyGMTを用いた震源プロット
# 2000年～2021年の日本全国のM3以上の地震をプロット

# データを整える
df_buf = df # 一時的なdfを作る
df_buf = df_buf[(df_buf['date'] >= datetime(2000,1,1)) & (df_buf['date'] <= datetime(2021, 12, 31))] # 指定範囲のデータを抽出
df_buf = df_buf[df_buf['mag1'] >= 3] # mag1列を3以上のデータを抽出


fig = pygmt.Figure()

# Define region of interest around Japan
region = [120, 150, 20 , 50]

title = f"Earthquakes in Japan 2000 - 2021"

# Load sample grid (5 arc-minutes global relief) in target area
# メッシュグリッド画像の解像度の指定(何度グリッドにするか, ここでは5分グリッド)
grid = pygmt.datasets.load_earth_relief(resolution="05m", region=region)

# Plot original grid
fig.basemap(region=region, projection="M12c", frame=["afg", f"+t{title}"])
fig.grdimage(grid=grid, cmap="geo", transparency = 50)
fig.coast(region=region, resolution='f', shorelines=True, area_thresh=50)

# 深さで色を指定するためのカラーマップを作成
pygmt.makecpt(
    cmap         = 'jet',        # カラーマップを選択
    series       = [1, 700, 10],  # min/max/increment
    background   = True,          # 値の上限・下限超過データの色を上限値・下限値と等しくする
    continuous   = True,
    reverse = True,
    transparency = 50,
)

fig.plot(
    x     = df_buf.lon,
    y     = df_buf.lat,
    size  = 0.003 * (2**df_buf['mag1']),
    style = 'c',
    cmap  = True, 
    fill  = df_buf.depth,

)
fig.colorbar(
    position = '+e',
    frame    = ['x+l"Depth [km]"'],
)

fig.show()
fig.savefig(DIR_PATH + f"epi2000-2021.png")

最初にデータを抽出します．今回は2000年(2000/1/1)から2021年(2021/12/31)の期間で，マグニチュード3以上を抽出しています．
マグニチュードの小さい地震の個数は指数関数的に増えることが経験的に知られており(グーテンベルク・リヒター則)，マグニチュード2などを20年規模で含めると描画時間が長くなるので注意が必要です．

続けて，#1同様，日本地図を描画し，fig.plotで震源をプロットします．
fig.grdimageで日本地図の描画しますが，transparency = 50を指定することで，淡くなり震源プロットが見やすくなります．
pygmt.makecptでカラーマップ(cmap)を作ります．深さ(km)に応じて色を指定したいので，範囲は1~700に指定し，データ量が多いので見やすくなるようtransparency = 50で半透明にしています．
fig.plotで震源をプロットします．sizeでマグニチュードに従った円の大きさになるように指定します．cmap=Trueでカラーマップを有効にし，fillで深さの値をもとに塗りつぶす様に指定します．
fig.colorbarでカラーバーを描画します．
fig.showでColab上に結果を表示させ，fig.savefigで最初に設定したDIR_PATHの中にepi2000-2021.pngを保存しています．

ところで， 1.準備 で，ghostscriptのバージョンを変える1行がありました．これは，本記事執筆時ではPyGMTと同時にインストールとされるghostscript 10.02.0ではtransparencyが使えなかったためです．

5.出力

epi2000-2021_edit.png

おわりに

PythonでGoogle Colabratory上でもGMTが使えるということは多くのメリットがあるなと思いました．ローカルの環境設定がいらないことやシェルスクリプトを書かなくて良いことは初学者にとってはハードルが下がるのではないでしょうか．また，Pythonではスクレイピングや数値計算なども行える点も便利だと思います．Pythonは多くのライブラリがあることで，やりたいことが少ない行数で書けることが優れた点の1つだと思っておりますが，今回の事例は特にその特徴が生きた様に感じています．

本記事は，記事投稿キャンペーン 「エンジニア×非エンジニアのコミュニケーション」に参加させていただきました．私の所属している組織で，GMTの使い方を学ぶために多くの時間を費やしている方を何名もお見受けしました．GMTは優れたツールですが，あくまで道具ですので何か勿体なく思ってしまいました．そんなタイミングで本記事の内容があまり周知されていないことを知り，執筆させていただきました．
プログラムが得意ではなかったり，初めて学ぶ地球科学分野に携わっている方の助けになれば幸いです．

本記事に関しまして，不備などがございましたらご指摘いただけると嬉しいです．また，ご質問等ございましたらお気軽にコメントいただければできる限り返信させていただきます．

謝辞

本記事の震源データは，気象庁より取得させていただきました．

参考文献



<- Google ColaboratoryでGMT(PyGMT)を使う #1

3
3
2

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
3
3

Delete article

Deleted articles cannot be recovered.

Draft of this article would be also deleted.

Are you sure you want to delete this article?