A comprehensive skill for astronomical data analysis using the Python astronomy ecosystem.
Covers reading and writing FITS files, World Coordinate System (WCS) transforms, remote
catalog queries, aperture photometry on imaging data, and CMB angular power spectrum
analysis using HEALPix pixelisation.
When to Use This Skill
Use this skill when you need to:
Read, write, or inspect FITS images and tables produced by telescopes or pipelines
Convert between pixel coordinates and sky coordinates (RA/Dec, Galactic, Ecliptic)
Cross-match source lists against Vizier or SIMBAD to identify objects
Perform aperture or PSF photometry on CCD images
Analyse CMB or large-scale-structure maps stored as HEALPix .fits files
Automate catalogue retrieval (Gaia, 2MASS, SDSS, NED) without manual downloads
Build reproducible astronomy pipelines that run on any machine
This skill is not appropriate for:
Real-time telescope control or instrument communication (use INDI / ASCOM instead)
Spectral extraction from 2-D spectra (use specutils skill)
N-body or hydrodynamical simulations (use yt or pynbody skill)
Background & Key Concepts
FITS — Flexible Image Transport System
FITS is the standard file format in astronomy. A FITS file contains one or more
Header/Data Units (HDUs). Each HDU has a text header of 80-character keyword cards
and an optional data array or binary table. Common HDU types:
HDU type
astropy class
Typical use
PrimaryHDU
fits.PrimaryHDU
2-D image, metadata-only
ImageHDU
fits.ImageHDU
Additional image planes
BinTableHDU
fits.BinTableHDU
Source catalogs, spectra
CompImageHDU
fits.CompImageHDU
Tile-compressed image
WCS — World Coordinate System
The FITS WCS standard maps pixel indices (x, y) to sky coordinates
(RA, Dec) via a projection (TAN, ZEA, CAR, …) plus a linear transform
encoded in the header keywords CRPIX, CRVAL, / .
reads these keywords and provides /
methods that return objects.
CD
CDELT
astropy.wcs.WCS
pixel_to_world
world_to_pixel
SkyCoord
SkyCoord and Angle
astropy.coordinates.SkyCoord represents one or more positions on the sky in
any supported reference frame. Angle arithmetic, frame conversions, and
on-sky separations are all handled automatically. Key frames: ICRS (RA/Dec),
Galactic (l, b), FK5, AltAz (horizontal coordinates).
Catalog Cross-Matching
astroquery provides uniform Python interfaces to dozens of online services:
Vizier — CDS catalogue service (Gaia, 2MASS, SDSS, …)
SIMBAD — Object identification database
NED — NASA/IPAC Extragalactic Database
ESASky / MAST / IRSA — Archive portals
All queries return astropy.table.Table objects directly.
HEALPix and healpy
HEALPix (Hierarchical Equal Area isoLatitude Pixelisation) divides the sphere
into 12 * nside**2 equal-area pixels. healpy wraps the C++ HEALPix
library and provides map I/O (read_map, write_map), spherical harmonic
transforms (map2alm, alm2cl), and visualisation (mollview, gnomview).
The angular power spectrum C_l is the key observable for CMB cosmology.
Aperture Photometry
photutils (companion to astropy) implements circular and elliptical apertures,
sky annuli for local background estimation, and source detection (DAOStarFinder,
IRAFStarFinder). The result is a flux in counts per second that can be converted
to magnitudes using a photometric zero-point.
# Fetch a small SDSS r-band cutout for testing
python - <<'EOF'
from astroquery.skyview import SkyView
imgs = SkyView.get_images("M51", survey=["SDSSr"], pixels=512)
imgs[0].writeto("m51_r.fits", overwrite=True)
print("Saved m51_r.fits")
EOF
Environment variables
No API keys are required for the public Vizier / SIMBAD / SkyView services.
If you use the ESO archive or proprietary data portals, store credentials as:
export ESO_USERNAME="<your-username>"export ESO_PASSWORD=$(cat ~/.eso_passwd) # read from file, never hardcode
Access them in Python with:
import os
username = os.getenv("ESO_USERNAME", "")
password = os.getenv("ESO_PASSWORD", "")
Core Workflow
Step 1 — Read and inspect a FITS file
from astropy.io import fits
import numpy as np
# Open without loading data into memorywith fits.open("m51_r.fits") as hdul:
hdul.info() # print HDU summary
header = hdul[0].header
data = hdul[0].data.astype(float) # pixel array (counts)print(f"Image shape : {data.shape}")
print(f"Instrument : {header.get('INSTRUME', 'unknown')}")
print(f"Filter : {header.get('FILTER', 'unknown')}")
print(f"Exposure : {header.get('EXPTIME', 'unknown')} s")
print(f"Min / Max : {data.min():.1f} / {data.max():.1f} counts")
# Mask NaN / Inf pixels that can corrupt photometry
data = np.where(np.isfinite(data), data, 0.0)
Step 2 — WCS coordinate transforms
from astropy.wcs import WCS
from astropy.coordinates import SkyCoord
import astropy.units as u
with fits.open("m51_r.fits") as hdul:
wcs = WCS(hdul[0].header)
data = hdul[0].data.astype(float)
ny, nx = data.shape
# Pixel centre of the image -> sky coordinates
cx, cy = nx / 2.0, ny / 2.0
sky_centre = wcs.pixel_to_world(cx, cy)
print(f"Image centre: RA={sky_centre.ra.deg:.4f} deg, "f"Dec={sky_centre.dec.deg:.4f} deg")
# Convert to Galactic coordinates
gal = sky_centre.galactic
print(f"Galactic : l={gal.l.deg:.4f} deg, b={gal.b.deg:.4f} deg")
# Sky -> pixel (useful for placing apertures on known objects)
m51_coord = SkyCoord(ra=202.4696 * u.deg, dec=47.1952 * u.deg)
pix = wcs.world_to_pixel(m51_coord)
print(f"M51 nucleus at pixel ({pix[0]:.1f}, {pix[1]:.1f})")
# Compute pixel scale
pixel_scale = wcs.proj_plane_pixel_scales()
print(f"Pixel scale : {pixel_scale[0].to(u.arcsec):.3f} / pixel")
Step 3 — Catalog query with astroquery (Vizier + SIMBAD)
from astroquery.vizier import Vizier
from astroquery.simbad import Simbad
from astropy.coordinates import SkyCoord
import astropy.units as u
# -- Gaia DR3 point sources within 5 arcmin of M51 --
coord = SkyCoord(ra=202.4696 * u.deg, dec=47.1952 * u.deg, frame="icrs")
radius = 5 * u.arcmin
v = Vizier(columns=["Source", "RA_ICRS", "DE_ICRS", "Gmag", "BP-RP"],
row_limit=500)
result = v.query_region(coord, radius=radius, catalog="I/355/gaiadr3")
if result:
gaia = result[0]
print(f"Found {len(gaia)} Gaia DR3 sources")
print(gaia["Source", "RA_ICRS", "DE_ICRS", "Gmag"][:5])
else:
print("No Gaia sources returned (check network)")
# -- SIMBAD identification of M51 --
simbad = Simbad()
simbad.add_votable_fields("distance", "flux(V)", "sptype")
result_s = simbad.query_object("M51")
if result_s:
print("\nSIMBAD result for M51:")
print(result_s["MAIN_ID", "RA", "DEC", "FLUX_V"])
