fitscube-rs 1.2.2

Combine single-frequency/single-time FITS images into a FITS cube (Rust port of fitscube)
Documentation 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114

---
file_format: mystnb
kernelspec:
  name: python3
  display_name: Python 3
---

# Python quickstart

This page is an executable notebook: every cell below is re-run on each docs
build, so the outputs are guaranteed to match the current release.

## Building a set of single-channel images

fitscube-rs combines many single-frequency (or single-time) FITS images into a
single FITS cube. To demonstrate, we first write out a handful of synthetic
single-channel images, each with its own frequency in the header.

```{code-cell} ipython3
import tempfile
from pathlib import Path

import numpy as np
from astropy.io import fits

workdir = Path(tempfile.mkdtemp())

freqs_hz = np.array([1.0e9, 1.1e9, 1.2e9, 1.3e9])  # evenly spaced in frequency
n = 64

paths = []
for i, freq in enumerate(freqs_hz):
    data = np.full((n, n), float(i), dtype=np.float32)
    header = fits.Header()
    header["CRVAL3"] = freq
    header["CTYPE3"] = "FREQ"
    header["CUNIT3"] = "Hz"
    path = workdir / f"chan_{i:04d}.fits"
    fits.writeto(path, data[np.newaxis, :, :], header, overwrite=True)
    paths.append(path)

sorted(p.name for p in paths)
```

## Combining images into a cube

{func}`~fitscube_rs.combine_fits` stacks the per-channel images along a new
spectral (or time) axis, building the WCS for that axis from the per-image
headers. The resulting cube is written to disk:

```{code-cell} ipython3
from fitscube_rs import combine_fits

out_cube = workdir / "cube.fits"
freqs = combine_fits(
    [str(p) for p in sorted(paths)],
    str(out_cube),
    overwrite=True,
)

with fits.open(out_cube) as hdul:
    print("cube shape:", hdul[0].data.shape)
    print("frequencies (Hz):", freqs)
```

The frequency axis is detected as evenly spaced, so the cube header carries a
linear `CRVAL3`/`CDELT3` description. When the spacing is uneven, fitscube-rs
instead writes an explicit per-plane frequency table so no information is lost
(see [Algorithm background](algorithms.md)).

## Extracting a plane

{func}`~fitscube_rs.extract_plane_from_cube` pulls a single plane back out of a
cube — the inverse of combining — which is handy for inspecting one channel or
feeding a downstream tool that expects a 2D image:

```{code-cell} ipython3
from fitscube_rs import extract_plane_from_cube

plane_path = workdir / "chan2.fits"
extract_plane_from_cube(
    str(out_cube), channel_index=2, output_path=str(plane_path), overwrite=True
)

with fits.open(plane_path) as hdul:
    print("plane shape:", hdul[0].data.shape)
    print("plane value:", float(np.nanmean(hdul[0].data)))
```

## Per-channel beams

Radio images often carry a restoring beam (`BMAJ`/`BMIN`/`BPA`) that varies
from channel to channel. When the input images have differing beams, the
combined cube records them in a CASA-style `BEAMS` binary-table extension and
sets `CASAMBM=T` in the primary header, matching the convention used by
`astropy` and CASA:

```{code-cell} ipython3
with fits.open(out_cube) as hdul:
    has_beams = any(h.name == "BEAMS" for h in hdul)
    print("primary HDUs:", [h.name or "PRIMARY" for h in hdul])
    print("multi-beam table present:", has_beams)
```

## Working with the CLI

For batch jobs (hundreds of channels, large mosaics) prefer the
[`fitscubers` CLI](cli.md) — it parallelises the read/stack across workers and
streams planes to disk so peak memory stays bounded:

```sh
fitscubers combine chan_*.fits cube.fits --overwrite
fitscubers extract cube.fits --channel-index 2 --output-path chan2.fits
```