Working with FITS Files

beginnerPythonastropyhealpymatplotlibnumpyFITSHDUImage HDUBinTableHEALPixWCSall-sky mapAitoff-HammerMollweideastropyhealpy

Working with FITS Files

FITS (Flexible Image Transport System) is the standard file format in astronomy. A FITS file is made up of one or more Header/Data Units (HDUs), each consisting of a plain-text header (keyword–value pairs) and an optional data block. The data block can be an N-dimensional array (Image HDU) or a structured table (BinTable or ASCII Table HDU). This playbook shows how to inspect an unfamiliar FITS file, identify what type of data it contains, and visualize it.

Requirements

pip install astropy matplotlib numpy healpy

Step 1: Inspect the HDU structure

Start by listing all HDUs without reading the full data into memory.

from astropy.io import fits

url = "https://..."  # local path or remote URL

fits.info(url)

fits.info prints a table of HDU index, name, type, dimensions, and format. This alone usually tells you whether the file contains images, tables, or both.

For a closer look at each HDU—especially to check WCS keywords or column names—open the file and iterate:

from astropy.io import fits

with fits.open(url) as hdul:
    for i, hdu in enumerate(hdul):
        hdr = hdu.header
        print(f"\n--- HDU {i}: {hdu.name} ({type(hdu).__name__}) ---")
        print(f"  NAXIS  = {hdr.get('NAXIS')}")
        for kw in ("NAXIS1", "NAXIS2", "NAXIS3",
                   "CTYPE1", "CTYPE2", "CUNIT1", "CUNIT2",
                   "BUNIT", "PIXTYPE", "ORDERING", "NSIDE", "COORDSYS"):
            if kw in hdr:
                print(f"  {kw:10s} = {hdr[kw]!r}")
        if hasattr(hdu, "columns"):
            print(f"  columns = {hdu.columns.names}")

After running this, you will typically fall into one of the three cases below.

Case A: Image HDU

An Image HDU has NAXIS >= 2 and no column definitions. The data is a NumPy array accessible via hdu.data.

import numpy as np
import matplotlib.pyplot as plt
from astropy.io import fits
from astropy.wcs import WCS

with fits.open(url) as hdul:
    hdu = hdul[0]          # adjust index as needed
    data = hdu.data
    wcs  = WCS(hdu.header)

print(data.shape, data.dtype)

# Replace non-positive values with NaN for log scaling
img = np.where(data > 0, data, np.nan)

fig, ax = plt.subplots(figsize=(8, 6))
im = ax.imshow(np.log10(img), origin="lower", cmap="afmhot", aspect="auto")
fig.colorbar(im, ax=ax, label="log₁₀ (value)")
plt.tight_layout()
plt.close()

All-sky projection images are a common sub-case. When CTYPE1 / CTYPE2 end in -AIT (Aitoff-Hammer) or -MOL (Mollweide), the 2D image is a full-sky map in that projection, stored in Galactic or Equatorial coordinates as indicated by CTYPE. The pixel layout directly matches the projection; no additional reprojection is needed to display it with imshow. See the MAXI All-Sky Map Visualization playbook for a concrete example with this structure.

Case B: HEALPix BinTable HDU

HEALPix stores full-sky data as a flat table: one row per sky pixel, ordered by the HEALPix pixel index. Recognition cues in the header:

Keyword	Typical value	Meaning
`PIXTYPE`	`'HEALPIX'`	Confirms this is a HEALPix map
`ORDERING`	`'RING'` or `'NESTED'`	Pixel ordering scheme
`NSIDE`	64, 128, 512, …	Resolution parameter (Npix = 12·NSIDE²)
`COORDSYS`	`'G'`, `'E'`, `'C'`	Galactic, Ecliptic, or Equatorial

Reading the table and plotting with healpy

import numpy as np
import healpy as hp
from astropy.io import fits
from astropy.table import Table

with fits.open(url) as hdul:
    # Find the HEALPix extension (often HDU 1)
    for i, hdu in enumerate(hdul):
        if hdu.header.get("PIXTYPE") == "HEALPIX":
            hdr   = hdu.header
            table = Table(hdu.data)
            break

nside    = hdr["NSIDE"]
ordering = hdr.get("ORDERING", "RING").upper()
nest     = (ordering == "NESTED")

print("Columns:", table.colnames)

# Pick the intensity column — adapt the name to your file
col   = table.colnames[0]
m     = table[col].data.astype(float)
m[m <= 0] = hp.UNSEEN    # healpy ignores UNSEEN pixels

hp.mollview(m, nest=nest, title=col, unit=hdr.get("TUNIT1", ""),
            coord=["G"] if hdr.get("COORDSYS") == "G" else None)
hp.graticule()
import matplotlib.pyplot as plt
plt.close()

Shortcut with `healpy.read_map`

healpy.read_map can read a standard HEALPix FITS table in one call. It returns the map array and automatically interprets the ORDERING keyword.

import healpy as hp

m, hdr = hp.read_map(url, h=True, dtype=float)
hdr = dict(hdr)

hp.mollview(m, title="HEALPix map", unit=hdr.get("TUNIT1", ""))
hp.graticule()
import matplotlib.pyplot as plt
plt.close()

Note: healpy.read_map reads the first column by default. Pass field=N to read column N.

Case C: BinTable HDU (events or catalogs)

A BinTable that is not a HEALPix map typically holds event lists (one row = one photon or particle) or source catalogs. Loading is straightforward:

from astropy.io import fits
from astropy.table import Table

with fits.open(url) as hdul:
    table = Table(hdul["EVENTS"].data)   # or use the HDU index

print(table.colnames)
print(table[:5])

For working examples of event-list FITS files, see the MAXI All-Sky Map from Event Data playbook and the Hinode/EIS Level-0 playbook.

Working with FITS Files