固定長データをPythonで頑張ってカンマ区切りにする(震源レコード編)

この記事は試行錯誤中の段階かつナンセンスなコードです。自己責任でお願いします。

本記事にて作成したデータセットの一部を配布しています。

データ配布
気象庁一元化震源データをカンマ区切りに変換したものを配布しています

目的

気象庁一元化震源データは未だに固定長データで運用されており、著者は大変ストレスに感じています。
そこで、今回はこの固定長データを、どうにかしてPython3でカンマ区切りにしてみます。

Fortranなら簡単にできるとか言わない。

震源レコードフォーマット

気象庁|震源レコードフォーマット
気象庁が提供するページです

面倒くさいですね。ゴリゴリカンマ区切りに変換してみます。

震源決定方法

A1 震源評価 震源を決定するにあたっての初期条件等
不明の場合は空白

1: 深さフリー
2: 深さ刻み条件(深さを一定の幅で変化させて計算)で最適解を求めた
3: 深さ固定等、人の判断による
4: Depth Phaseを用いた
5: S-Pを用いた
7: 参考(2016年3月まで)
8: 決定不能または不採用
9: 震源固定(最も早い相の情報を検測した観測点の緯度、経度)
M: Matched Filter法を用いた

面倒くさくなる諸悪の根源です。 これをもとに場合分けが必要になります。
震源固定の場合は後々面倒くさいことになります。深さ固定も同じみたいです。

固定長分割

pandas.read_fwf()の使い方 - sohatach's blog
www.oreilly.co.jp Think Stats 第2版を読んでいます。冒頭でサンプルデータを読み込んでいる import nsfg df = nsfg.ReadFemPreg() の中で実行されている pandas.read_fwf() の意味が公式ドキュメントを見てもさっぱり分からなかったので実際に動作...
pandas.read_fwf — pandas 1.2.2 documentation

ここら辺を読みました。pandas.read_fwfを使えば行けそう。
その結果がこちら。マグニチュードの分割だけアクロバティックなことをしています。

def kotei(hypo):
    # 固定長に分割する
    colspecs = [(0,1),(1,5),(5,7),(7,9),(9,11),(11,13),(13,17),(17,21),(21,24),(24,28),(28,32),(32,36),(36,40),(40,44),(44,49),(49,52),(52,53),(53,54), (54,55),(55,57),(57,58),(58,59),(59,60),(60,61),(61,62),(62,63),(63,64),(64,65),(65,68),(68,92),(92,95),(95,96)]
    names = ['recode','year','month','day','hour','minute','second','time_SE','latitude_deg','latitude_min','lat_SD','longitude_deg','longitude_min','lon_SD','depth','dep_SD','mag1_1', 'mag1_2','mag1_type','mag2','mag2_type','travel_time','hypo_judge','hypo_type','max_coefficient','scale_victim','tsunami_victim','region_main','region_sub','hypo_locale','obs_number','determination_way']
    df = pd.read_fwf(hypo, colspecs=colspecs, names=names)
    return df

おまじない

とりあえずFloat型に変換しておきます。

def convert(input):
    if input.latitude_deg.dtypes == 'O' :
        input['latitude_deg'] = input['latitude_deg'].str.replace(' ', '')
    if input.longitude_deg.dtypes == 'O' :
        input['longitude_deg'] = input['longitude_deg'].str.replace(' ', '')
    if input.depth.dtypes == 'O' :
        input['depth'] = input['depth'].str.replace(' ', '')
    if input.hypo_judge.dtypes == 'O' :
        input.replace({'hypo_judge' : {'M' : 10}}, inplace=True)
    input.hypo_judge = input.hypo_judge.astype("float64")
    input.latitude_deg = input.latitude_deg.astype("float64")
    input.latitude_min = input.latitude_min.astype("float64")
    input.longitude_deg = input.longitude_deg.astype("float64")
    input.longitude_min = input.longitude_min.astype("float64")
    input.depth = input.depth.astype("float64")
    return input

度分を直す

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

場合分けが必要なので面倒くさいですね。

