astropy-astronomy

Astropy — 天文学・天体物理学ツールキット

概要

Astropy は、天文学向けの中核的な Python パッケージであり、天文学研究に不可欠な機能を提供します。単位を考慮した計算、天球座標変換、FITS ファイルの入出力、宇宙論計算、高精度な時刻処理、表形式データ操作、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.184 秒）	地心計算
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') を使用します
座標変換が遅い	座標を 1 つずつループ処理している	配列形式の 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 の「主要な概念」（概要表）と 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/