astropy-astronomy

Astropy — 天文学与天体物理学工具包

概述

Astropy 是用于天文学的核心 Python 软件包，为天文研究提供必要功能：支持单位的计算、天球坐标转换、FITS 文件 I/O、宇宙学计算、精确时间处理、表格数据操作，以及 WCS 图像坐标映射。

适用场景

• 在不同天球坐标系（ICRS、Galactic、FK5、AltAz）之间转换
• 处理物理量和单位（Jy→mJy、秒差距→km、光谱等效关系）
• 读取、写入或操作 FITS 文件（图像和表格）
• 宇宙学计算（光度距离、回望时间、共动体积）
• 使用多种时间尺度（UTC、TAI、TT、TDB）和格式（JD、MJD、ISO）进行精确时间处理
• 根据天球位置对天文星表进行交叉匹配
• 在像素坐标和世界坐标之间进行 WCS 转换
• 对于通用表格数据：改用 pandas 或 polars
• 对于射电干涉测量：改用 CASA

前置条件

pip install astropy           # Core package
pip install astropy[all]      # With optional dependencies (regions, photutils, etc.)
pip install pytz              # For timezone conversions

快速入门

import astropy.units as u
from astropy.coordinates import SkyCoord
from astropy.time import Time
from astropy.io import fits
from astropy.table import Table
from astropy.cosmology import Planck18

# Units and quantities
distance = 100 * u.pc
print(f"{distance.to(u.km):.3e}")  # 3.086e+15 km

# Coordinates
coord = SkyCoord(ra=10.5*u.degree, dec=41.2*u.degree, frame='icrs')
print(f"Galactic: l={coord.galactic.l:.2f}, b={coord.galactic.b:.2f}")

# Cosmology
d_L = Planck18.luminosity_distance(z=1.0)
print(f"Luminosity distance at z=1: {d_L:.1f}")  # ~6780 Mpc

# Time
t = Time('2023-01-15 12:30:00')
print(f"JD: {t.jd:.6f}, MJD: {t.mjd:.6f}")

核心 API

1. 单位与物理量（astropy.units）

import astropy.units as u
import numpy as np

# Create quantities
distance = 10 * u.kpc
flux = 3.5e-15 * u.erg / u.s / u.cm**2
wavelength = 6563 * u.Angstrom

# Unit conversions
distance_ly = distance.to(u.lyr)
flux_jy = flux.to(u.Jy, equivalencies=u.spectral_density(wavelength))
print(f"Distance: {distance_ly:.2f}")

# Arithmetic with automatic unit tracking
velocity = 300 * u.km / u.s
time = 1 * u.Gyr
distance_traveled = (velocity * time).to(u.Mpc)
print(f"Distance traveled: {distance_traveled:.2f}")

# Equivalencies for domain-specific conversions
freq = wavelength.to(u.Hz, equivalencies=u.spectral())
energy = wavelength.to(u.eV, equivalencies=u.spectral())
parallax_dist = (0.1 * u.arcsec).to(u.pc, equivalencies=u.parallax())
print(f"Frequency: {freq:.3e}, Parallax distance: {parallax_dist:.1f}")

# Logarithmic units (magnitudes)
mag = -2.5 * u.mag
flux_ratio = mag.to(u.dimensionless_unscaled)

# Performance: pre-compute composite units
flux_unit = u.erg / u.s / u.cm**2 / u.Angstrom
fluxes = np.array([1e-15, 2e-15, 3e-15]) * flux_unit

# Custom units
bbl = u.def_unit('bbl', 158.987 * u.liter)

2. 坐标系（astropy.coordinates）

from astropy.coordinates import SkyCoord, EarthLocation, AltAz
from astropy.time import Time
import astropy.units as u

# Create coordinates (multiple formats)
c = SkyCoord(ra='05h23m34.5s', dec='-69d45m22s', frame='icrs')
c = SkyCoord(ra=10.5*u.degree, dec=41.2*u.degree)
c = SkyCoord(l=280*u.degree, b=-30*u.degree, frame='galactic')

