Astrophysics Data Guide

A skill for processing and analyzing astronomical data using standard astrophysics tools. Covers FITS file handling, coordinate transformations, photometric analysis, spectral analysis, catalog cross-matching, and accessing major sky survey archives.

Astronomical Data Formats

FITS Files

FITS (Flexible Image Transport System) is the standard data format in astronomy:

from astropy.io import fits
import numpy as np

def inspect_fits(filepath: str) -> dict:
    """
    Inspect the structure of a FITS file.
    Returns information about each HDU (Header/Data Unit).
    """
    with fits.open(filepath) as hdul:
        info = []
        for i, hdu in enumerate(hdul):
            entry = {
                "index": i,
                "name": hdu.name,
                "type": type(hdu).__name__,
            }
            if hdu.data is not None:
                entry["shape"] = hdu.data.shape
                entry["dtype"] = str(hdu.data.dtype)
            if hasattr(hdu, "columns") and hdu.columns is not None:
                entry["columns"] = [c.name for c in hdu.columns]
            info.append(entry)
        return {"filename": filepath, "n_hdus": len(hdul), "hdus": info}

def read_fits_image(filepath: str, hdu_index: int = 0) -> tuple:
    """Read a FITS image and its WCS (World Coordinate System)."""
    from astropy.wcs import WCS

    with fits.open(filepath) as hdul:
        data = hdul[hdu_index].data
        header = hdul[hdu_index].header
        wcs = WCS(header)

    return data, wcs, header

Working with FITS Tables

from astropy.table import Table

def read_fits_catalog(filepath: str, hdu: int = 1) -> Table:
    """Read a FITS binary table extension as an Astropy Table."""
    catalog = Table.read(filepath, hdu=hdu)
    print(f"Catalog: {len(catalog)} objects, {len(catalog.columns)} columns")
    print(f"Columns: {catalog.colnames}")
    return catalog

Coordinate Systems

Astronomical Coordinate Transformations

from astropy.coordinates import SkyCoord, EarthLocation, AltAz
from astropy.time import Time
import astropy.units as u

def coordinate_transforms(ra_deg: float, dec_deg: float) -> dict:
    """
    Transform between astronomical coordinate systems.
    ra_deg, dec_deg: right ascension and declination in degrees (ICRS/J2000)
    """
    coord = SkyCoord(ra=ra_deg * u.degree, dec=dec_deg * u.degree, frame="icrs")

    return {
        "icrs": {
            "ra": coord.ra.to_string(unit=u.hourangle, precision=2),
            "dec": coord.dec.to_string(unit=u.degree, precision=2),
        },
        "galactic": {
            "l": round(coord.galactic.l.degree, 4),
            "b": round(coord.galactic.b.degree, 4),
        },
        "ecliptic": {
            "lon": round(coord.geocentricmeanecliptic.lon.degree, 4),
            "lat": round(coord.geocentricmeanecliptic.lat.degree, 4),
        },
    }

def compute_altaz(ra_deg: float, dec_deg: float,
                   obs_time: str, location: tuple) -> dict:
    """
    Compute altitude and azimuth for a target from a given location and time.
    location: (latitude_deg, longitude_deg, elevation_m)
    """
    target = SkyCoord(ra=ra_deg * u.degree, dec=dec_deg * u.degree)
    time = Time(obs_time)
    loc = EarthLocation(
        lat=location[0] * u.degree,
        lon=location[1] * u.degree,
        height=location[2] * u.m,
    )
    altaz_frame = AltAz(obstime=time, location=loc)
    altaz = target.transform_to(altaz_frame)

    return {
        "altitude_deg": round(altaz.alt.degree, 2),
        "azimuth_deg": round(altaz.az.degree, 2),
        "airmass": round(altaz.secz.value, 3) if altaz.alt.degree > 0 else None,
        "is_observable": altaz.alt.degree > 10,
    }