はじめに

MAXI（Monitor of All-sky X-ray Image) は、2009年7月23日に国際宇宙ステーション「きぼう」実験棟船外実験プラットフォームに設置されたX線観測装置です。宇宙ステーションが地球を周回する動きを利用して、全天のX線を観測することができます。

本記事では、MAXIの光度曲線データをPythonで処理し、matplotlibとplotlyを使用して可視化するスクリプトの使い方を詳しく解説します。

コードを読めばわかる方向けには、google colab での実行例を用意してますので、こちらも参考にしてください。

MAXIの quick なライトカーブのプロット方法について紹介する記事ですので、専門的な解析については、

や、

http://maxi.riken.jp/top/readme.html の MAXI on-demand data を読んで、

から on-demand でスペクトルなどの詳細なデータ解析にトライしてください。

Swift/BAT については、

を参考にしてください。

目的

MAXIのデータを取得し、Pythonで読み込む
データをフィルタリングし、必要な計算を行う
matplotlibまたはplotlyを使用してデータを可視化する
表示範囲やハイライト期間を指定して、特定の期間に焦点を当てる

前提条件

Python 3.xがインストールされていることと、最低限の Python の知識
必要なライブラリがインストールされていること

必要なライブラリのインストール

以下のライブラリを使用します。インストールされていない場合は、以下のコマンドを実行してください。

pip install numpy matplotlib plotly astropy requests

スクリプトの概要

紹介するスクリプトは以下の機能を持っています。

データのダウンロード
指定したURLからMAXIのデータファイルをダウンロードします。
データの読み込みとフィルタリング
データを読み込み、相対誤差のしきい値に基づいてデータをフィルタリングします。
ハードネス比の計算
フィルタリングされたデータから、異なるエネルギーバンド間のフラックス比を計算します。
データの可視化
matplotlibまたはplotlyを使用して、光度曲線やハードネス比をプロットします。
表示範囲とハイライト期間の指定
ユーザーが指定した期間で表示範囲を制限し、特定の期間をハイライトします。

スクリプトの使い方

1. スクリプトの準備

この記事では、Cyg X-1 を例として、プロット方法を紹介します。
まずは、
http://maxi.riken.jp/star_data/J1958+352/J1958+352.html
を確認して、注意事項や README に目を通しましょう。

一日平均のライトカーブをプロットする例を紹介します。

このファイルを python でダウンロードしましょう。

データをダウンロードするスクリプトをmaxi_download_data.pyという名前で保存します。

あるいは、

からダウンロードしてください。

maxi_download_data.py

#!/usr/bin/env python

import requests
import argparse
import os

# argparse設定
def get_args():
    parser = argparse.ArgumentParser(description="Download a file from a specified URL and save it locally.")
    parser.add_argument("url", help="URL of the file to download")
    parser.add_argument("--output", help="Output filename (default is derived from the URL)", default=None)
    return parser.parse_args()

# ファイルをダウンロードして保存する関数
def download_file(url, output_filename):
    try:
        # ファイルをダウンロード
        response = requests.get(url)

        # ダウンロードが成功したか確認
        if response.status_code == 200:
            # ファイルに保存
            with open(output_filename, "wb") as file:
                file.write(response.content)
            print(f"ファイルが '{output_filename}' として保存されました。")
        else:
            print("ファイルのダウンロードに失敗しました。ステータスコード:", response.status_code)
    except requests.exceptions.RequestException as e:
        print(f"エラーが発生しました: {e}")

def main():
    # 引数を取得
    args = get_args()

    # 出力ファイル名を設定
    if args.output:
        output_filename = args.output
    else:
        # URLからファイル名を取得
        output_filename = os.path.basename(args.url)

    # ファイルをダウンロード
    download_file(args.url, output_filename)

if __name__ == "__main__":
    main()

このスクリプトを使って、

python maxi_download_data.py http://maxi.riken.jp/star_data/J1958+352/J1958+352_g_lc_1day_all.dat

とするとファイルが保存されます。将来的に自動化をしたい方のために、このようにスクリプトでダウンロードできる方法を確認しておきます。

