sunpy

SunPy - 太陽物理学解析

SunPy は Astropy エコシステムを基盤としながら、太陽物理学者に必要な機能を拡張して提供します。主な強みは、2D 画像データ用の Map オブジェクト、光度曲線用の TimeSeries、複数のリポジトリを横断して統一的にデータを検索する Fido インターフェースです。

使用する場面

• VSO、JSOC、その他のアーカイブから太陽データを検索してダウンロードする。
• 太陽画像（EUV、磁場図、白色光）を可視化する。
• 異なる太陽座標系間で座標を変換する（例：Stonyhurst から Helioprojective）。
• 太陽時系列データ（GOES X 線フラックス、黒点数）を解析する。
• 太陽の差動回転を補正する。
• 異なる観測装置のデータを重ね合わせる（例：AIA と HMI）。
• 太陽円盤のリムフィッティングと特徴追跡を行う。

リファレンスドキュメント

公式ドキュメント: https://docs.sunpy.org/
SunPy ギャラリー: https://docs.sunpy.org/en/stable/generated/gallery/index.html
検索パターン: sunpy.map.Map、sunpy.net.Fido、sunpy.coordinates.frames、sunpy.timeseries

基本原則

Map オブジェクト

2D 空間データの主要なデータ構造です。画像配列をメタデータ（FITS ヘッダー）および座標参照系（WCS）と組み合わせることで、物理的に正確なプロットや切り抜きを可能にします。

太陽座標

太陽物理学では、非直交座標系が多く使われます。SunPy は、太陽に対する地球の位置を考慮しながら、3D の太陽位置を2D の望遠鏡面に投影する計算を処理します。

Fido（Federated IDO）

データを検索するための統一インターフェースです。必要なデータ（時刻、観測装置、波長）を指定すると、Fido が異なるプロバイダーを自動的に検索します。

クイックリファレンス

インストール

pip install "sunpy[all]"

標準的なインポート

import sunpy.map
import sunpy.coordinates
from sunpy.net import Fido, attrs as a
import astropy.units as u
from astropy.coordinates import SkyCoord

基本パターン - 検索、ダウンロード、プロット

from sunpy.net import Fido, attrs as a
import sunpy.map

# 1. Search for data
result = Fido.search(a.Time("2011/06/07 06:30", "2011/06/07 06:35"),
                     a.Instrument.aia,
                     a.Wavelength(171*u.angstrom))

# 2. Download (returns list of file paths)
# downloaded_files = Fido.fetch(result)

# 3. Create and plot map
# aiayamap = sunpy.map.Map(downloaded_files[0])
# aiayamap.peek()

重要なルール

✅ 推奨事項

• sunpy.map.Map を使用する - 太陽画像を正しい座標メタデータとともに扱う、信頼できる唯一の方法です。
• 物理単位を使用する - 波長、距離、時刻には必ず astropy.units を使用します。
• 観測者情報を確認する - 太陽座標は観測者（通常は地球または衛星）の位置に依存します。マップのヘッダーに有効な観測者メタデータがあることを確認します。
• マップを正規化する - 画像を比較する際、露出時間の違いを考慮するために map.plot(norm=...) または aiamap.exposure_time を使用します。
• 座標変換をベクトル化する - SkyCoord 配列は個々の点よりもはるかに高速です。
• ダウンロードには Fido を使用する - 再試行、ローカルキャッシュ、複数のプロバイダーをシームレスに処理します。

❌ 非推奨事項

• 生の plt.imshow でプロットしない - これは World Coordinate System（WCS）を無視します。代わりに ax = plt.subplot(projection=my_map) を使用します。
• 時刻系を無視しない - Astropy と同様に、UTC と TAI の違いに注意します。
• 波長をハードコードしない - 検索で単位が考慮されるように、a.Wavelength(193 * u.AA) を使用します。
• 手動で切り抜かない - numpy 配列をスライスせず、座標メタデータを保持するために my_map.submap(bottom_left, top_right) を使用します。

アンチパターン（絶対に避ける）

# ❌ BAD: Manual cropping by pixel indices
# cropped_data = my_map.data[100:200, 300:400] # Coordinates lost!

# ✅ GOOD: Submap using SkyCoord
from astropy.coordinates import SkyCoord
bottom_left = SkyCoord(-500*u.arcsec, -500*u.arcsec, frame=my_map.coordinate_frame)
top_right = SkyCoord(500*u.arcsec, 500*u.arcsec, frame=my_map.coordinate_frame)
submap = my_map.submap(bottom_left, top_right=top_right)

# ❌ BAD: Hardcoding solar radius or distance
# distance = 1.496e11 # meters

# ✅ GOOD: Get from map or constants
from sunpy.coordinates import sun
dist = my_map.dsun # Distance to Sun at observation time

太陽マップ（sunpy.map）

操作と可視化

import sunpy.map
import matplotlib.pyplot as plt

# Load map
aia_map = sunpy.map.Map("aia_data.fits")

# Coordinate-aware plotting
fig = plt.figure()
ax = fig.add_subplot(projection=aia_map)
aia_map.plot(ax=ax)
aia_map.draw_limb(ax=ax, color='white')
aia_map.draw_grid(ax=ax, color='white', linestyle='dotted')

