astropy

Astropy - 天文学・天体物理学

Astropyは単なるライブラリではなく、コミュニティの標準です。数値配列を物理単位や天球座標系と結び付けることで、天文計算の再現性と物理的な整合性を確保します。

使用する場面

• 計算で物理単位や定数を扱う場合。
• 天球座標（赤経、赤緯）や座標系の変換を扱う場合。
• FITS（Flexible Image Transport System）、ASDF、VO（Virtual Observatory）ファイルを読み書きする場合。
• 天文時系（UTC、TAI、TDB、ユリウス日）を管理する場合。
• World Coordinate System（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、つまり単位と組み合わせた数値または配列として扱います。これにより、メートルとフィートを加算するようなエラーを防げます。

天球座標

座標は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を優先する - 列をQuantityオブジェクトとして保持するため、Tableの代わりにQTableを使ってください。
• ロジックにはSkyCoordを使う - 角距離を自分で計算せず、c1.separation(c2)を使ってください。
• 座標をベクトル化する - SkyCoordには位置の配列を格納できます。ループよりもはるかに高速です。
• 時系を指定する - Timeオブジェクトの作成時には、必ず時系（utc、tdb）を指定してください。
• FITSにはwithを使う - ファイルハンドルを確実に閉じるため、必ずコンテキストマネージャーを使ってください。

❌ しないこと

• 早い段階で単位を取り除かない - 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は天文学における再現可能な研究の基盤です。単位を通じて物理的な整合性を徹底し、座標の扱いを標準化することで、かつて天文ソフトウェアを悩ませていた、気付かれにくいエラーを防ぎます。