Astropy Guide: FITS Images and Coordinate Systems

Astropy is the foundational Python library for professional astronomy, providing robust subpackages to handle specialized data formats and complex physical calculations. This guide demonstrates how Astropy manipulates Flexible Image Transport System (FITS) files using astropy.io.fits and astropy.wcs, and how it performs precise celestial coordinate transformations using astropy.coordinates and astropy.units.

Manipulating FITS Images with Astropy

The FITS format organizes astronomical data into Header Data Units (HDUs), each containing metadata headers and data arrays. Astropy interacts with these files through the astropy.io.fits module.

Opening and Inspecting FITS Files

To open and inspect the contents of a FITS file:

from astropy.io import fits

# Open the FITS file
with fits.open('sample_image.fits') as hdul:
    # Print a summary of the HDU list
    hdul.info()
    
    # Access the primary HDU
    primary_hdu = hdul[0]
    
    # Read the header and image data
    header = primary_hdu.header
    image_data = primary_hdu.data

image_data is returned as a standard NumPy array, allowing direct mathematical operations, slicing, and statistical analysis (such as calculating background noise, mean flux, or cropping specific regions).

Modifying and Writing FITS Data

Headers function like Python dictionaries, mapping standard astronomical keywords to values and comments. You can update headers and modify pixel values before exporting:

# Modify a header keyword
header['OBSERVER'] = 'Jane Doe'

# Perform an operation on image pixels (e.g., flat-fielding or thresholding)
modified_data = image_data - 100.0

# Write the modified data to a new FITS file
new_hdu = fits.PrimaryHDU(data=modified_data, header=header)
new_hdu.writeto('processed_image.fits', overwrite=True)

Mapping Pixels with World Coordinate System (WCS)

Astropy's astropy.wcs module translates raw pixel coordinates into standard sky coordinates using the projection parameters stored in the FITS header:

from astropy.wcs import WCS

# Initialize WCS from the header
wcs = WCS(header)

# Convert pixel coordinates (X=100, Y=200) to Right Ascension and Declination
ra, dec = wcs.pixel_to_world_values(100, 200)

Calculating Celestial Coordinate Conversions

The astropy.coordinates subpackage, integrated with astropy.units, provides a framework for representing and transforming spherical coordinates between different astronomical reference frames.

Representing Coordinates with SkyCoord

The core object is SkyCoord, which accepts coordinates in various formats, including decimal degrees, radians, or sexagesimal strings:

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

# Define a coordinate using decimal degrees
coord_icrs = SkyCoord(ra=280.0 * u.deg, dec=-15.0 * u.deg, frame='icrs')

# Define a coordinate using sexagesimal format
coord_hms = SkyCoord('18h40m00s', '-15d00m00s', frame='icrs')

Frame Transformations

Astropy supports transformations across standard reference systems, including ICRS (International Celestial Reference System), Galactic, FK5, and horizontal (AltAz) systems.

To convert equatorial coordinates to Galactic coordinates:

# Transform ICRS to the Galactic coordinate system
galactic_coord = coord_icrs.galactic

print(f"Galactic Longitude (l): {galactic_coord.l.deg:.4f} deg")
print(f"Galactic Latitude (b): {galactic_coord.b.deg:.4f} deg")

Topocentric (AltAz) Conversions

Converting to the local observer frame (Altitude-Azimuth) requires specifying an observation time and Earth location:

from astropy.coordinates import EarthLocation, AltAz
from astropy.time import Time

# Define observer location and time
location = EarthLocation(lat=31.9583 * u.deg, lon=-111.5967 * u.deg, height=2090 * u.m)
time = Time('2026-03-31 04:00:00')

# Create the AltAz reference frame
altaz_frame = AltAz(obstime=time, location=location)

# Transform the sky coordinate to AltAz
local_coord = coord_icrs.transform_to(altaz_frame)

print(f"Altitude: {local_coord.alt.deg:.2f} deg")
print(f"Azimuth: {local_coord.az.deg:.2f} deg")

Measuring Angular Separation

Calculating the angular distance between two sky positions handles spherical geometry automatically:

target1 = SkyCoord(ra=10.684 * u.deg, dec=41.269 * u.deg, frame='icrs')
target2 = SkyCoord(ra=10.700 * u.deg, dec=41.300 * u.deg, frame='icrs')

separation = target1.separation(target2)
print(f"Separation: {separation.arcmin:.3f} arcminutes")