天文数据分析工具Skill astropy

天文数据分析工具Astropy是一个Python库,用于天文学数据处理、坐标转换、宇宙学计算和统计分析,支持FITS文件操作、天体坐标变换、单位管理和模型拟合,是天文学研究和量化金融中数据科学的核心工具。

数据分析 0 次安装 0 次浏览 更新于 3/16/2026

name: astropy description: “天文学工具包。FITS文件输入/输出、天体坐标转换、宇宙学计算、时间系统、世界坐标系统、单位、天文表格,用于天文数据分析和成像。”

Astropy

概述

Astropy 是天文社区标准的Python库,为天文数据分析和计算提供核心功能。本技能提供了全面的指南和工具,涵盖坐标系、文件输入/输出、单位和量值、时间系统、宇宙学、建模等多个方面的功能。

何时使用此技能

在以下情况下应使用此技能:

  • 处理FITS文件(读取、写入、检查、修改)
  • 执行天文参考坐标系之间的坐标变换
  • 计算宇宙学距离、年龄或其他量值
  • 处理天文时间系统和转换
  • 处理物理单位和维度分析
  • 使用专业列类型处理天文数据表格
  • 将模型拟合到天文数据
  • 转换像素坐标和世界坐标(WCS)
  • 对天文数据进行稳健的统计分析
  • 以适当的缩放和拉伸可视化天文图像

核心功能

1. FITS文件操作

FITS(灵活图像传输系统)是天文学中的标准文件格式。Astropy提供全面的FITS支持。

快速FITS检查: 使用包含的 scripts/fits_info.py 脚本进行快速文件检查:

python scripts/fits_info.py observation.fits
python scripts/fits_info.py observation.fits --detailed
python scripts/fits_info.py observation.fits --ext 1

常见的FITS工作流程

from astropy.io import fits

# 读取FITS文件
with fits.open('image.fits') as hdul:
    hdul.info()  # 显示结构
    data = hdul[0].data
    header = hdul[0].header

# 写入FITS文件
fits.writeto('output.fits', data, header, overwrite=True)

# 快速访问(对于多次操作效率较低)
data = fits.getdata('image.fits', ext=0)
header = fits.getheader('image.fits', ext=0)

# 更新特定的头关键字
fits.setval('image.fits', 'OBJECT', value='M31')

多扩展FITS

from astropy.io import fits

# 创建多扩展FITS
primary = fits.PrimaryHDU(primary_data)
image_ext = fits.ImageHDU(science_data, name='SCI')
error_ext = fits.ImageHDU(error_data, name='ERR')

hdul = fits.HDUList([primary, image_ext, error_ext])
hdul.writeto('multi_ext.fits', overwrite=True)

二进制表格

from astropy.io import fits

# 读取二进制表格
with fits.open('catalog.fits') as hdul:
    table_data = hdul[1].data
    ra = table_data['RA']
    dec = table_data['DEC']

# 更好:使用 astropy.table 进行表格操作(见第5节)

2. 坐标系和变换

Astropy支持约25个坐标系,并可进行无缝变换。

快速坐标转换: 使用包含的 scripts/coord_convert.py 脚本:

python scripts/coord_convert.py 10.68 41.27 --from icrs --to galactic
python scripts/coord_convert.py --file coords.txt --from icrs --to galactic --output sexagesimal

基本坐标操作

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

# 创建坐标(支持多种输入格式)
c = SkyCoord(ra=10.68*u.degree, dec=41.27*u.degree, frame='icrs')
c = SkyCoord('00:42:44.3 +41:16:09', unit=(u.hourangle, u.deg))
c = SkyCoord('00h42m44.3s +41d16m09s')

# 在坐标系之间变换
c_galactic = c.galactic
c_fk5 = c.fk5

print(f"银河系坐标: l={c_galactic.l.deg:.3f}, b={c_galactic.b.deg:.3f}")

处理坐标数组

import numpy as np
from astropy.coordinates import SkyCoord
import astropy.units as u

# 坐标数组
ra = np.array([10.1, 10.2, 10.3]) * u.degree
dec = np.array([40.1, 40.2, 40.3]) * u.degree
coords = SkyCoord(ra=ra, dec=dec, frame='icrs')

# 计算分离
sep = coords[0].separation(coords[1])
print(f"分离: {sep.to(u.arcmin)}")

# 位置角
pa = coords[0].position_angle(coords[1])

目录匹配

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

catalog1 = SkyCoord(ra=[10, 11, 12]*u.degree, dec=[40, 41, 42]*u.degree)
catalog2 = SkyCoord(ra=[10.01, 11.02, 13]*u.degree, dec=[40.01, 41.01, 43]*u.degree)