# Basic Map Math
# Note: Adding/subtracting maps is possible if they share the same WCS
diff_map = aia_map_t2 - aia_map_t1

太陽座標（sunpy.coordinates）

座標系の変換

from sunpy.coordinates import frames

# Helioprojective (Arcsec from Sun center)
coord = SkyCoord(100*u.arcsec, 200*u.arcsec, frame=frames.Helioprojective, 
                 obstime="2023-01-01", observer="earth")

# Convert to Heliographic Stonyhurst (Lat/Lon on the Sun)
hgs_coord = coord.transform_to(frames.HeliographicStonyhurst)
print(f"Lat: {hgs_coord.lat}, Lon: {hgs_coord.lon}")

# Rotating a coordinate with the Sun's differential rotation
from sunpy.physics.differential_rotation import solar_rotate_coordinate
new_coord = solar_rotate_coordinate(coord, time="2023-01-02")

データ検索（sunpy.net.Fido）

複雑なクエリ

from sunpy.net import Fido, attrs as a

# Search for GOES X-ray flares AND AIA images simultaneously
query = Fido.search(
    a.Time("2012-08-31 19:00", "2012-08-31 20:00"),
    (a.Instrument.goes & a.XRayRange("low")) | 
    (a.Instrument.aia & a.Wavelength(171 * u.AA))
)

print(query)
# files = Fido.fetch(query)

TimeSeries（sunpy.timeseries）

光度曲線の処理

import sunpy.timeseries as ts

# Load GOES data
goes = ts.TimeSeries("goes_data.nc")

# Plotting
goes.peek()

# Export to Pandas
df = goes.to_dataframe()
# Perform rolling mean in pandas
df_smooth = df.rolling('1min').mean()

太陽物理学のアルゴリズム

リムフィッティングと表面の計算

# Finding the solar limb in an image
from sunpy.map.maputils import all_coordinates_is_on_disk

# Mask pixels off-disk
mask = all_coordinates_is_on_disk(aia_map)
aia_map_masked = sunpy.map.Map(aia_map.data * mask, aia_map.meta)

実践的なワークフロー

1. AIA と HMI の重ね合わせ（複数観測装置による解析）

def plot_overlay(aia_file, hmi_file):
    aia = sunpy.map.Map(aia_file)
    hmi = sunpy.map.Map(hmi_file)
    
    # Re-project HMI to AIA's viewpoint
    hmi_reprojected = hmi.reproject_to(aia.wcs)
    
    fig = plt.figure()
    ax = fig.add_subplot(projection=aia)
    aia.plot(ax=ax)
    # Draw HMI contours (magnetic fields) over AIA (hot plasma)
    hmi_reprojected.draw_contours(ax=ax, levels=[-100, 100]*u.G, colors='white')

2. フレアのフラックス変化の測定

def get_flare_evolution(timeseries_file):
    goes = ts.TimeSeries(timeseries_file)
    # Filter for the X-ray flare window
    flare_data = goes.truncate("2017-09-06 11:50", "2017-09-06 12:30")
    
    # Get peak time and flux
    peak_time = flare_data.to_dataframe()['xrsa'].idxmax()
    return peak_time, flare_data.max()

3. 一連の画像に対する差動回転補正

def rotate_map_series(maps, reference_time):
    rotated_maps = []
    for m in maps:
        # Rotate each map to the reference time
        m_rot = m.rotate(reproject=True, order=3) # Basic rotation
        # For actual solar surface rotation, use differential_rotation modules
        rotated_maps.append(m_rot)
    return rotated_maps

パフォーマンスの最適化

プレビュー用のダウンサンプリング

大きな AIA 画像（4096x4096）はプロットに時間がかかります。

# Create a lower resolution map for fast display
aia_resampled = aia_map.resample([512, 512] * u.pix)

並列ダウンロード

Fido.fetch は複数の接続を使用して、ダウンロードを高速化できます。

files = Fido.fetch(result, path='./data/', max_conn=5)

よくある落とし穴と解決策

「観測者が見つからない」エラー

SunPy が衛星の位置を把握できないため、計算に失敗します。

# ❌ Problem: SkyCoord transformation fails
# ✅ Solution: Manually set observer if missing in FITS header
from sunpy.coordinates import get_body_heliographic_stonyhurst
obs_coord = get_body_heliographic_stonyhurst('earth', aia_map.date)
aia_map.meta['hgln_obs'] = obs_coord.lon.value
aia_map.meta['hglt_obs'] = obs_coord.lat.value
aia_map.meta['dsun_obs'] = obs_coord.radius.to(u.m).value

方位の問題（北が上ではありません）

太陽望遠鏡は回転します。

# ❌ Problem: Solar North is tilted in the raw image
# ✅ Solution: Rotate to align Solar North with the Y-axis
aia_rotated = aia_map.rotate() # This aligns the CROTA or PCi_j matrix

検索結果が空の場合

プロバイダーがダウンしていると、Fido が何も返さないことがあります。

# ✅ Solution: Always check result length
result = Fido.search(...)
if len(result) == 0:
    raise ValueError("No data found for the given criteria.")

SunPy は、私たちに最も近い恒星のダイナミクスを探究するための決定版ツールキットです。物理単位と太陽に特化した座標系を統合することで、堅牢で再現性の高い太陽地球物理学研究を可能にします。