[Python] ds9 and python: how to create RGB healpix fits image?

This python function creates 3 fits files, one per RGB channel:

# HealpixFullArray is a numpy array of shape (N,3) and dtype uint8.
# N: number of healpix pixels (in my case corresponding to nside 1024)
def save_fits(healpix_map: HealpixFullArray, output_base: Path) -> None:
for i, channel in enumerate(["red", "green", "blue"]):
output_path = output_base.parent / f"{output_base.stem}_{channel}.fits"
healpy.write_map(
output_path,
healpix_map[:, i],
nest=True,
dtype=healpix_map.dtype,
overwrite=True,
)
logger.info(f"Saved {channel} channel to {output_path}")


Following ds9 instructions (https://ds9.si.edu/doc/ref/file.html#FITSRGB), I can open a corresponding RGB image:

ds9 -rgb -red foo.fits -green bar.fits -blue wow.fits # rgb image from 3 fits images

But I ideally I would like to create a single file containing all RGB data, i.e. that I could open via on of these commands:

ds9 -rgbimage rgb.fits # load rgb image consisting of one fits file with 3 image exts

or

ds9 -rgbcube cube.fits # load rgb image consisting of one fits data cube

No code I tried so far had any success.

Continue reading...