# Transform between frames
c_gal = c.galactic
c_fk5 = c.fk5
print(f"Galactic: l={c_gal.l:.4f}, b={c_gal.b:.4f}")

# Observer-dependent AltAz (requires time + location)
location = EarthLocation(lat=40*u.deg, lon=-120*u.deg, height=1000*u.m)
obstime = Time('2023-06-15 23:00:00')
altaz = c.transform_to(AltAz(obstime=obstime, location=location))
print(f"Alt={altaz.alt:.2f}, Az={altaz.az:.2f}")

# Angular separation and matching
c1 = SkyCoord(ra=10*u.deg, dec=20*u.deg)
c2 = SkyCoord(ra=10.1*u.deg, dec=20.05*u.deg)
sep = c1.separation(c2)
print(f"Separation: {sep.arcsec:.2f} arcsec")

# Catalog matching
from astropy.coordinates import match_coordinates_sky
idx, sep, _ = coords1.match_to_catalog_sky(coords2)
matches = sep < 1 * u.arcsec

# Named object lookup
m31 = SkyCoord.from_name('M31')

# 3D coordinates with distance
c3d = SkyCoord(ra=10*u.deg, dec=20*u.deg, distance=50*u.kpc)
print(f"Cartesian: {c3d.cartesian}")

# Velocity information
c_vel = SkyCoord(ra=10*u.deg, dec=20*u.deg,
                 pm_ra_cosdec=5*u.mas/u.yr, pm_dec=-3*u.mas/u.yr,
                 radial_velocity=100*u.km/u.s)

3. FITS 文件处理（astropy.io.fits）

from astropy.io import fits
import numpy as np

# Read FITS file
with fits.open('observation.fits') as hdul:
    hdul.info()                    # Show HDU structure
    data = hdul[0].data            # Image data as NumPy array
    header = hdul[0].header        # Header as dict-like object

# Access header values
exptime = header['EXPTIME']
header['OBSERVER'] = 'Smith'       # Modify
header.add_history('Processed with astropy')

# Convenience functions
data = fits.getdata('image.fits')
header = fits.getheader('image.fits')
value = fits.getval('image.fits', 'EXPTIME')

# Create new FITS file
hdu_primary = fits.PrimaryHDU(data=np.zeros((100, 100)))
hdu_primary.header['OBJECT'] = 'M31'

# Multi-extension file
hdu_image = fits.ImageHDU(data=np.random.random((256, 256)), name='SCI')
hdu_table = fits.BinTableHDU.from_columns([
    fits.Column(name='ID', format='J', array=np.arange(100)),
    fits.Column(name='FLUX', format='E', array=np.random.random(100)),
    fits.Column(name='NAME', format='20A', array=['star']*100)
])
hdul = fits.HDUList([hdu_primary, hdu_image, hdu_table])
hdul.writeto('output.fits', overwrite=True)

# Large file handling with memory mapping
hdul = fits.open('huge.fits', memmap=True)
cutout = hdul[0].section[100:200, 100:200]  # Read only a slice

4. 表格操作（astropy.table）

from astropy.table import Table, QTable
import astropy.units as u
import numpy as np

# Create tables
t = Table({'ra': [10.0, 20.0, 30.0], 'dec': [41.0, 42.0, 43.0],
           'mag': [15.2, 16.1, 14.8]})

# Read from file (auto-detect format)
t = Table.read('catalog.fits')
t = Table.read('data.csv', format='csv')
t = Table.read('catalog.vot', format='votable')

# Unit-aware QTable
qt = QTable({'distance': [10, 20, 30] * u.kpc,
             'flux': [1e-15, 2e-15, 3e-15] * u.erg / u.s / u.cm**2})

# Filter, sort, column operations
bright = t[t['mag'] < 15.5]
t.sort('mag')
t['abs_mag'] = t['mag'] - 5 * np.log10(100)
print(f"Rows: {len(t)}, Columns: {t.colnames}")

# Joins and grouping
from astropy.table import join, vstack, hstack

# Database-style join
merged = join(t1, t2, keys='id', join_type='inner')