# 查找最近邻居
idx, sep2d, dist3d = catalog1.match_to_catalog_sky(catalog2)

# 按分离阈值过滤
max_sep = 1 * u.arcsec
matched = sep2d < max_sep

水平坐标(高度/方位角)

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

location = EarthLocation(lat=40*u.deg, lon=-70*u.deg, height=300*u.m)
obstime = Time('2023-01-01 03:00:00')
target = SkyCoord(ra=10*u.degree, dec=40*u.degree, frame='icrs')

altaz_frame = AltAz(obstime=obstime, location=location)
target_altaz = target.transform_to(altaz_frame)

print(f"高度: {target_altaz.alt.deg:.2f}°, 方位角: {target_altaz.az.deg:.2f}°")

可用的坐标系

  • icrs - 国际天体参考系(默认,首选)
  • fk5, fk4 - 第五/第四基本星表
  • galactic - 银河系坐标
  • supergalactic - 超级银河系坐标
  • altaz - 水平(高度-方位角)坐标
  • gcrs, cirs, itrs - 地球基础系统
  • 黄道坐标系:BarycentricMeanEcliptic, HeliocentricMeanEcliptic, GeocentricMeanEcliptic

3. 单位和量值

物理单位是天文学计算的基础。Astropy的单位系统提供维度分析和自动转换。

基本单位操作

import astropy.units as u

# 创建量值
distance = 5.2 * u.parsec
velocity = 300 * u.km / u.s
time = 10 * u.year

# 转换单位
distance_ly = distance.to(u.lightyear)
velocity_mps = velocity.to(u.m / u.s)

# 带单位的算术运算
wavelength = 500 * u.nm
frequency = wavelength.to(u.Hz, equivalencies=u.spectral())

处理数组

import numpy as np
import astropy.units as u

wavelengths = np.array([400, 500, 600]) * u.nm
frequencies = wavelengths.to(u.THz, equivalencies=u.spectral())

fluxes = np.array([1.2, 2.3, 1.8]) * u.Jy
luminosities = 4 * np.pi * (10*u.pc)**2 * fluxes

重要的等价关系

  • u.spectral() - 转换波长 ↔ 频率 ↔ 能量
  • u.doppler_optical(rest) - 光学多普勒速度
  • u.doppler_radio(rest) - 无线电多普勒速度
  • u.doppler_relativistic(rest) - 相对论多普勒
  • u.temperature() - 温度单位转换
  • u.brightness_temperature(freq) - 亮度温度

物理常数

from astropy import constants as const

print(const.c)      # 光速
print(const.G)      # 引力常数
print(const.M_sun)  # 太阳质量
print(const.R_sun)  # 太阳半径
print(const.L_sun)  # 太阳光度

性能提示:使用 << 操作符为数组快速分配单位:

# 快速
result = large_array << u.m

# 较慢
result = large_array * u.m

4. 时间系统

天文时间系统需要高精度和多种时间尺度。

创建时间对象

from astropy.time import Time
import astropy.units as u

# 各种输入格式
t1 = Time('2023-01-01T00:00:00', format='isot', scale='utc')
t2 = Time(2459945.5, format='jd', scale='utc')
t3 = Time(['2023-01-01', '2023-06-01'], format='iso')

# 转换格式
print(t1.jd)    # 儒略日
print(t1.mjd)   # 修正儒略日
print(t1.unix)  # Unix时间戳
print(t1.iso)   # ISO格式

# 转换时间尺度
print(t1.tai)   # 国际原子时
print(t1.tt)    # 地球时
print(t1.tdb)   # 太阳系质心力学时

时间算术

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

t1 = Time('2023-01-01T00:00:00')
dt = TimeDelta(1*u.day)

t2 = t1 + dt
diff = t2 - t1
print(diff.to(u.hour))

# 时间数组
times = t1 + np.arange(10) * u.day

天文时间计算

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

location = EarthLocation(lat=40*u.deg, lon=-70*u.deg)
t = Time('2023-01-01T00:00:00')

# 本地恒星时
lst = t.sidereal_time('apparent', longitude=location.lon)

# 太阳系质心修正
target = SkyCoord(ra=10*u.deg, dec=40*u.deg)
ltt = t.light_travel_time(target, location=location)
t_bary = t.tdb + ltt

可用的时间尺度

  • utc - 协调世界时
  • tai - 国际原子时
  • tt - 地球时
  • tcb, tcg - 太阳系质心/地球质心坐标时
  • tdb - 太阳系质心力学时
  • ut1 - 世界时