2. プロットするコード `maxi_plot_lc.py` のコマンドライン引数の説明

をダウンロードしておきます。

スクリプトはコマンドラインから以下の引数を受け取ります。

filename：データファイルのパスまたはURL
--plotter：使用するプロットライブラリ（matplotlibまたはplotly）、デフォルトはmatplotlib
--start_date：表示範囲の開始日時（ISOT形式、例：2023-10-15T00:00:00）
--end_date：表示範囲の終了日時（ISOT形式、例：2024-05-01T00:00:00）
--highlight：ハイライト期間（ISOT形式で開始と終了を指定、複数指定可能）

3. データの確認

プロットの前に、データのダウンロードが中途半端で終わっていることがあると、プログラムが動かないこともあるので、初めての方は、データの中身をぜひ確認してから進みましょう。

4. スクリプトの実行方法

4.1 `matplotlib` を使用する場合

python maxi_plot_lc.py J1958+352_g_lc_1day_all.dat

4.2 plotlyを使用する場合

python maxi_plot_lc.py J1958+352_g_lc_1day_all.dat --plotter plotly

ハイライトする時刻を指定したい場合

python maxi_plot_lc.py J1958+352_g_lc_1day_all.dat --highlight 2023-11-03T23:51:04 2023-12-05T14:01:04 --plotter plotly

開始時刻と終了時刻を指定したい場合

python maxi_plot_lc.py J1958+352_g_lc_1day_all.dat --start_date 2023-11-03T23:51:04 --end_date 2024-9-05T14:01:04 --plotter plotly

4.3 引数の詳細

--start_dateと--end_dateは省略可能です。指定しない場合、全期間が表示されます。
--highlightは複数回指定可能です。ハイライトしたい期間を開始日時と終了日時のペアで指定します。

5. 出力結果

5.1 matplotlibの場合

スクリプトを実行すると、光度曲線とハードネス比のプロットが表示され、同時に<ファイル名>_matplotlib.pngとして保存されます。

5.2 plotlyの場合

インタラクティブなプロットがブラウザで表示され、<ファイル名>_plotly.htmlとして保存されます。

スクリプトの詳細解説

1. データのダウンロード

def download_file(url, output_filename):
    try:
        response = requests.get(url)
        if response.status_code == 200:
            with open(output_filename, "wb") as file:
                file.write(response.content)
            print(f"ファイルが '{output_filename}' として保存されました。")
        else:
            print("ファイルのダウンロードに失敗しました。ステータスコード:", response.status_code)
    except requests.exceptions.RequestException as e:
        print(f"エラーが発生しました: {e}")

指定したURLからデータをダウンロードし、ローカルに保存します。

2. データの読み込みとフィルタリング

def load_and_filter_data(filename, relative_error_threshold=10):
    data = np.loadtxt(filename)
    # データの抽出とフィルタリング処理

データをNumPy配列として読み込みます。
相対誤差が指定したしきい値を超えるデータを除去します。(デフォルトで相対誤差が10%以上のデータはカットするようにしています。天体によって最適化が必要です。)

3. ハードネス比の計算

def calculate_hardness_ratios(data):
    # 各エネルギーバンド間のフラックス比とその誤差を計算

異なるエネルギーバンド間のフラックスの比率（ハードネス比）を計算します。
誤差伝播を用いて、ハードネス比の誤差も計算します。

4. データの可視化

4.1 matplotlibを使用する場合

def create_matplotlib_subplots(dates, flux_data, ratio_data, base_filename, start_date=None, end_date=None, highlight_periods=None):
    # matplotlibを用いたプロットの作成

3つのサブプロットを作成し、光度曲線とハードネス比を表示します。
表示範囲やハイライト期間を適用します。

4.2 plotlyを使用する場合

def create_plotly_subplots(dates, flux_data, ratio_data, base_filename, start_date=None, end_date=None, highlight_periods=None):
    # plotlyを用いたプロットの作成

インタラクティブなプロットを作成します。
ハイライト期間はオレンジ色で表示されます。

5 ハイライト部の値の内挿の例