# Stack tables
combined = vstack([t1, t2, t3])       # Vertical (row-append)
combined = hstack([t_coords, t_phot])  # Horizontal (column-append)

# Group and aggregate
grouped = t.group_by('field')
stats = grouped.groups.aggregate(np.mean)

# Write
t.write('output.fits', format='fits', overwrite=True)
t.write('output.ecsv', format='ascii.ecsv')  # Preserves units + metadata

5. 时间处理（astropy.time）

from astropy.time import Time, TimeDelta
import astropy.units as u
import numpy as np

# Create from various formats
t = Time('2023-01-15 12:30:45', format='iso', scale='utc')
t = Time(2460000.0, format='jd')
t = Time(59945.0, format='mjd')
t = Time(1673785845.0, format='unix')

# Convert between formats and scales
print(f"ISO: {t.iso}")
print(f"JD: {t.jd}, MJD: {t.mjd}")
print(f"TAI: {t.tai.iso}")      # UTC → TAI (includes leap seconds)
print(f"TDB: {t.tdb.iso}")      # UTC → Barycentric Dynamical Time

# Time arithmetic
dt = TimeDelta(7, format='jd')
t_future = t + dt
t_future = t + 1 * u.hour
duration = Time('2024-01-01') - Time('2023-01-01')
print(f"Duration: {duration.jd:.1f} days")

# Array of times
times = Time('2023-01-01') + np.arange(365) * u.day

# Observing features
from astropy.coordinates import SkyCoord, EarthLocation

location = EarthLocation.of_site('Keck Observatory')
t = Time('2023-06-15 23:00:00', location=location)

# Sidereal time
lst = t.sidereal_time('apparent')
print(f"LST: {lst}")

# Barycentric correction
target = SkyCoord(ra='23h23m08.55s', dec='+18d24m59.3s')
ltt = t.light_travel_time(target, kind='barycentric')
t_bary = t.tdb + ltt
print(f"Barycentric correction: {ltt.sec:.3f} seconds")

6. 宇宙学计算（astropy.cosmology）

from astropy.cosmology import Planck18, FlatLambdaCDM
import astropy.units as u
import numpy as np

# Built-in cosmologies: Planck18, Planck15, Planck13, WMAP9, WMAP7
z = 1.5

# Distance calculations
d_L = Planck18.luminosity_distance(z)
d_A = Planck18.angular_diameter_distance(z)
d_C = Planck18.comoving_distance(z)
dm = Planck18.distmod(z)  # Distance modulus
print(f"d_L={d_L:.1f}, d_A={d_A:.1f}, d_C={d_C:.1f}")

# Time calculations
age = Planck18.age(z)
lookback = Planck18.lookback_time(z)
print(f"Age at z={z}: {age.to(u.Gyr):.2f}")
print(f"Lookback time: {lookback.to(u.Gyr):.2f}")

# Scale and volume
scale = Planck18.kpc_proper_per_arcmin(z)
vol = Planck18.comoving_volume(z)
print(f"Scale: {scale:.2f}")

# Inverse calculations — find z for given property
from astropy.cosmology import z_at_value

z_10gyr = z_at_value(Planck18.lookback_time, 10 * u.Gyr)
z_1gpc = z_at_value(Planck18.comoving_distance, 1 * u.Gpc)
print(f"z at lookback 10 Gyr: {z_10gyr:.4f}")

# Custom cosmology
cosmo = FlatLambdaCDM(H0=70, Om0=0.3, Tcmb0=2.725)
d_L_custom = cosmo.luminosity_distance(z=1.0)

# Array operations (all methods accept arrays)
z_array = np.linspace(0.1, 3.0, 100)
distances = Planck18.luminosity_distance(z_array)
print(f"Distance array shape: {distances.shape}")  # (100,)

7. WCS 与图像处理

from astropy.wcs import WCS
from astropy.io import fits
import astropy.units as u

# Read WCS from FITS
with fits.open('image.fits') as hdul:
    wcs = WCS(hdul[0].header)

# Pixel ↔ world transformations
world = wcs.pixel_to_world(100, 200)  # Returns SkyCoord
print(f"RA: {world.ra:.6f}, Dec: {world.dec:.6f}")