5. 数据表格

Astropy表格提供天文学特定的增强功能,优于pandas。

创建和操作表格

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

# 创建表格
t = Table()
t['name'] = ['Star1', 'Star2', 'Star3']
t['ra'] = [10.5, 11.2, 12.3] * u.degree
t['dec'] = [41.2, 42.1, 43.5] * u.degree
t['flux'] = [1.2, 2.3, 0.8] * u.Jy

# 列元数据
t['flux'].description = '在1.4 GHz的流量'
t['flux'].format = '.2f'

# 添加计算列
t['flux_mJy'] = t['flux'].to(u.mJy)

# 过滤和排序
bright = t[t['flux'] > 1.0 * u.Jy]
t.sort('flux')

表格输入/输出

from astropy.table import Table

# 读取(格式自动检测扩展名)
t = Table.read('data.fits')
t = Table.read('data.csv', format='ascii.csv')
t = Table.read('data.ecsv', format='ascii.ecsv')  # 保留单位!
t = Table.read('data.votable', format='votable')

# 写入
t.write('output.fits', overwrite=True)
t.write('output.ecsv', format='ascii.ecsv', overwrite=True)

高级操作

from astropy.table import Table, join, vstack, hstack

# 连接表格(类似SQL)
joined = join(table1, table2, keys='id')

# 堆叠表格
combined_rows = vstack([t1, t2])
combined_cols = hstack([t1, t2])

# 分组和聚合
t.group_by('category').groups.aggregate(np.mean)

包含天文对象的表格

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

coords = SkyCoord(ra=[10, 11, 12]*u.deg, dec=[40, 41, 42]*u.deg)
times = Time(['2023-01-01', '2023-01-02', '2023-01-03'])

t = Table([coords, times], names=['coords', 'obstime'])
print(t['coords'][0].ra)  # 访问坐标属性

6. 宇宙学计算

使用标准模型进行快速的宇宙学计算。

使用宇宙学计算器

python scripts/cosmo_calc.py 0.5 1.0 1.5
python scripts/cosmo_calc.py --range 0 3 0.5 --cosmology Planck18
python scripts/cosmo_calc.py 0.5 --verbose
python scripts/cosmo_calc.py --convert 1000 --from luminosity_distance

编程使用

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

cosmo = Planck18

# 计算距离
z = 1.5
d_L = cosmo.luminosity_distance(z)
d_A = cosmo.angular_diameter_distance(z)
d_C = cosmo.comoving_distance(z)

# 时间计算
age = cosmo.age(z)
lookback = cosmo.lookback_time(z)

# 哈勃参数
H_z = cosmo.H(z)

print(f"在z={z}时:")
print(f"  光度距离: {d_L:.2f}")
print(f"  宇宙年龄: {age:.2f}")

转换观测值

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

cosmo = Planck18
z = 1.5

# 角大小到物理大小
d_A = cosmo.angular_diameter_distance(z)
angular_size = 1 * u.arcsec
physical_size = (angular_size.to(u.radian) * d_A).to(u.kpc)

# 流量到光度
flux = 1e-17 * u.erg / u.s / u.cm**2
d_L = cosmo.luminosity_distance(z)
luminosity = flux * 4 * np.pi * d_L**2

# 查找给定距离的红移
from astropy.cosmology import z_at_value
z = z_at_value(cosmo.luminosity_distance, 1000*u.Mpc)

可用的宇宙学模型

  • Planck18, Planck15, Planck13 - 普朗克卫星参数
  • WMAP9, WMAP7, WMAP5 - WMAP卫星参数
  • 自定义:FlatLambdaCDM(H0=70*u.km/u.s/u.Mpc, Om0=0.3)

7. 模型拟合

将数学模型拟合到天文数据。

1D拟合示例

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

# 生成数据
x = np.linspace(0, 10, 100)
y_data = 10 * np.exp(-0.5 * ((x - 5) / 1)**2) + np.random.normal(0, 0.5, x.shape)

# 创建和拟合模型
g_init = models.Gaussian1D(amplitude=8, mean=4.5, stddev=0.8)
fitter = fitting.LevMarLSQFitter()
g_fit = fitter(g_init, x, y_data)

# 结果
print(f"振幅: {g_fit.amplitude.value:.3f}")
print(f"平均值: {g_fit.mean.value:.3f}")
print(f"标准差: {g_fit.stddev.value:.3f}")

# 评估拟合模型
y_fit = g_fit(x)

