【Python】FITSファイル内のコラムの差分を計算して新しいコラムを追加する方法

はじめに

FITS形式のデータは天文学や物理学でよく利用されるファイル形式で、特に観測データの保存や解析に多用されています。今回は、Pythonを使ってFITSファイル内のあるコラムの差分を計算し、その結果を新しいコラムとして追加する方法を紹介します。

具体的には、astropy.io.fitsライブラリを使ってデータを操作し、指定したコラム（例えば、時刻データ）に対して1つ前の行との差分を計算し、FITSファイルに新しいコラムを追加するスクリプトを作成します。また、コマンドライン引数を用いて柔軟にコラム名を指定できるスクリプトも紹介します。

---

データの生成と差分の計算方法

まずは、FITSファイルを使ってデータを操作するための基礎として、時刻データとランダムなデータを生成し、その差分を計算する方法を紹介します。

from astropy.io import fits
import numpy as np

# ダミーデータの生成
n_rows = 100
time_data = np.linspace(0, 100, n_rows)  # 0から100までの時刻データ
value_data = np.random.rand(n_rows)  # ランダムなデータ

# 差分を計算（1行前の値との差）
time_diff = np.diff(time_data, prepend=time_data[0])  # 先頭に最初の値をプレースホルダーとして追加

# FITSファイル用のカラムを作成
cols = [
    fits.Column(name='TIME', format='E', array=time_data),
    fits.Column(name='VALUE', format='E', array=value_data),
    fits.Column(name='TIME_DIFF', format='E', array=time_diff)  # 差分列を追加
]

# カラム定義をテーブルHDUにパック
hdu = fits.BinTableHDU.from_columns(cols)

# FITSファイルに書き込み
output_filename = 'output_time_diff.fits'
hdu.writeto(output_filename, overwrite=True)

print(f'FITS file with TIME_DIFF created: {output_filename}')

解説

• np.linspaceで0から100までの時刻データを生成し、ランダムな値を持つデータ列を作成します。
• np.diffを使って、指定した列（ここでは TIME）の差分を計算します。prependを使うことで、最初の値をそのまま残し、差分計算を行います。
• 新しい列 TIME_DIFF を追加した上で、FITSファイルとして保存します。

このコードを実行することで、output_time_diff.fits というFITSファイルが作成され、元の時刻データ列 (TIME) に基づいて差分を計算した TIME_DIFF 列が追加されます。

---

コラムを追加するスクリプト

次に、コマンドライン引数を使って指定したFITSファイルとコラムに基づいて差分を計算し、新しいファイルに保存するPythonスクリプトを紹介します。このスクリプトは、コラム名をオプションで指定できるため、柔軟に利用できます。

diffadd.py

#!/usr/bin/env python 

import argparse
from astropy.io import fits
import numpy as np
import os

# コマンドライン引数の設定
def parse_args():
    parser = argparse.ArgumentParser(description='Generate a new column with differences in a FITS file.')
    parser.add_argument('fits_file', type=str, help='Input FITS file name')
    parser.add_argument('--cname', '-c', type=str, default='TIME', help='Column name to calculate differences (default: TIME)')
    return parser.parse_args()

# FITSファイルの差分を計算し、新しい列として追加
def calculate_diff_column(fits_file, cname):
    # FITSファイルを開く
    with fits.open(fits_file) as hdul:
        data = hdul[1].data
        # 指定したコラムから差分を計算
        column_data = data[cname]
        diff_data = np.diff(column_data, prepend=column_data[0])  # 最初の行はプレースホルダー
        
        # 新しいカラムを作成
        new_col_name = f"{cname}_DIFF"
        cols = hdul[1].columns
        new_col = fits.ColDefs([fits.Column(name=new_col_name, format='E', array=diff_data)])
        new_hdu = fits.BinTableHDU.from_columns(cols + new_col)
        
        # 新しいファイル名を作成
        output_filename = f"{cname}_DIFF_{os.path.basename(fits_file)}"
        new_hdu.writeto(output_filename, overwrite=True)
        print(f'FITS file with {new_col_name} created: {output_filename}')

# メイン関数
if __name__ == '__main__':
    args = parse_args()
    calculate_diff_column(args.fits_file, args.cname)

解説

• argparse を使ってコマンドライン引数を処理しています。fits_fileにはFITSファイル名を指定し、--cnameまたは-cオプションでコラム名を指定します（デフォルトは TIME）。
• FITSファイルを読み込み、指定されたコラム（cname）に基づいて差分を計算し、新しいコラムを追加します。新しいコラムは、cname_DIFF という名前になります。
• 出力ファイル名は cname_DIFF_元のファイル名.fits という形式になります。

---

使い方の例

例えば、VALUE というコラムの差分を計算したい場合、以下のようにコマンドを実行します。

%python diffadd.py output_time_diff.fits -c VALUE
FITS file with VALUE_DIFF created: VALUE_DIFF_output_time_diff.fits

このコマンドを実行すると、output_time_diff.fits ファイルの VALUE コラムに基づいて差分を計算し、新しいコラム VALUE_DIFF を追加した VALUE_DIFF_output_time_diff.fits というFITSファイルが生成されます。

---

まとめ

この記事では、Pythonを使ってFITSファイルのコラムに対して差分を計算し、新しいコラムを追加する方法を紹介しました。特に、argparseを使うことで、コマンドラインから柔軟にコラム名を指定できるスクリプトを作成し、実行例を示しました。

この方法を応用すれば、複数のFITSファイルや異なるコラムに対しても簡単に処理を行うことができます。天文データの解析やデータ処理の自動化に役立てていただければ幸いです。

---