---

name: astropy-astronomy description: "Core Python library for astronomy/astrophysics: units with dimensional analysis, celestial coordinate transforms (ICRS/Galactic/AltAz/FK5), FITS I/O, tables (FITS/HDF5/VOTable/CSV), cosmology (Planck18, distance/age), precise time (UTC/TAI/TT/TDB, Julian, barycentric), WCS pixel-world mapping, model fitting. For general tables use pandas/polars; for radio interferometry use CASA." license: BSD-3-Clause

Astropy — Astronomy & Astrophysics Toolkit

Overview

Astropy is the core Python package for astronomy, providing essential functionality for astronomical research: unit-aware calculations, celestial coordinate transformations, FITS file I/O, cosmological calculations, precise time handling, tabular data operations, and WCS image coordinate mapping.

When to Use

• Converting between celestial coordinate systems (ICRS, Galactic, FK5, AltAz)
• Working with physical quantities and units (Jy→mJy, parsec→km, spectral equivalencies)
• Reading, writing, or manipulating FITS files (images and tables)
• Cosmological calculations (luminosity distance, lookback time, comoving volume)
• Precise time handling with multiple scales (UTC, TAI, TT, TDB) and formats (JD, MJD, ISO)
• Cross-matching astronomical catalogs by sky position
• WCS transformations between pixel and world coordinates
• For general tabular data: use pandas or polars instead
• For radio interferometry: use CASA instead

Prerequisites

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

Quick Start

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}")

Core API

1. Units & Quantities (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. Coordinate Systems (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 File Handling (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. Table Operations (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. Time Handling (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. Cosmological Calculations (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 & Image Processing

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}")

Key Concepts

Unit Equivalency System

Astropy's equivalencies parameter enables domain-specific conversions that are not dimensionally equivalent:

| Equivalency | Converts Between | Example | |-------------|-----------------|---------| | u.spectral() | Wavelength ↔ frequency