常见的1D模型

  • Gaussian1D - 高斯轮廓
  • Lorentz1D - 洛伦兹轮廓
  • Voigt1D - 福格特轮廓
  • Moffat1D - 莫法特轮廓(点扩散函数建模)
  • Polynomial1D - 多项式
  • PowerLaw1D - 幂律
  • BlackBody - 黑体光谱

常见的2D模型

  • Gaussian2D - 2D高斯
  • Moffat2D - 2D莫法特(恒星点扩散函数)
  • AiryDisk2D - 艾里盘(衍射模式)
  • Disk2D - 圆盘

带约束的拟合

from astropy.modeling import models, fitting

g = models.Gaussian1D(amplitude=10, mean=5, stddev=1)

# 设置边界
g.amplitude.bounds = (0, None)  # 仅正数
g.mean.bounds = (4, 6)          # 约束中心

# 固定参数
g.stddev.fixed = True

# 复合模型
model = models.Gaussian1D() + models.Polynomial1D(degree=1)

可用的拟合器

  • LinearLSQFitter - 线性最小二乘(快速,适用于线性模型)
  • LevMarLSQFitter - 列文伯格-马夸尔特(最常用)
  • SimplexLSQFitter - 下山单纯形法
  • SLSQPLSQFitter - 带约束的顺序最小二乘

8. 世界坐标系统(WCS)

在图像中转换像素坐标和世界坐标。

基本WCS使用

from astropy.io import fits
from astropy.wcs import WCS

# 读取带WCS的FITS
hdu = fits.open('image.fits')[0]
wcs = WCS(hdu.header)

# 像素到世界
ra, dec = wcs.pixel_to_world_values(100, 200)

# 世界到像素
x, y = wcs.world_to_pixel_values(ra, dec)

# 使用SkyCoord(更强大)
from astropy.coordinates import SkyCoord
import astropy.units as u

coord = SkyCoord(ra=150*u.deg, dec=-30*u.deg)
x, y = wcs.world_to_pixel(coord)

使用WCS绘图

from astropy.io import fits
from astropy.wcs import WCS
from astropy.visualization import ImageNormalize, LogStretch, PercentileInterval
import matplotlib.pyplot as plt

hdu = fits.open('image.fits')[0]
wcs = WCS(hdu.header)
data = hdu.data

# 创建带WCS投影的图形
fig = plt.figure()
ax = fig.add_subplot(111, projection=wcs)

# 使用坐标网格绘图
norm = ImageNormalize(data, interval=PercentileInterval(99.5),
                     stretch=LogStretch())
ax.imshow(data, norm=norm, origin='lower', cmap='viridis')

# 坐标标签和网格
ax.set_xlabel('赤经')
ax.set_ylabel('赤纬')
ax.coords.grid(color='white', alpha=0.5)

9. 统计和数据处处理

稳健的统计工具用于天文数据。

Sigma裁剪(去除异常值):

from astropy.stats import sigma_clip, sigma_clipped_stats

# 去除异常值
clipped = sigma_clip(data, sigma=3, maxiters=5)

# 获取清理后数据的统计信息
mean, median, std = sigma_clipped_stats(data, sigma=3)

# 使用裁剪数据
background = median
signal = data - background
snr = signal / std

其他统计函数

from astropy.stats import mad_std, biweight_location, biweight_scale

# 稳健标准差
std_robust = mad_std(data)

# 稳健中心位置
center = biweight_location(data)

# 稳健尺度
scale = biweight_scale(data)

10. 可视化

以适当的缩放显示天文图像。

图像归一化和拉伸

from astropy.visualization import (ImageNormalize, MinMaxInterval,
                                   PercentileInterval, ZScaleInterval,
                                   SqrtStretch, LogStretch, PowerStretch,
                                   AsinhStretch)
import matplotlib.pyplot as plt

# 常见组合:百分位间隔 + 平方根拉伸
norm = ImageNormalize(data,
                     interval=PercentileInterval(99),
                     stretch=SqrtStretch())

plt.imshow(data, norm=norm, origin='lower', cmap='gray')
plt.colorbar()

可用的间隔(确定最小/最大值):

  • MinMaxInterval() - 使用实际最小/最大值
  • PercentileInterval(percentile) - 裁剪到百分位(例如99%)
  • ZScaleInterval() - IRAF的zscale算法
  • ManualInterval(vmin, vmax) - 手动指定

可用的拉伸(非线性缩放):

  • LinearStretch() - 线性(默认)
  • SqrtStretch() - 平方根(图像常用)
  • LogStretch() - 对数(适用于高动态范围)
  • PowerStretch(power) - 幂律
  • AsinhStretch() - 反双曲正弦(适用于宽范围)

