はじめに
MAXI(Monitor of All-sky X-ray Image) は、2009年7月23日に国際宇宙ステーション「きぼう」実験棟船外実験プラットフォームに設置されたX線観測装置です。宇宙ステーションが地球を周回する動きを利用して、全天のX線を観測することができます。
本記事では、MAXIの光度曲線データをPythonで処理し、matplotlibとplotlyを使用して可視化するスクリプトの使い方を詳しく解説します。
コードを読めばわかる方向けには、google colab での実行例を用意してますので、こちらも参考にしてください。
MAXIのライトカーブについては、
のページから好きな天体を選び、
画像を選択すると、個々の天体のページに飛びます。個々の天体のページの下の方にある、1 day LCなどを選択すると、1日ごとの light curve の fits ファイルが取得できます。
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
という名前で保存します。
あるいは、
からダウンロードしてください。
#!/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時間であり、重心時刻を得るにはさらに補正が必要です。