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操作、坐标变换、宇宙学、建模和完整分析管道的完整代码示例。
最佳实践
-
对FITS文件使用上下文管理器:
with fits.open('file.fits') as hdul: # 处理文件 -
对于更好的单位/元数据支持,优先使用 astropy.table 而不是原始FITS表格
-
使用 SkyCoord 处理坐标(高级接口),而不是低级框架类
-
尽可能为量值附加单位,以确保维度安全
-
保存表格时使用 ECSV 格式,以保留单位和元数据
-
为性能起见,对坐标操作进行向量化,而不是循环
-
打开大型FITS文件时使用 memmap=True,以节省内存
-
安装 Bottleneck 包 以加速统计操作
-
对重复操作预计算复合单位,以提高性能
-
查询
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] # 包括可选依赖项
额外资源
- 官方文档:https://docs.astropy.org
- 教程:https://learn.astropy.org
- API参考:查询本技能中的
references/module_overview.md - 工作示例:查询本技能中的
references/common_workflows.md