def degree(input) :
    # 度分秒を度のみに直す
    # 震源固定の場合は小数点以下空白
    # 震源固定は震源評価9番。
    input.latitude_deg.mask(input.hypo_judge == 9, input.latitude_deg + (input.latitude_min / 60.0), inplace=True)
    input.longitude_deg.mask(input.hypo_judge == 9, input.longitude_deg + (input.longitude_min / 60.0), inplace=True)
    ## それ以外はF4.2なので、
    input.latitude_deg.where(input.hypo_judge == 9, input.latitude_deg + (input.latitude_min / 60.0 * 0.01), inplace=True)
    input.longitude_deg.where(input.hypo_judge == 9, input.longitude_deg + (input.longitude_min / 60.0 * 0.01), inplace=True)
    return input

マグニチュード1を直す

※便宜上、マグニチュード2は無視しています。

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

めちゃくちゃ複雑です。上のコードではマグニチュードを桁ごとに分けて処理していますが、自分でも正しいのか怪しいです。

def magnitude(input) :
    # マグニチュードを直す
    # 0未満の場合は以下のように表記する
    # -0.1: -1   -0.9: -9   -1.0: A0
    # -1.9: A9   -2.0: B0   -3.0: C0
    if input.mag1_1.dtypes == 'O' :
        input.replace({'mag1_1' : {'-' : -0}}, inplace=True)
        input.replace({'mag1_1' : {'A' : -1}}, inplace=True)
        input.replace({'mag1_1' : {'B' : -2}}, inplace=True)
        input.replace({'mag1_1' : {'C' : -3}}, inplace=True)
    input.mag1_1 = input.mag1_1.astype("float64")
    input.mag1_2 = input.mag1_2.astype("float64")
    input.mag1_1 = input["mag1_1"].where(input["mag1_1"] >= 0 , input["mag1_1"] - ( input["mag1_2"] / 10.0 ))
    input.mag1_1 = input["mag1_1"].where(input["mag1_1"] < 0 , input["mag1_1"] + ( input["mag1_2"] / 10.0 ))
    return input

発生時刻を直す

14-17 F4.2 秒 オリジンタイムの秒
震源固定の場合は小数点以下空白

これは簡単です。

def minute(input) :
    # 秒を直す。
    # 震源固定の場合は小数点以下空白
    # 震源固定は震源評価9番。
    input.second.where(input.hypo_judge == 9, input.second * 0.01, inplace=True )
    return input

深さを直す

F5.2 深さ(km) 深さフリーの条件で計算した時の震源の深さ(km)
I3, 2X 深さ固定または下記の刻みの条件で計算した時の震源の深さ(km)
10km(1926年~1960年、1967年~1982年)
20km(1961年~1966年)
1km(1983年~)
1982年以前の地震については適宜再調査され、深さフリーまたは1km刻みの震源に置き換えられる

ここは謎です。合っているのか分かりませんが、取り敢えず下記で動いています。

def depth(input) :
    # 深さを直す
    # 深さフリーならF5.2
    # 深さ固定ならI3
    # 深さフリーは震源評価1番
    input.depth.mask(input.hypo_judge == 1, input.depth * 0.01, inplace=True )
    input.depth.mask(input.hypo_judge == 10, input.depth * 0.01, inplace=True )
    return input

これ、
input.depth.where(input.hypo_judge == 3, input.depth * 0.01, inplace=True )
でもいいんですかね?よく分かりません。

整理

あとはお掃除です。
FAR FIELDは整理する際に厄介になるので削除しています。

    # FAR FIELDは削除する
    df = df[df.hypo_locale != "FAR FIELD"]
    # 出力
    df2 = df.round(4)
    df2 = df2.drop(columns=['latitude_min','longitude_min', 'mag1_2'])
    df2.to_csv(str(input) + ".csv", index=None)

とすれば、良い感じになると思います。

まとめ

こんな感じで組み合わせればうまく行くと思います。

def main() :
    input = args[1]
    # 固定長データに分割
    df = kotei(input)
    # 型を変更する
    df = convert(df)
    # 度分を直す
    df = degree(df)
    # マグニチュードを直す
    df = magnitude(df)
    # 秒を直す
    df = minute(df)
    # 深さを直す
    df = depth(df)
    # FAR FIELDは削除する
    df = df[df.hypo_locale != "FAR FIELD"]
    # 出力
    df2 = df.round(4)
    df2 = df2.drop(columns=['latitude_min','longitude_min', 'mag1_2'])
    #
    # df2.to_csv(str(input) + ".csv", index=None, header=None)
    df2.to_csv(str(input) + ".csv", index=None)

本当に面倒くさいのでカンマ区切りの震源レコードも作成してほしいです。。
USGSは作成しているのに…。

今日はここまで。

コメント

タイトルとURLをコピーしました