【令和最新版】F-netの発震機構解をモーメントテンソルで描きたい（仮）

はじめに

F-netのメカニズム検索でGMT形式が選べなくなってから、はや数年。
もう終わりだよこの国、と思いましたが、chatGPTのおかげで取り合えず作成できたので共有します。

F-net形式とGMT形式の違い

座標系が違うのは知っていたんですが、大きさが違うことに最近気づいたんですよね。

座標系が違う。F-netは北,東,下、GMTは上,南,東。
単位が違う。F-netはNm、GMTは10^ exponent dynes-cm。　←気づかなかった。。

たとえば mxx = 0.0362 で Unit(Nm)=1e+15 のとき、実際のモーメント成分は 0.0362 × 1e+15 = 3.62e+13 Nm
これを dynes·cm に直すにはさらに 10^7 を掛けるので 3.62e+20 dyn·cm となる。

これをAIに指摘されて首を垂れる稲穂になりました。

課題

そうなると、GMTは10のべき乗表記で指数を表記してやらなければなりません。

確かに。つまり強引に合わせるので若干精度は落ちるということですね。

スクリプト

というわけで、スクリプトです。
元ファイルはF-netのカタログファイルを指定してやればOK.
検索システムの結果では動きません。

#!/usr/bin/env python3
import sys
import math
import numpy as np

def rotate_tensor_fnet_to_cmt(mxx, mxy, mxz, myy, myz, mzz):
    """
    F-net座標(N=x, E=y, D=z) --> GlobalCMT座標(r, t, f)への変換
      R = [[ 0,  0, -1 ],
           [-1,  0,  0 ],
           [ 0,  1,  0 ]]
      M' = R * M * R^T
    """
    R = np.array([[0.,  0., -1.],
                  [-1., 0.,  0.],
                  [0.,  1.,  0.]])
    # 元のテンソルM（対称行列）
    M = np.array([[mxx, mxy, mxz],
                  [mxy, myy, myz],
                  [mxz, myz, mzz]])
    # 回転後のテンソル M'
    M_prime = R @ M @ R.T
    # M'の6成分を展開
    #   M' = [[mrr, mrt, mrf],
    #         [mrt, mtt, mtf],
    #         [mrf, mtf, mff]]
    mrr = M_prime[0, 0]
    mtt = M_prime[1, 1]
    mff = M_prime[2, 2]
    mrt = M_prime[0, 1]
    mrf = M_prime[0, 2]
    mtf = M_prime[1, 2]
    return mrr, mtt, mff, mrt, mrf, mtf

def main():
    if len(sys.argv) < 2:
        sys.exit(f"Usage: {sys.argv[0]} input_file.txt")

    input_file = sys.argv[1]
    # 出力ファイル名を自動生成
    if input_file.lower().endswith('.txt'):
        output_file = input_file[:-4] + '_out.txt'
    else:
        output_file = input_file + '_out.txt'

    # 入力ファイルと出力ファイルをオープン
    with open(input_file, 'r', encoding='utf-8') as f_in, open(output_file, 'w', encoding='utf-8') as f_out:
        first_line = True
        for line in f_in:
            line = line.strip()
            if not line:
                continue

            # ヘッダ行（"Origin_Time(UT)" などを含む）をスキップ
            if first_line and ("Origin_Time" in line or "Latitude" in line):
                first_line = False
                continue
            first_line = False

            tokens = line.split()
            if len(tokens) < 21:
                # 列数が足りなければスキップ
                continue

            # 列対応（0-based）
            #   0: Origin_Time(UT)
            #   1: Latitude(deg)
            #   2: Longitude(deg)
            #   3: JMA_Depth(km)
            #   4: JMA_Magnitude(Mj)
            #   5: Region_Name
            #   6: Strike
            #   7: Dip
            #   8: Rake
            #   9: Mo(Nm)
            #  10: MT_Depth(km)
            #  11: MT_Magnitude(Mw)
            #  12: Var._Red.
            #  13: mxx
            #  14: mxy
            #  15: mxz
            #  16: myy
            #  17: myz
            #  18: mzz
            #  19: Unit(Nm)
            #  20: Number_of_Stations
            origin_time = tokens[0]            # "2025/01/01,04:02:33.41"など
            lat = float(tokens[1])            # 緯度
            lon = float(tokens[2])            # 経度
            mt_depth = float(tokens[10])      # メカニズム推定深さ
            mxx_ = float(tokens[13])
            mxy_ = float(tokens[14])
            mxz_ = float(tokens[15])
            myy_ = float(tokens[16])
            myz_ = float(tokens[17])
            mzz_ = float(tokens[18])

            # Unit(Nm)
            unit_str = tokens[19]  # "1e+15" など
            try:
                unit_val = float(unit_str)
            except ValueError:
                unit_val = 1.0  # デフォルト

            # --------------------------------------------
            # (1) F-net公表値 × unit_val → Nm 実値
            # (2) Nm → dyn.cm (× 1e7)
            # (3) F-net座標 → GlobalCMT座標 へ変換
            # (4) exponent 調整
            # --------------------------------------------
            # (1)
            mxx_nm = mxx_ * unit_val
            mxy_nm = mxy_ * unit_val
            mxz_nm = mxz_ * unit_val
            myy_nm = myy_ * unit_val
            myz_nm = myz_ * unit_val
            mzz_nm = mzz_ * unit_val

            # (2)
            factor = 1.0e7
            mxx_dyn = mxx_nm * factor
            mxy_dyn = mxy_nm * factor
            mxz_dyn = mxz_nm * factor
            myy_dyn = myy_nm * factor
            myz_dyn = myz_nm * factor
            mzz_dyn = mzz_nm * factor

            # (3)
            mrr, mtt, mff, mrt, mrf, mtf = rotate_tensor_fnet_to_cmt(
                mxx_dyn, mxy_dyn, mxz_dyn, myy_dyn, myz_dyn, mzz_dyn
            )

            # (4) 6成分の絶対値最大を見て 10^exponent を決める
            abs_vals = [abs(mrr), abs(mtt), abs(mff), abs(mrt), abs(mrf), abs(mtf)]
            max_val = max(abs_vals) if abs_vals else 1e-30
            if max_val < 1e-30:
                exponent = 0
            else:
                exponent = int(math.floor(math.log10(max_val)))
            scale = 10.0 ** (-exponent)  # 各成分をこの値で乗じる

            mrr_ = mrr * scale
            mtt_ = mtt * scale
            mff_ = mff * scale
            mrt_ = mrt * scale
            mrf_ = mrf * scale
            mtf_ = mtf * scale

            # 出力: lon lat depth_km mrr mtt mff mrt mrf mtf exponent 0 0 text
            # text に origin_time を入れる(必要に応じて変更してください)
            out_line = (f"{lon:.4f} {lat:.4f} {mt_depth:.2f} "
                        f"{mrr_:.6g} {mtt_:.6g} {mff_:.6g} "
                        f"{mrt_:.6g} {mrf_:.6g} {mtf_:.6g} "
                        f"{exponent:d} 0 0 {origin_time}")
            f_out.write(out_line + "\n")

    print(f"Converted '{input_file}' -> '{output_file}'.")

if __name__ == "__main__":
    main()

結果

載せたら怒られるので載せません。

今回はここまで。