from astropy.coordinates import SkyCoord
coord = SkyCoord(ra=10.5*u.degree, dec=41.2*u.degree)
x, y = wcs.world_to_pixel(coord)

# WCS properties
print(f"Ref pixel: {wcs.wcs.crpix}")
print(f"Ref value: {wcs.wcs.crval}")
print(f"Pixel scale: {wcs.proj_plane_pixel_scales()}")
footprint = wcs.calc_footprint()  # Corner coordinates

# Image visualization
from astropy.visualization import simple_norm, ZScaleInterval, AsinhStretch, ImageNormalize
import matplotlib.pyplot as plt

data = fits.getdata('image.fits')
norm = simple_norm(data, 'sqrt', percent=99)
plt.imshow(data, norm=norm, cmap='gray', origin='lower')
plt.colorbar()

# Advanced normalization
interval = ZScaleInterval()
stretch = AsinhStretch()
norm = ImageNormalize(data, interval=interval, stretch=stretch)

# Sigma clipping for robust statistics
from astropy.stats import sigma_clipped_stats
mean, median, std = sigma_clipped_stats(data, sigma=3.0)
print(f"Background: {median:.2f} ± {std:.2f}")

核心概念

单位等效关系系统

Astropy 的 equivalencies 参数支持不满足量纲等价关系的特定领域转换：

等效关系	转换内容	示例
u.spectral()	波长 ↔ 频率 ↔ 能量	(500*u.nm).to(u.THz, u.spectral())
u.spectral_density(wav)	流量密度（Fλ ↔ Fν ↔ Jy）	flux.to(u.Jy, u.spectral_density(wav))
u.parallax()	视差角 ↔ 距离	(10*u.mas).to(u.pc, u.parallax())
u.doppler_optical(rest)	速度 ↔ 波长（光学）	vel.to(u.Angstrom, u.doppler_optical(rest))
u.brightness_temperature(freq)	流量 ↔ 温度	用于射电天文学

时间尺度

尺度	说明	适用场景
UTC	协调世界时（含闰秒）	默认；民用时间
TAI	国际原子时（UTC + 闰秒）	连续计时
TT	地球时（TAI + 32.184s）	地心计算
TDB	太阳系质心力学时	太阳系动力学、星历
UT1	地球自转角	恒星时、AltAz 转换

通过以下属性访问：t.utc、t.tai、t.tt、t.tdb、t.ut1

坐标系层级

• ICRS — 国际天球参考系（默认，约为 J2000）
• FK5 / FK4 — 历史赤道坐标系（FK4 需要春分点）
• Galactic — 银河坐标系（l、b）
• AltAz — 取决于观测者（需要 obstime + location）
• Ecliptic — 太阳系平面
• Galactocentric — 以星系为中心的笛卡尔坐标系

常见工作流

工作流 1：坐标转换流程

from astropy.coordinates import SkyCoord, EarthLocation, AltAz
from astropy.time import Time
import astropy.units as u

# Load catalog of sources
from astropy.table import Table
cat = Table.read('sources.fits')
coords = SkyCoord(ra=cat['RA']*u.degree, dec=cat['DEC']*u.degree)

# Transform to galactic
gal = coords.galactic
print(f"Galactic l range: {gal.l.min():.1f} to {gal.l.max():.1f}")

# Check observability (AltAz)
location = EarthLocation.of_site('Paranal Observatory')
obstime = Time('2023-06-15 23:00:00')
altaz = coords.transform_to(AltAz(obstime=obstime, location=location))
observable = altaz.alt > 30 * u.deg
print(f"Observable (alt>30°): {observable.sum()} of {len(coords)}")

工作流 2：FITS 图像分析

from astropy.io import fits
from astropy.wcs import WCS
from astropy.stats import sigma_clipped_stats
from astropy.visualization import simple_norm
import numpy as np

# Load image and WCS
with fits.open('science_image.fits') as hdul:
    data = hdul[0].data.astype(float)
    wcs = WCS(hdul[0].header)

# Background statistics
mean, median, std = sigma_clipped_stats(data, sigma=3.0)
print(f"Background: {median:.2f} ± {std:.2f}")