ハイライトした時刻の対して、np.interp を用いた台形補完で値を出すコードの例もつけています。スパースなデータに対しては、np.interp を用いたシンプルな台形補完が初手としてはおすすめです。

    if highlight_periods:
        # ハイライト部分の値を内挿で推定する例
        # datetimeのリストを数値に変換
        dates_num = np.array([date.timestamp() for date in flux_data["dates"]])
        # highlight_periodsをtimestampに変換
        highlight_times = np.array([date.timestamp() for period in highlight_periods for date in period])
        # numpy.interpを用いてhighlight_periodsの時間で内挿したflux_2_20keVの値を計算
        highlight_flux_2_20keV = np.interp(highlight_times, dates_num, flux_data["flux_2_20keV"])
        # 内挿結果の表示
        for i, period in enumerate(highlight_periods):
            print(f"Period {period[0]} to {period[1]}:")
            print(f"Interpolated flux at {period[0]}: {highlight_flux_2_20keV[2*i]}")
            print(f"Interpolated flux at {period[1]}: {highlight_flux_2_20keV[2*i+1]}")

Google Colabでの実行

Google Colabでもこのスクリプトを実行できます。

に実行例を公開しています。

run_in_colab()関数の中身を書き換えて使ってください。
下記の例は、matplotlib で表示する例です。

# Google Colab向け: `argparse`を回避して直接値を指定できるようにする
def run_in_colab():

    # ファイルを指定
    filename = "J1958+352_g_lc_1day_all.dat"  # Colabにアップロードされたファイルを指定

    # プロッターを指定
    plot_library = "matplotlib" 
#    plot_library = "plotly"

    # 描画範囲の時刻を指定
#    start_date = None
#    end_date = None
    start_date = "2023-11-01T00:00:00"
    end_date = "2024-06-01T00:00:00"
    
    # ハイライトしたい時刻範囲を指定
    highlight_periods = [
        ("2023-11-03T23:51:04", "2023-11-05T14:01:04"),
        ("2024-04-07T16:55:04", "2024-04-10T13:41:04")
    ]
    main(filename, plot_library, start_date, end_date, highlight_periods)

if __name__ == "__main__":
    # Google Colab環境の確認
    if "google.colab" in sys.modules:
        print("Running in Google Colab environment. Using predefined settings.")
        run_in_colab()
    else:
        args = get_args()  # スタンドアローンの場合、引数をパース
        main(args.filename, args.plotter, args.start_date, args.end_date, args.highlight)

まとめ

MAXIのデータをPythonで処理し、matplotlibやplotlyを使用して光度曲線やハードネス比を可視化する方法を解説しました。宇宙の研究やデータ解析に役立てていただければ幸いです。

補足情報

Lomb-Scargle ペリオドグラムの作成

MAXIのデータは非一様な時系列データのため、単純なFFTによる解析は難しいです。非一様サンプルな時系列データに対しては、Lomb-Scargle ペリオドグラムを使うことができます。

などを参考にして試してみるのも一興です。

光度曲線データのREADMEの内容

データのフォーマット：

MJDcenter 2-20keV[ph/s/cm2] err 2-4keV err 4-10keV err 10-20keV err

MJDcenterは、1日光度曲線の場合はその日の中心時刻、1軌道光度曲線の場合は軌道の中心時刻です。
MJD列は、1日のデータでは24時間、6時間データでは6時間、軌道データでは1.5時間で等間隔に区切られています。
フラックスデータは、該当期間中に行われた1つまたは複数のスキャンの平均値です。
これはISS（国際宇宙ステーション）の時間であり、地心座標時や太陽系座標時ではありません。
1回のMAXIスキャンの継続時間は約60秒であり、92分ごとに行われます。
正確なスキャン時刻（スキャン停止時間の中心時刻）を知りたい場合は、v6L、v6m、およびv7Lページの「スキャンデータ」をダウンロードしてください。または、MAXIオンデマンド（http://maxi.riken.jp のオンデマンドタブをクリック）を使用してください。光度曲線ビンに0.02（日）（つまり<90分）を指定すると、MJD列でスキャンの中心時刻が得られます。
いずれの場合も、MJDはISS時間であり、重心時刻を得るにはさらに補正が必要です。

参考文献