捆绑资源

scripts/

fits_info.py - 全面的FITS文件检查工具

python scripts/fits_info.py observation.fits
python scripts/fits_info.py observation.fits --detailed
python scripts/fits_info.py observation.fits --ext 1

coord_convert.py - 批量坐标转换实用程序

python scripts/coord_convert.py 10.68 41.27 --from icrs --to galactic
python scripts/coord_convert.py --file coords.txt --from icrs --to galactic

cosmo_calc.py - 宇宙学计算器

python scripts/cosmo_calc.py 0.5 1.0 1.5
python scripts/cosmo_calc.py --range 0 3 0.5 --cosmology Planck18

references/

module_overview.md - 所有astropy子包、类和方法的全面参考。查询此文件以获取详细的API信息、可用函数和模块功能。

common_workflows.md - 常见天文数据分析任务的完整工作示例。包含FITS操作、坐标变换、宇宙学、建模和完整分析管道的完整代码示例。

最佳实践

  1. 对FITS文件使用上下文管理器

    with fits.open('file.fits') as hdul:
        # 处理文件
    
  2. 对于更好的单位/元数据支持,优先使用 astropy.table 而不是原始FITS表格

  3. 使用 SkyCoord 处理坐标(高级接口),而不是低级框架类

  4. 尽可能为量值附加单位,以确保维度安全

  5. 保存表格时使用 ECSV 格式,以保留单位和元数据

  6. 为性能起见,对坐标操作进行向量化,而不是循环

  7. 打开大型FITS文件时使用 memmap=True,以节省内存

  8. 安装 Bottleneck 包 以加速统计操作

  9. 对重复操作预计算复合单位,以提高性能

  10. 查询 references/module_overview.md 以获取详细的模块信息,以及 references/common_workflows.md** 以获取完整的工作示例

常见模式

模式:FITS → 处理 → FITS

from astropy.io import fits
from astropy.stats import sigma_clipped_stats

# 读取
with fits.open('input.fits') as hdul:
    data = hdul[0].data
    header = hdul[0].header

    # 处理
    mean, median, std = sigma_clipped_stats(data, sigma=3)
    processed = (data - median) / std

    # 写入
    fits.writeto('output.fits', processed, header, overwrite=True)

模式:目录匹配

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

# 加载目录
cat1 = Table.read('catalog1.fits')
cat2 = Table.read('catalog2.fits')

# 创建坐标对象
coords1 = SkyCoord(ra=cat1['RA'], dec=cat1['DEC'], unit=u.degree)
coords2 = SkyCoord(ra=cat2['RA'], dec=cat2['DEC'], unit=u.degree)

# 匹配
idx, sep2d, dist3d = coords1.match_to_catalog_sky(coords2)

# 按分离过滤
max_sep = 1 * u.arcsec
matched_mask = sep2d < max_sep

# 创建匹配目录
matched_cat1 = cat1[matched_mask]
matched_cat2 = cat2[idx[matched_mask]]

模式:时间序列分析

from astropy.time import Time
from astropy.timeseries import TimeSeries
import astropy.units as u

# 创建时间序列
times = Time(['2023-01-01', '2023-01-02', '2023-01-03'])
flux = [1.2, 2.3, 1.8] * u.Jy

ts = TimeSeries(time=times)
ts['flux'] = flux

# 按周期折叠
from astropy.timeseries import aggregate_downsample
period = 1.5 * u.day
folded = ts.fold(period=period)

模式:带WCS的图像显示

from astropy.io import fits
from astropy.wcs import WCS
from astropy.visualization import ImageNormalize, SqrtStretch, PercentileInterval
import matplotlib.pyplot as plt

hdu = fits.open('image.fits')[0]
wcs = WCS(hdu.header)
data = hdu.data

fig = plt.figure(figsize=(10, 10))
ax = fig.add_subplot(111, projection=wcs)

norm = ImageNormalize(data, interval=PercentileInterval(99),
                     stretch=SqrtStretch())
im = ax.imshow(data, norm=norm, origin='lower', cmap='viridis')

ax.set_xlabel('赤经')
ax.set_ylabel('赤纬')
ax.coords.grid(color='white', alpha=0.5, linestyle='solid')
plt.colorbar(im, ax=ax)

安装说明

确保Python环境中安装了astropy:

pip install astropy

如需额外的性能和功能:

pip install astropy[all]  # 包括可选依赖项

额外资源