# Find bright pixels (simple threshold detection)
threshold = median + 5 * std
sources = np.where(data > threshold)
print(f"Pixels above 5σ: {len(sources[0])}")

# Convert pixel positions to sky coordinates
sky_coords = wcs.pixel_to_world(sources[1], sources[0])
print(f"RA range: {sky_coords.ra.min():.4f} to {sky_coords.ra.max():.4f}")

工作流 3：星表交叉匹配

from astropy.table import Table
from astropy.coordinates import SkyCoord
import astropy.units as u

# Read two catalogs
cat1 = Table.read('catalog1.fits')
cat2 = Table.read('catalog2.fits')

coords1 = SkyCoord(ra=cat1['RA']*u.degree, dec=cat1['DEC']*u.degree)
coords2 = SkyCoord(ra=cat2['RA']*u.degree, dec=cat2['DEC']*u.degree)

# Match
idx, sep, _ = coords1.match_to_catalog_sky(coords2)
max_sep = 1 * u.arcsec
matches = sep < max_sep

cat1_matched = cat1[matches]
cat2_matched = cat2[idx[matches]]
print(f"Matched: {matches.sum()} of {len(cat1)} (within {max_sep})")

关键参数

参数	模块	默认值	说明
frame	SkyCoord	'icrs'	坐标参考系
scale	Time	'utc'	时间尺度（utc、tai、tt、tdb、ut1）
format	Time	auto	时间格式（iso、jd、mjd、unix 等）
equivalencies	.to()	None	特定领域的单位转换规则
memmap	fits.open	True	对大型文件使用内存映射
join_type	join()	'inner'	联接类型（inner、outer、left、right）
sigma	sigma_clip	3.0	以标准差表示的裁剪阈值
stretch	simple_norm	'linear'	图像拉伸方式（linear、sqrt、log、asinh）
percent	simple_norm	100	归一化范围的百分位数

最佳实践

1. 始终附加单位 — 使用 Quantity 对象（例如 10 * u.kpc），以避免量纲错误。裸数值会悄无声息地产生错误结果。

2. 对 FITS 使用上下文管理器 — with fits.open(...) as hdul: 可确保正确关闭文件并清理内存映射。

3. 处理数组，而不是循环 — 所有 astropy 操作都接受数组。一次处理 10,000 个坐标，而不是逐个循环处理。

4. 明确指定时间尺度 — Time('2023-01-15', scale='utc') 可避免歧义。对高精度计时而言，UTC↔TDB 的差异很重要。

5. 使用 QTable 保存带单位的列 — QTable 会在输入/输出过程中保留单位；普通 Table 仅将单位作为元数据存储。

6. 使用 ECSV 确保往返保真度 — t.write('file.ecsv') 会保留单位、数据类型和元数据。CSV/FITS 会丢失部分元数据。

7. 反模式 — 使用错误的宇宙学模型：务必指定所用的宇宙学模型。不同模型在相同红移下会给出不同距离。Planck18 是当前标准。

8. 反模式 — 忽略 WCS 原点约定：astropy 使用从 0 开始的像素坐标。FITS 标准使用从 1 开始的坐标。请使用 wcs.pixel_to_world()（会自动处理此问题），而不要手动计算。

常见示例

示例：包含中微子的自定义宇宙学模型

from astropy.cosmology import FlatLambdaCDM
import astropy.units as u

cosmo = FlatLambdaCDM(
    H0=67.66, Om0=0.3111, Tcmb0=2.7255,
    Neff=3.046, m_nu=[0, 0, 0.06] * u.eV
)
print(f"Age of universe: {cosmo.age(0).to(u.Gyr):.3f}")

示例：模型拟合

from astropy.modeling import models, fitting
import numpy as np

# Generate noisy Gaussian data
x = np.linspace(0, 10, 100)
y = 10 * np.exp(-0.5 * ((x - 5) / 1.0)**2) + np.random.normal(0, 0.5, 100)

