## 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.

```python
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:

```python
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`.

```python
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](https://darts.isas.jaxa.jp/analysis-playbooks/maxi-allskymap-products) 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

```python
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.

```python
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:

```python
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](https://darts.isas.jaxa.jp/analysis-playbooks/maxi-allskymap-from-events) and the [Hinode/EIS Level-0 playbook](https://darts.isas.jaxa.jp/analysis-playbooks/hinode-eis-level0).

---

### References

- [FITS Standard (NASA/GSFC)](https://fits.gsfc.nasa.gov/fits_standard.html)
- [astropy.io.fits documentation](https://docs.astropy.org/en/stable/io/fits/)
- [astropy.wcs documentation](https://docs.astropy.org/en/stable/wcs/)
- [HEALPix primer](https://healpix.sourceforge.io/documentation.php)
- [healpy documentation](https://healpy.readthedocs.io/)