Astropy for Astronomy — Deep Dive

Master Astropy's coordinate transformations, WCS projections, unit algebra, FITS manipulation, and integration with observatory pipelines.

Architecture and design philosophy

Astropy follows a layered architecture. The core package provides foundational infrastructure that affiliated packages extend for specialized tasks. This separation keeps the core stable and well-tested while allowing rapid innovation at the edges.

Key design principles:

Units are first-class citizens. Quantities propagate through all operations, eliminating silent unit errors.
Interoperability over reinvention. Astropy wraps proven C libraries (ERFA for astrometry, cfitsio-inspired I/O) rather than reimplementing them in Python.
Community governance. An Astropy Project Team coordinates development across packages, ensuring compatibility and preventing ecosystem fragmentation.

Units and quantities in depth

import astropy.units as u

# Basic quantity creation
distance = 4.22 * u.lightyear  # Proxima Centauri
speed = 300_000 * u.km / u.s

# Unit conversion
distance_pc = distance.to(u.pc)
print(f"Proxima Centauri: {distance_pc:.3f}")
# 1.293 pc

# Arithmetic with units
travel_time = distance / speed
print(f"Travel time at light speed: {travel_time.to(u.year):.2f}")

# Decompose to SI
energy = 13.6 * u.eV  # hydrogen ionization energy
print(f"In joules: {energy.to(u.J):.3e}")

Unit equivalencies

Some conversions are not simple scaling — they require physics. Astropy handles this with equivalencies:

import astropy.units as u

# Wavelength ↔ frequency ↔ energy (for photons)
wavelength = 656.3 * u.nm  # hydrogen-alpha line

frequency = wavelength.to(u.GHz, equivalencies=u.spectral())
energy = wavelength.to(u.eV, equivalencies=u.spectral())

print(f"H-alpha: {wavelength} = {frequency:.1f} = {energy:.3f}")

# Temperature ↔ energy (thermal)
temp = 5778 * u.K  # solar surface temperature
thermal_energy = temp.to(u.eV, equivalencies=u.temperature_energy())
print(f"Solar surface kT: {thermal_energy:.4f}")

Custom units

import astropy.units as u

# Define a custom unit
solar_luminosity = u.def_unit("Lsun", 3.828e26 * u.W)

star_luminosity = 100 * solar_luminosity
print(f"Star luminosity: {star_luminosity.to(u.W):.3e}")

Coordinate transformations

SkyCoord basics

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

# Define a position (Betelgeuse)
betelgeuse = SkyCoord(ra=88.7929 * u.deg, dec=7.4070 * u.deg, frame="icrs")

# Convert to Galactic coordinates
galactic = betelgeuse.galactic
print(f"Galactic: l={galactic.l:.2f}, b={galactic.b:.2f}")

# Convert to AltAz (needs observer location and time)
from astropy.coordinates import EarthLocation, AltAz
from astropy.time import Time

location = EarthLocation(lat=31.9583 * u.deg, lon=-111.5967 * u.deg, height=2097 * u.m)  # Kitt Peak
time = Time("2026-01-15 04:00:00", scale="utc")

altaz_frame = AltAz(obstime=time, location=location)
altaz = betelgeuse.transform_to(altaz_frame)
print(f"Alt: {altaz.alt:.1f}, Az: {altaz.az:.1f}")

Separation and matching

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

# Angular separation between two objects
m31 = SkyCoord(ra=10.6847 * u.deg, dec=41.2687 * u.deg)
m33 = SkyCoord(ra=23.4621 * u.deg, dec=30.6602 * u.deg)

sep = m31.separation(m33)
print(f"M31 to M33: {sep.deg:.2f} degrees")

# Cross-match two catalogs
import numpy as np

catalog1 = SkyCoord(ra=[10.0, 20.0, 30.0] * u.deg, dec=[5.0, 15.0, 25.0] * u.deg)
catalog2 = SkyCoord(ra=[10.01, 30.02, 50.0] * u.deg, dec=[5.01, 25.01, 45.0] * u.deg)

idx, sep2d, _ = catalog1.match_to_catalog_sky(catalog2)
print(f"Matches: {idx}")
print(f"Separations: {sep2d.arcsec}")

Time handling

from astropy.time import Time

# Create from different formats
t1 = Time("2026-03-28T00:00:00", format="isot", scale="utc")
t2 = Time(2460760.5, format="jd")  # Julian Date

# Convert between scales
print(f"UTC: {t1.utc.iso}")
print(f"TAI: {t1.tai.iso}")
print(f"TDB: {t1.tdb.iso}")
print(f"MJD: {t1.mjd}")

# Time arithmetic
from astropy.time import TimeDelta

dt = TimeDelta(365.25, format="jd")
t_next_year = t1 + dt
print(f"One year later: {t_next_year.iso}")

# Sidereal time (local apparent)
from astropy.coordinates import EarthLocation
import astropy.units as u

location = EarthLocation(lon=-118.0 * u.deg)
lst = t1.sidereal_time("apparent", longitude=location.lon)
print(f"Local sidereal time: {lst}")

FITS file operations

Reading images

from astropy.io import fits
import numpy as np

# Open a FITS file
hdul = fits.open("observation.fits")
hdul.info()  # Lists all HDUs (Header Data Units)

# Access the primary image
header = hdul[0].header
data = hdul[0].data  # NumPy array

print(f"Image shape: {data.shape}")
print(f"Telescope: {header.get('TELESCOP', 'unknown')}")
print(f"Exposure: {header.get('EXPTIME', 0)} seconds")
print(f"Filter: {header.get('FILTER', 'unknown')}")