# Fit
fitter = fitting.LevMarLSQFitter()
model = models.Gaussian1D(amplitude=8, mean=4, stddev=1.5)
fitted = fitter(model, x, y)
print(f"Amplitude: {fitted.amplitude.value:.2f}")
print(f"Mean: {fitted.mean.value:.2f}")
print(f"Stddev: {fitted.stddev.value:.2f}")

示例：NDData 和 CCDData

from astropy.nddata import CCDData, StdDevUncertainty
import astropy.units as u
import numpy as np

# Create CCDData with uncertainty
data = np.random.random((256, 256))
uncertainty = StdDevUncertainty(np.sqrt(np.abs(data)))
ccd = CCDData(data, unit=u.adu, uncertainty=uncertainty,
              meta={'OBJECT': 'M31', 'EXPTIME': 300.0})

# Read/write
ccd.write('processed.fits', overwrite=True)
ccd2 = CCDData.read('processed.fits', unit=u.adu)

故障排除

问题	原因	解决方案
UnitConversionError	单位不兼容且未指定等价关系	添加 equivalencies=u.spectral() 或适当的等价关系
转换后坐标参考系错误	AltAz 缺少 obstime/location	同时提供二者：AltAz(obstime=t, location=loc)
ErfaWarning: dubious year	UT1 时间超出 1960–2040 年范围	对极端日期使用 scale='tt' 或 scale='tdb'
IERS 数据出现 FileNotFoundError	尚未下载闰秒表	运行 from astropy.utils.iers import IERS_Auto; IERS_Auto.open()
FITS 头出现 VerifyError	使用了非标准 FITS 关键字	使用 fits.open(f, ignore_missing_end=True) 或 hdul.verify('fix')
坐标转换缓慢	对单个坐标进行循环处理	使用数组形式的 SkyCoord：SkyCoord(ra=ra_array, dec=dec_array)
QTable 写入时丢失单位	使用了 CSV 格式	使用 ECSV 格式：qt.write('file.ecsv')
访问 FITS 扩展时出现 KeyError	扩展索引错误	使用 hdul.info() 查看结构；按名称访问：hdul['SCI']
重心修正不准确	时间尺度错误或缺少位置	使用 scale='utc'，并在 Time 对象上设置 location
ModelFit 未收敛	初始参数不佳	提供更接近的初始猜测值；尝试 SimplexLSQFitter

捆绑资源

• references/data_io_guide.md — 详细介绍 FITS 操作（标头、多扩展、二进制表、列格式代码、内存映射、远程访问）、表操作（创建、I/O 格式、联接、分组、索引、QTable、掩码数据、显示、性能）以及集合转换模式。由原始 fits.md 和 tables.md 合并而成。
• references/coordinates_time_cosmology.md — 完整介绍坐标系参考（所有坐标系、3D 坐标、自行、表示形式、目录匹配）、时间处理（所有格式、所有时间刻度、TimeDelta、恒星时、光行时、重心修正、精度）和宇宙学模型（所有内置模型、自定义模型、距离、体积、逆向计算、中微子效应）。由原始 coordinates.md、time.md 和 cosmology.md 合并而成。
• references/auxiliary_modules.md — WCS 详细操作、NDData/CCDData、建模框架（1D/2D 模型、拟合、复合模型）、图像可视化（拉伸、区间、归一化）、常数目录、卷积、稳健统计、SAMP 互操作性、数据下载工具。由原始 wcs_and_other_modules.md 和 units.md（等价关系详情）合并而成。

未迁移为单独文件：原先有 7 个参考文件。现已合并为 3 个主题参考文件，涵盖所有功能。units.md 中的等价关系内容分为 SKILL.md Key Concepts（摘要表）和 auxiliary_modules.md（详细代码）。

相关技能

• matplotlib-scientific-plotting — 用于 astropy 图像显示和绘图的可视化库
• zarr-python — 分块数组存储；与 astropy 配合用于大型天文数据集

参考资料

• 官方文档：https://docs.astropy.org/en/stable/
• 教程：https://learn.astropy.org/
• GitHub：https://github.com/astropy/astropy
• 坐标系参考：https://docs.astropy.org/en/stable/coordinates/