sunpy

SunPy - 太阳物理分析

SunPy 基于 Astropy 生态系统构建，并扩展了其中的功能，以满足太阳物理学家的需求。它的核心优势包括用于二维成像数据的 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 对象

用于二维空间数据的主要数据结构。它将图像数组与元数据（FITS 标头）和坐标参考系统（WCS）结合起来，从而支持物理上准确的绘图和裁剪。

太阳坐标

太阳物理学会使用许多非笛卡尔坐标系。SunPy 会处理将三维太阳位置投影到二维望远镜平面上的计算，并考虑地球相对于太阳的位置。

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 绘图 - 这会忽略世界坐标系统（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 是探索离我们最近的恒星动态的权威工具包。通过整合物理单位和太阳专用坐标系，它能够支持稳健且可复现的日球层物理研究。