hdul.close()

World Coordinate System (WCS)

WCS maps pixel coordinates to sky coordinates. It is encoded in the FITS header.

from astropy.wcs import WCS
from astropy.io import fits
from astropy.coordinates import SkyCoord
import astropy.units as u

header = fits.getheader("image.fits")
wcs = WCS(header)

# Pixel → sky
sky = wcs.pixel_to_world(512, 512)
print(f"Center pixel maps to: RA={sky.ra.deg:.4f}, Dec={sky.dec.deg:.4f}")

# Sky → pixel
target = SkyCoord(ra=150.0 * u.deg, dec=2.0 * u.deg)
px, py = wcs.world_to_pixel(target)
print(f"Target at pixel: ({px:.1f}, {py:.1f})")

Creating FITS files

from astropy.io import fits
import numpy as np

# Create a synthetic image
image_data = np.random.poisson(lam=100, size=(1024, 1024)).astype(np.float32)

# Build header
header = fits.Header()
header["TELESCOP"] = "SimTelescope"
header["EXPTIME"] = 300.0
header["FILTER"] = "R"
header["DATE-OBS"] = "2026-03-28T03:15:00"

# Write
hdu = fits.PrimaryHDU(data=image_data, header=header)
hdu.writeto("synthetic.fits", overwrite=True)

Tables and catalogs

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

# Create a catalog
catalog = Table()
catalog["name"] = ["Sirius", "Betelgeuse", "Rigel"]
catalog["ra"] = [101.287, 88.793, 78.634] * u.deg
catalog["dec"] = [-16.716, 7.407, -8.202] * u.deg
catalog["vmag"] = [-1.46, 0.42, 0.13] * u.mag
catalog["distance"] = [2.64, 168.1, 264.0] * u.pc

# Filter bright stars
bright = catalog[catalog["vmag"] < 0.5 * u.mag]
print(bright)

# Write to FITS table
catalog.write("star_catalog.fits", format="fits", overwrite=True)

# Read from VOTable (Virtual Observatory standard)
# vot = Table.read("catalog.xml", format="votable")

Cosmological calculations

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

# Lookback time and distances for a galaxy at redshift z=1
z = 1.0

lookback = cosmo.lookback_time(z)
luminosity_dist = cosmo.luminosity_distance(z)
comoving_dist = cosmo.comoving_distance(z)
angular_size = cosmo.angular_diameter_distance(z)

print(f"Redshift z={z}")
print(f"  Lookback time: {lookback:.2f}")
print(f"  Luminosity distance: {luminosity_dist:.1f}")
print(f"  Comoving distance: {comoving_dist:.1f}")
print(f"  Angular diameter distance: {angular_size:.1f}")
print(f"  Age of universe at z={z}: {cosmo.age(z):.2f}")

Integration with the ecosystem

Querying remote databases

# Using astroquery (affiliated package)
from astroquery.simbad import Simbad

result = Simbad.query_object("M31")
print(result["RA", "DEC", "FLUX_V"])

Photometry workflow sketch

from astropy.io import fits
from astropy.wcs import WCS
from astropy.coordinates import SkyCoord
from astropy.table import Table
import astropy.units as u

# 1. Load image
data = fits.getdata("science_frame.fits")
wcs = WCS(fits.getheader("science_frame.fits"))

# 2. Source detection (using photutils)
# from photutils.detection import DAOStarFinder
# sources = DAOStarFinder(fwhm=3.0, threshold=5*std)(data - median)

# 3. Convert pixel positions to sky coordinates
# sky_positions = wcs.pixel_to_world(sources['xcentroid'], sources['ycentroid'])

# 4. Aperture photometry
# from photutils.aperture import CircularAperture, aperture_photometry
# apertures = CircularAperture(zip(sources['xcentroid'], sources['ycentroid']), r=5.0)
# phot_table = aperture_photometry(data, apertures)

# 5. Save catalog
# phot_table.write("photometry_catalog.fits", format="fits")

Performance considerations

Operation	Typical scale	Tips
FITS image I/O	100 MB – 10 GB	Use `memmap=True` (default) for large files
Coordinate transforms	Millions of positions	Vectorized — pass arrays, not loops
Cross-matching	10M × 10M catalogs	Use `search_around_sky` for cone searches
Time conversions	Microsecond precision	Specify `precision` parameter for sub-second work
Table operations	Millions of rows	Use `Table.read` with `memmap=True` for FITS tables

Lazy loading

# Memory-mapped FITS — data stays on disk until accessed
hdul = fits.open("huge_image.fits", memmap=True)
# Only the accessed slice loads into RAM
cutout = hdul[0].data[100:200, 100:200]

Common pitfalls

Forgetting the unit on bare numbers. SkyCoord(ra=10.5, dec=41.2) fails because Astropy does not know if those are degrees or radians. Always attach u.deg or u.rad.
Mixing coordinate frames. Adding ICRS and Galactic coordinates without explicit transformation gives wrong results silently if you bypass Astropy’s type system.
Ignoring leap seconds in time comparisons. UTC has leap seconds; TAI does not. Subtracting UTC times across a leap second boundary requires Astropy’s awareness of the leap second table.
Assuming FITS images are oriented north-up. Many instruments rotate the image. Always use the WCS to determine sky orientation.

The one thing to remember: Astropy is the bedrock infrastructure for Python astronomy — its unit-aware quantities, coordinate transforms, time systems, and FITS tools prevent entire categories of bugs and let researchers focus on the science rather than the plumbing.

pythonastronomyscience