全天X線監視装置(MAXI)を用いた光度曲線のプロット方法：matplotlibとplotlyを用いたPythonスクリプト解説

# はじめに

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

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

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

https://colab.research.google.com/drive/1axbmskPdxg5CrGwBfTHt0wgT3PcaiLZs?usp=sharing


MAXIのライトカーブについては、

http://maxi.riken.jp/top/lc_bh.html

のページから好きな天体を選び、

![スクリーンショット 2024-11-29 19.20.50.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/179652/3686291f-b0de-fc8f-435f-2ff6d031edf5.png)

画像を選択すると、個々の天体のページに飛びます。個々の天体のページの下の方にある、1 day LCなどを選択すると、１日ごとの light curve の fits ファイルが取得できます。

![スクリーンショット 2024-11-29 19.22.39.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/179652/76de5abd-30cb-64d8-0af1-0b8858b9d8bb.png)

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

https://www.darts.isas.jaxa.jp/astro/maxi/analysis/index.html

や、

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

http://maxi.riken.jp/mxondem/

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

Swift/BAT については、

https://qiita.com/yamadasuzaku/items/dc1f1a999cabed46fd63

を参考にしてください。

## 目的

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

## 前提条件

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

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

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

```bash
pip install numpy matplotlib plotly astropy requests
```

## スクリプトの概要

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

1. **データのダウンロード**
指定したURLからMAXIのデータファイルをダウンロードします。

2. **データの読み込みとフィルタリング**
データを読み込み、相対誤差のしきい値に基づいてデータをフィルタリングします。

3. **ハードネス比の計算**
フィルタリングされたデータから、異なるエネルギーバンド間のフラックス比を計算します。

4. **データの可視化**
matplotlibまたはplotlyを使用して、光度曲線やハードネス比をプロットします。

5. **表示範囲とハイライト期間の指定**
ユーザーが指定した期間で表示範囲を制限し、特定の期間をハイライトします。

---

## スクリプトの使い方

### 1. スクリプトの準備

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

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

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

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

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

あるいは、

https://github.com/yamadasuzaku/rksysoft/blob/main/astro/maxi/maxi_download_data.py

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

```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()
```

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

``` bash:
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` のコマンドライン引数の説明

https://github.com/yamadasuzaku/rksysoft/blob/main/astro/maxi/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` を使用する場合

```bash
python maxi_plot_lc.py J1958+352_g_lc_1day_all.dat
```

#### 4.2 plotlyを使用する場合

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

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

``` bash:
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
```

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

``` bash:
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`として保存されます。

![J1958+352_g_lc_1day_all_matplotlib.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/179652/083ed731-8ae1-58a2-506b-5eae574e168a.png)

#### 5.2 plotlyの場合

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


![スクリーンショット 2024-10-07 1.18.57.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/179652/815fd030-947f-d7fe-e649-8b08c75ca055.png)


---

## スクリプトの詳細解説

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

```python
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. データの読み込みとフィルタリング

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

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

### 3. ハードネス比の計算

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

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

### 4. データの可視化

#### 4.1 matplotlibを使用する場合

```python
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を使用する場合

```python
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` を用いたシンプルな台形補完が初手としてはおすすめです。

``` python:
    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でもこのスクリプトを実行できます。

https://colab.research.google.com/drive/1m-44bPAgTCkXkOw7_fDFAc4kvAaoGcbr?usp=sharing

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

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

```python
# 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 ペリオドグラムを使うことができます。

https://qiita.com/yamadasuzaku/items/f1dddccb29561302a62a

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

### 光度曲線データの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](http://maxi.riken.jp) のオンデマンドタブをクリック）を使用してください。光度曲線ビンに0.02（日）（つまり<90分）を指定すると、MJD列でスキャンの中心時刻が得られます。
- いずれの場合も、MJDはISS時間であり、重心時刻を得るにはさらに補正が必要です。

### 参考文献

- [MAXI公式サイト](http://maxi.riken.jp/)
- [Astropy Documentation](https://docs.astropy.org/en/stable/)
- [Plotly Documentation](https://plotly.com/python/)
- https://qiita.com/yamadasuzaku/items/dc1f1a999cabed46fd63