astropy

Astropy - 天文学与天体物理学

Astropy 不只是一个库，更是社区标准。它将数值数组与物理单位和天球参考系关联起来，确保天文计算可复现且符合物理规律。

适用场景

• 在计算中处理物理单位和常数。
• 使用天球坐标（赤经、赤纬）并进行参考系转换。
• 读取/写入 FITS（灵活图像传输系统）、ASDF 和 VO（虚拟天文台）文件。
• 管理天文时间尺度（UTC、TAI、TDB、儒略日）。
• 世界坐标系统（WCS）映射（像素到天球坐标）。
• 宇宙学计算（距离、年龄、哈勃参数）。
• 处理带有单位和元数据的表格数据（astropy.table）。

参考文档

官方文档：https://docs.astropy.org/
Astropy 入门：https://learn.astropy.org/
搜索模式：astropy.units.Quantity、astropy.coordinates.SkyCoord、astropy.io.fits、astropy.time.Time

核心原则

单位和量

在 Astropy 中，所有内容都应是一个 Quantity——即与 Unit 配对的数字或数组。这样可以避免将米与英尺相加之类的错误。

天球坐标

坐标由 SkyCoord 对象表示，它们会自动处理球面三角学和参考系进动中的复杂计算。

带单位的表格

Pandas 非常适合处理数据框，但对于物理计算，astropy.table.QTable 更胜一筹，因为它在运算过程中会保留单位。

快速参考

安装

pip install astropy

标准导入

import numpy as np
from astropy import units as u
from astropy import constants as const
from astropy.coordinates import SkyCoord
from astropy.time import Time
from astropy.table import QTable
from astropy.io import fits
from astropy.wcs import WCS

基本模式 - 物理计算

# Calculate Schwarzschild radius of the Sun
mass_sun = 1.0 * u.M_sun
rs = (2 * const.G * mass_sun) / (const.c**2)

print(f"Radius: {rs.to(u.km):.2f}")
# Output: Radius: 2.95 km

关键规则

✅ 应当

• 使用 astropy.units - 从一开始就为原始数字附加单位。
• 优先使用 QTable - 使用 QTable 而非 Table，以便将列保留为 Quantity 对象。
• 使用 SkyCoord 进行计算 - 切勿手动计算角距离；请使用 c1.separation(c2)。
• 对坐标进行向量化处理 - SkyCoord 可以包含位置数组，比循环快得多。
• 指定时间尺度 - 创建 Time 对象时，始终定义时间尺度（utc、tdb）。
• 使用 with 处理 FITS - 始终使用上下文管理器，确保文件句柄关闭。

❌ 不应

• 过早剥离单位 - 除非确实需要将数据传递给非 Astropy 函数，否则不要使用 .value。
• 硬编码常数 - 不要使用 3e8 表示光速；请使用 const.c。
• 假定为 J2000 - 明确指定参考系（ICRS、FK5、Galactic）。
• 忽略 WCS - 不要假设像素坐标是线性的；请使用 WCS 对象进行映射。

反模式（切勿采用）

# ❌ BAD: Manual unit conversion (Error prone)
distance_pc = 10.0
distance_m = distance_pc * 3.086e16 

# ✅ GOOD: Automatic conversion
dist = 10.0 * u.pc
dist_m = dist.to(u.m)

# ❌ BAD: Stripping units to use with NumPy functions
# (Astropy Quantities are NumPy-compatible, no need to strip!)
import numpy as np
vals = [1, 2, 3] * u.m
res = np.mean(vals.value) # ❌ Unit info lost

# ✅ GOOD: Keep units
res = np.mean(vals) # ✅ Returns <Quantity 2. m>

# ❌ BAD: Manual RA/Dec math
# sep = np.sqrt((ra1-ra2)**2 + (dec1-dec2)**2) # ❌ Incorrect for spherical geometry

# ✅ GOOD: Proper separation
c1 = SkyCoord(ra=10*u.deg, dec=20*u.deg)
c2 = SkyCoord(ra=11*u.deg, dec=21*u.deg)
sep = c1.separation(c2)

单位和常数（astropy.units）

复杂转换

# Equivalencies (e.g., converting wavelength to frequency)
wav = 500 * u.nm
freq = wav.to(u.THz, equivalencies=u.spectral())

# Flux conversions
flux_density = 15 * u.Jy
obs_freq = 1.4 * u.GHz
brightness_temp = flux_density.to(u.K, equivalencies=u.brightness_temperature(obs_freq, beam_area=1*u.sr))

# Solar units
lum = 3.8e26 * u.W
print(lum.to(u.L_sun))

表格（astropy.table）

管理天文星表

# Create a QTable (Quantity Table)
tbl = QTable()
tbl['source_id'] = ['Star A', 'Star B', 'Star C']
tbl['flux'] = [1.2, 2.5, 0.8] * u.mJy
tbl['distance'] = [10, 25, 100] * u.pc

