astropy-astronomyAstronomical data analysis with astropy and astroquery — FITS I/O, WCS transforms, catalog cross-matching, aperture photometry, and CMB power spectra.
Independent GitHub snapshot · Automated classification · Original documentation governs
name: astropy-astronomy description: Astronomical data analysis with astropy and astroquery — FITS I/O, WCS transforms, catalog cross-matching, aperture photometry, and CMB power spectra. tags:
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.
Use this skill when you need to:
.fits filesThis skill is not appropriate for:
specutils skill)yt or pynbody skill)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 |
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, CD / CDELT.
astropy.wcs.WCS reads these keywords and provides pixel_to_world /
world_to_pixel methods that return SkyCoord objects.
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).
astroquery provides uniform Python interfaces to dozens of online services:
All queries return astropy.table.Table objects directly.
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.
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.
pip install "astropy>=5.3" "astroquery>=0.4.6" "healpy>=1.16" \
"matplotlib>=3.7" "numpy>=1.24" photutils
On Apple Silicon / conda:
conda install -c conda-forge astropy astroquery healpy matplotlib numpy photutils
# 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
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", "")
from astropy.io import fits
import numpy as np
# Open without loading data into memory
with 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)
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")
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"])
import numpy as np
from astropy.io import fits
from astropy.wcs import WCS
from astropy.coordinates import SkyCoord
import astropy.units as u
from photutils.detection import DAOStarFinder
from photutils.aperture import (CircularAperture, CircularAnnulus,
aperture_photometry)
from photutils.background import MMMBackground
with fits.open("m51_r.fits") as hdul:
data = hdul[0].data.astype(float)
header = hdul[0].header
wcs = WCS(header)
# Estimate and subtract sky background
bkg_estimator = MMMBackground()
bkg = bkg_estimator(data)
data_sub = data - bkg
print(f"Background level: {bkg:.2f} counts")
# Detect sources with DAOStarFinder
from astropy.stats import sigma_clipped_stats
_, median, std = sigma_clipped_stats(data, sigma=3.0)
finder = DAOStarFinder(fwhm=3.5, threshold=5.0 * std)
sources = finder(data_sub)
print(f"Detected {len(sources)} sources")
# Perform circular aperture photometry (r=5 px, sky annulus 7-10 px)
positions = np.column_stack([sources["xcentroid"], sources["ycentroid"]])
aperture = CircularAperture(positions, r=5.0)
annulus = CircularAnnulus(positions, r_in=7.0, r_out=10.0)
phot_table = aperture_photometry(data_sub, [aperture, annulus])
# Local background correction
annulus_area = annulus.area
aperture_area = aperture.area
sky_per_pixel = phot_table["aperture_sum_1"] / annulus_area
sky_in_aper = sky_per_pixel * aperture_area
flux_corrected = phot_table["aperture_sum_0"] - sky_in_aper
# Instrumental magnitude (zero-point = 25.0 as example)
ZERO_POINT = 25.0
exptime = header.get("EXPTIME", 1.0)
flux_rate = flux_corrected / exptime
mag = ZERO_POINT - 2.5 * np.log10(np.abs(flux_rate))
phot_table["flux_corrected"] = flux_corrected
phot_table["mag_inst"] = mag
phot_table.sort("mag_inst")
print(phot_table["id", "xcenter", "ycenter", "flux_corrected", "mag_inst"][:10])
import numpy as np
import healpy as hp
import matplotlib.pyplot as plt
# ---- Simulate a Gaussian CMB temperature map ----
nside = 512 # HEALPix resolution (npix = 12*nside^2)
lmax = 3 * nside - 1 # maximum multipole
npix = hp.nside2npix(nside)
# Planck 2018-like flat LCDM power spectrum (toy model)
ells = np.arange(lmax + 1)
Cl = np.zeros(lmax + 1)
# Simple Sachs-Wolfe plateau + first acoustic peak shape
Cl[2:] = 1e-10 * (1.0 / (ells[2:] * (ells[2:] + 1))) * (
1.0 + 6.0 * np.exp(-((ells[2:] - 200) / 40) ** 2)
)
# Generate a realisation
np.random.seed(42)
cmb_map = hp.synfast(Cl, nside=nside, lmax=lmax, verbose=False)
print(f"Map shape : {cmb_map.shape}, RMS = {cmb_map.std():.3e} K")
# Recover power spectrum from the map
alm = hp.map2alm(cmb_map, lmax=lmax)
Cl_rec = hp.alm2cl(alm)
ell_rec = np.arange(len(Cl_rec))
# Visualise
fig, axes = plt.subplots(1, 2, figsize=(14, 5))
hp.mollview(cmb_map, title="Simulated CMB (K)", unit="K",
hold=False, fig=fig.number, sub=(1, 2, 1))
ax = axes[1]
ax.plot(ell_rec[2:], ell_rec[2:] * (ell_rec[2:] + 1) * Cl_rec[2:] / (2 * np.pi),
lw=1.5, label=r"Recovered $C_\ellquot;)
ax.set_xlabel(r"Multipole $\ellquot;)
ax.set_ylabel(r"$\ell(\ell+1)C_\ell / 2\pi$ [K$^2$]")
ax.set_title("Angular Power Spectrum")
ax.legend()
fig.tight_layout()
plt.savefig("cmb_power_spectrum.png", dpi=150)
print("Saved cmb_power_spectrum.png")
from astropy.coordinates import SkyCoord, match_coordinates_sky
import astropy.units as u
import numpy as np
# Suppose cat_a and cat_b are astropy Tables with RA/Dec columns
# Here we mock them with random positions for illustration
rng = np.random.default_rng(0)
n_a, n_b = 200, 5000
ra_a = rng.uniform(202.0, 203.0, n_a)
dec_a = rng.uniform(46.8, 47.5, n_a)
ra_b = rng.uniform(201.5, 203.5, n_b)
dec_b = rng.uniform(46.5, 47.8, n_b)
coords_a = SkyCoord(ra=ra_a * u.deg, dec=dec_a * u.deg, frame="icrs")
coords_b = SkyCoord(ra=ra_b * u.deg, dec=dec_b * u.deg, frame="icrs")
# For each source in A, find the nearest in B
idx, sep2d, _ = match_coordinates_sky(coords_a, coords_b)
# Keep only matches within 1 arcsec
match_mask = sep2d < 1.0 * u.arcsec
n_match = match_mask.sum()
print(f"Matched {n_match} / {n_a} sources within 1 arcsec")
matched_a_idx = np.where(match_mask)[0]
matched_b_idx = idx[match_mask]
separations = sep2d[match_mask].to(u.arcsec)
print(f"Median separation : {np.median(separations):.3f}")
import numpy as np
from astropy.io import fits
from astropy.wcs import WCS
ny, nx = 512, 512
wcs = WCS(naxis=2)
wcs.wcs.ctype = ["RA---TAN", "DEC--TAN"]
wcs.wcs.crpix = [nx / 2 + 0.5, ny / 2 + 0.5] # reference pixel (1-based)
wcs.wcs.crval = [202.4696, 47.1952] # RA, Dec of reference (deg)
wcs.wcs.cdelt = [-0.000277778, 0.000277778] # 1 arcsec/pixel
wcs.wcs.cunit = ["deg", "deg"]
# Create a synthetic Gaussian source
rng = np.random.default_rng(1)
image = rng.normal(100.0, 10.0, (ny, nx)) # sky background + noise
yy, xx = np.mgrid[:ny, :nx]
# Place a point source at pixel (256, 256) with peak 5000 counts
sigma = 2.5
image += 5000.0 * np.exp(-((xx - 256) ** 2 + (yy - 256) ** 2) / (2 * sigma ** 2))
# Write to FITS
header = wcs.to_header()
header["EXPTIME"] = (300.0, "Exposure time [s]")
header["FILTER"] = ("r", "Filter name")
hdu = fits.PrimaryHDU(data=image.astype(np.float32), header=header)
hdu.writeto("synthetic_field.fits", overwrite=True)
print("Saved synthetic_field.fits")
from astroquery.gaia import Gaia
# Download Gaia DR3 stars brighter than G=15 around M51
query = """
SELECT TOP 200
source_id, ra, dec, parallax, pmra, pmdec,
phot_g_mean_mag, phot_bp_mean_mag, phot_rp_mean_mag
FROM gaiadr3.gaia_source
WHERE CONTAINS(
POINT('ICRS', ra, dec),
CIRCLE('ICRS', 202.4696, 47.1952, 0.1)
) = 1
AND phot_g_mean_mag < 15
ORDER BY phot_g_mean_mag ASC
"""
job = Gaia.launch_job(query)
result = job.get_results()
print(f"Retrieved {len(result)} Gaia DR3 stars")
print(result["source_id", "ra", "dec", "phot_g_mean_mag"][:5])
from astropy.io import fits
from astropy.table import Table
import numpy as np
# Example: Chandra source catalogue (or any BinTable FITS)
# hdul = fits.open("chandra_sources.fits")
# For demonstration, create a mock BinTable
n = 100
rng = np.random.default_rng(2)
col_ra = fits.Column(name="RA", format="D", array=rng.uniform(0, 360, n))
col_dec = fits.Column(name="DEC", format="D", array=rng.uniform(-90, 90, n))
col_flux = fits.Column(name="FLUX", format="E",
array=rng.exponential(1e-13, n).astype(np.float32))
col_name = fits.Column(name="NAME", format="20A",
array=[f"CXO J{i:04d}" for i in range(n)])
hdu_table = fits.BinTableHDU.from_columns([col_ra, col_dec, col_flux, col_name])
hdu_table.name = "SOURCES"
primary = fits.PrimaryHDU()
hdul = fits.HDUList([primary, hdu_table])
hdul.writeto("mock_sources.fits", overwrite=True)
# Read back as an astropy Table (cleanest interface)
tbl = Table.read("mock_sources.fits", hdu="SOURCES")
print(tbl[:5])
print(f"Brightest source: {tbl[np.argmax(tbl['FLUX'])]['NAME']}")
astroquery times out or returns empty resultsVizier.TIMEOUT = 60 before the query.from astropy.table import Table
import os
CACHE = "gaia_cache.fits"
if os.path.exists(CACHE):
gaia = Table.read(CACHE)
else:
# ... run Vizier query ...
gaia.write(CACHE, overwrite=True)
WCS warnings about non-standard keywordsimport warnings
from astropy.wcs import FITSFixedWarning
warnings.filterwarnings("ignore", category=FITSFixedWarning)
photutils source detection finds too many / too few sourcesAdjust the detection threshold. A typical starting point:
from astropy.stats import sigma_clipped_stats
_, _, std = sigma_clipped_stats(data, sigma=3.0)
finder = DAOStarFinder(fwhm=3.0, threshold=5.0 * std) # raise threshold to reduce
finder = DAOStarFinder(fwhm=3.0, threshold=3.0 * std) # lower to detect fainter
Also ensure the PSF FWHM parameter matches the actual image seeing in pixels.
healpy crashes on import (macOS)On Apple Silicon, install via conda-forge rather than pip:
conda install -c conda-forge healpy
data = hdul[0].data
if data.dtype.byteorder not in ("=", "|"):
data = data.byteswap().newbyteorder()
"""
End-to-end aperture photometry pipeline:
1. Download SDSS r-band image of M51 via SkyView
2. Read FITS, solve WCS
3. Detect sources with DAOStarFinder
4. Perform aperture photometry
5. Cross-match detected sources against Gaia DR3
6. Compute instrumental zero-point from Gaia G magnitudes
7. Save a calibrated catalog to FITS
"""
import numpy as np
import astropy.units as u
from astropy.io import fits
from astropy.wcs import WCS
from astropy.table import Table
from astropy.coordinates import SkyCoord
from astropy.stats import sigma_clipped_stats
from astroquery.skyview import SkyView
from astroquery.vizier import Vizier
from photutils.detection import DAOStarFinder
from photutils.aperture import CircularAperture, CircularAnnulus, aperture_photometry
# --- 1. Download image ---
print("Downloading SDSS r-band image …")
imgs = SkyView.get_images("M51", survey=["SDSSr"], pixels=512)
hdul = imgs[0]
data = hdul[0].data.astype(float)
header = hdul[0].header
wcs = WCS(header)
exptime = float(header.get("EXPTIME", 53.91)) # SDSS typical 53.91 s
hdul.writeto("m51_sdssr.fits", overwrite=True)
# --- 2. Background subtraction ---
_, bkg_med, bkg_std = sigma_clipped_stats(data, sigma=3.0)
data_sub = data - bkg_med
# --- 3. Source detection ---
finder = DAOStarFinder(fwhm=3.5, threshold=5.0 * bkg_std)
sources = finder(data_sub)
print(f"Detected {len(sources)} sources")
# --- 4. Aperture photometry ---
positions = np.column_stack([sources["xcentroid"], sources["ycentroid"]])
aper = CircularAperture(positions, r=5.0)
ann = CircularAnnulus(positions, r_in=7.0, r_out=10.0)
phot = aperture_photometry(data_sub, [aper, ann])
sky_bkg = (phot["aperture_sum_1"] / ann.area) * aper.area
flux = phot["aperture_sum_0"] - sky_bkg
flux_rate = flux / exptime
# --- 5. WCS -> sky coordinates for each source ---
sky_pos = wcs.pixel_to_world(sources["xcentroid"], sources["ycentroid"])
sources["ra"] = sky_pos.ra.deg
sources["dec"] = sky_pos.dec.deg
# --- 6. Cross-match against Gaia DR3 ---
field_centre = wcs.pixel_to_world(data.shape[1] / 2, data.shape[0] / 2)
v = Vizier(columns=["Source", "RA_ICRS", "DE_ICRS", "Gmag"], row_limit=2000)
gaia_result = v.query_region(field_centre, radius=5 * u.arcmin,
catalog="I/355/gaiadr3")
gaia = gaia_result[0] if gaia_result else None
zp = 25.0 # fallback zero-point
if gaia is not None:
det_coords = SkyCoord(ra=sources["ra"] * u.deg, dec=sources["dec"] * u.deg)
gaia_coords = SkyCoord(ra=gaia["RA_ICRS"] * u.deg, dec=gaia["DE_ICRS"] * u.deg)
idx, sep, _ = det_coords.match_to_catalog_sky(gaia_coords)
good = sep < 1.0 * u.arcsec
if good.sum() > 5:
inst_mag = -2.5 * np.log10(np.abs(flux_rate[good]))
gaia_mag = gaia["Gmag"][idx[good]].data.astype(float)
zp = float(np.median(gaia_mag - inst_mag))
print(f"Photometric zero-point (Gaia G): {zp:.3f} mag")
# --- 7. Build calibrated catalog ---
mag_cal = zp - 2.5 * np.log10(np.abs(flux_rate))
catalog = Table()
catalog["id"] = sources["id"]
catalog["ra"] = sources["ra"]
catalog["dec"] = sources["dec"]
catalog["flux"] = flux_rate
catalog["mag"] = mag_cal
catalog.sort("mag")
catalog.write("m51_catalog.fits", overwrite=True)
print("Saved m51_catalog.fits with", len(catalog), "sources")
print(catalog["id", "ra", "dec", "mag"][:10])
"""
CMB power spectrum analysis:
1. Simulate a CMB temperature map with realistic Planck beam
2. Apply a Galactic mask (|b| < 20 deg)
3. Compute pseudo-Cl power spectrum
4. Correct for the beam and pixel window functions
5. Plot the final power spectrum
"""
import numpy as np
import healpy as hp
import matplotlib.pyplot as plt
# --- Parameters ---
nside = 256
lmax = 2 * nside
fwhm_deg = 0.5 # Planck 143 GHz beam ~ 7 arcmin; use 30 arcmin here
fwhm_rad = np.radians(fwhm_deg)
npix = hp.nside2npix(nside)
ells = np.arange(lmax + 1)
# --- 1. Build input power spectrum (Planck 2018-like toy model) ---
Cl_in = np.zeros(lmax + 1)
Cl_in[2:] = 2.7e-9 / (ells[2:] * (ells[2:] + 1) / (2 * np.pi)) * (
1 + 5 * np.exp(-((ells[2:] - 200) / 50) ** 2)
+ 2 * np.exp(-((ells[2:] - 540) / 40) ** 2)
)
# --- 2. Simulate map with beam ---
np.random.seed(7)
beam_window = hp.gauss_beam(fwhm_rad, lmax=lmax)
alm_true = hp.synalm(Cl_in, lmax=lmax)
alm_beam = hp.almxfl(alm_true, beam_window)
cmb_map = hp.alm2map(alm_beam, nside=nside, verbose=False)
# --- 3. Build Galactic mask (cut |b| < 20 deg) ---
theta, phi = hp.pix2ang(nside, np.arange(npix))
b_rad = np.pi / 2 - theta # colatitude -> Galactic latitude proxy
mask = np.abs(b_rad) > np.radians(20)
fsky = mask.mean()
print(f"Sky fraction after masking: f_sky = {fsky:.3f}")
masked_map = cmb_map * mask
# --- 4. Pseudo-Cl estimation ---
alm_masked = hp.map2alm(masked_map, lmax=lmax)
Cl_pseudo = hp.alm2cl(alm_masked)
# --- 5. Beam and pixel window correction ---
pixel_win = hp.pixwin(nside)[:lmax + 1]
correction = (beam_window ** 2) * (pixel_win ** 2) * fsky
correction = np.where(correction > 1e-30, correction, 1e-30)
Cl_corr = Cl_pseudo / correction
# --- 6. Plot ---
fig, ax = plt.subplots(figsize=(10, 5))
norm = ells * (ells + 1) / (2 * np.pi)
ax.plot(ells[2:], norm[2:] * Cl_in[2:], lw=2, label="Input $C_\\ellquot;, color="k")
ax.plot(ells[2:], norm[2:] * Cl_corr[2:], lw=1.5, label="Recovered (beam-corrected)",
color="tab:blue", alpha=0.8)
ax.set_xlabel(r"Multipole $\ellquot;)
ax.set_ylabel(r"$\ell(\ell+1)C_\ell / 2\pi$ [$\mu$K$^2$]")
ax.set_xlim(2, lmax)
ax.set_yscale("log")
ax.legend()
ax.set_title(f"CMB Power Spectrum (nside={nside}, mask f_sky={fsky:.2f})")
fig.tight_layout()
plt.savefig("cmb_analysis.png", dpi=150)
print("Saved cmb_analysis.png")