# Filtering
nearby = tbl[tbl['distance'] < 50 * u.pc]

# New column based on existing ones
tbl['luminosity'] = 4 * np.pi * tbl['distance']**2 * tbl['flux']

# Exporting
tbl.write('catalog.fits', overwrite=True)
tbl.write('catalog.csv', format='ascii.csv')

坐标（astropy.coordinates）

参考系转换

# Create coordinate from string
c = SkyCoord("05h35m17.3s", "-05d23m28s", frame='icrs')

# Convert to Galactic coordinates
gal = c.galactic
print(f"L: {gal.l:.2f}, B: {gal.b:.2f}")

# Transform to AltAz (Horizontal) for a specific location and time
from astropy.coordinates import EarthLocation, AltAz
loc = EarthLocation(lat=45*u.deg, lon=-110*u.deg, height=2000*u.m)
time = Time.now()

altaz = c.transform_to(AltAz(obstime=time, location=loc))
print(f"Altitude: {altaz.alt:.2f}")

时间（astropy.time）

精确计时

# Julian Date vs ISO
t = Time("2023-12-01T12:00:00", scale='utc')
print(t.jd)
print(t.mjd) # Modified Julian Date

# Time deltas
dt = 24 * u.hour
future_t = t + dt

# Siderial Time
print(t.sidereal_time('mean', longitude=-110*u.deg))

FITS 和 WCS（astropy.io.fits、astropy.wcs）

图像数据处理

# Reading FITS
with fits.open("galaxy.fits") as hdul:
    data = hdul[0].data
    header = hdul[0].header
    
    # Initialize WCS from header
    w = WCS(header)

# Convert Pixel to World (Sky)
pix_coords = [[100, 100], [200, 200]]
world_coords = w.pixel_to_world(pix_coords[0][0], pix_coords[0][1])

# Convert Sky to Pixel
target = SkyCoord(ra=150.0*u.deg, dec=2.0*u.deg, frame='icrs')
px, py = w.world_to_pixel(target)

实用工作流

1. 根据视差计算距离

def get_absolute_magnitude(mag, parallax_mas):
    """Calculate absolute magnitude given apparent mag and parallax in mas."""
    dist = (parallax_mas * u.mas).to(u.pc, equivalencies=u.parallax())
    # Absolute mag M = m - 5*log10(d/10pc)
    abs_mag = mag - 5 * np.log10(dist / (10 * u.pc))
    return abs_mag

# Star with mag 15 and 2 mas parallax
M = get_absolute_magnitude(15, 2)

2. 修正孔径测光单位

def flux_to_mag(flux, zero_point=25.0):
    """Convert raw flux to magnitudes with unit checking."""
    if not isinstance(flux, u.Quantity):
        raise ValueError("Flux must have units (e.g., u.electron / u.s)")
    
    # Assuming zero_point is defined for 1 electron/s
    mag = -2.5 * np.log10(flux / (1 * u.electron / u.s)) + zero_point
    return mag

3. 宇宙学距离演化

from astropy.cosmology import FlatLambdaCDM

# Define cosmology (Planck 2018 defaults)
cosmo = FlatLambdaCDM(H0=67.4, Om0=0.315)

z = np.linspace(0, 5, 100)
dist = cosmo.luminosity_distance(z) # Returns distance in Mpc

print(f"Lookback time at z=1: {cosmo.lookback_time(1).to(u.Gyr):.2f}")

性能优化

使用 astropy.utils.console 处理长时间任务

from astropy.utils.console import ProgressBar

with ProgressBar(len(large_catalog)) as bar:
    for source in large_catalog:
        # Complex calculation
        bar.update()

使用 C 引擎实现快速 I/O

# When reading huge CSV catalogs
tbl = QTable.read('massive_data.csv', format='ascii.fast_csv')

常见问题与解决方案

“剥离单位”错误

# ❌ Problem: Using a library that doesn't support Quantities
# solution = some_external_lib(data * u.m) # Crashes!

# ✅ Solution: Wrap the function or convert explicitly
data_to_pass = (data * u.m).to_value(u.cm) # Strip unit but ensure it's in CM

循环中的 SkyCoord 性能

# ❌ Problem: Creating SkyCoord objects inside a loop is very slow
for r, d in zip(ra_list, dec_list):
    c = SkyCoord(r, d, unit='deg')

# ✅ Solution: Create an array of coordinates
c_array = SkyCoord(ra=ra_list*u.deg, dec=dec_list*u.deg)

分解单位

# ❌ Problem: Result has messy units like 'm / km'
res = (10 * u.m) / (2 * u.km)

# ✅ Solution: Use .decompose() or .simplify()
print(res.decompose()) # Result: 0.005 (dimensionless)

Astropy 是天文学可复现研究的基石。它通过单位确保物理一致性，并规范坐标处理，从而避免过去困扰天文软件的隐蔽错误。