FITSCUBE

From the wsclean docs:

WSClean does not output these images in a normal “imaging cube” like CASA does, i.e., a single fits file with several images in it. For now I’ve decided not to implement this (one of the reasons for this is that information about the synthesized beam is not properly stored in a multi-frequency fits file). One has of course the option to combine the output manually, e.g. with a simple Python script.

This is a simple Python script to combine (single-frequency) FITS images manually.

Current assumptions:

• All files have the same WCS
• All files have the same shape / pixel grid
• Frequency is either a WCS axis or in the REFFREQ header keyword
• All the relevant information is in the first header of the first image

Installation

Install from PyPI (stable):

pip install fitscube

Or, install from this git repo (latest):

pip install git+https://github.com/AlecThomson/fitscube.git

To install and use the uvloop Async runner do:

pip install fitscube[uvloop]

Usage

Command line:

❯ fitscube -h
usage: fitscube [-h] {combine,extract} ...

Tooling to create fitscubes

positional arguments:
  {combine,extract}
    combine          Combine FITS images together into a cube
    extract          Extract a plane from an existing cube

options:
  -h, --help         show this help message and exit

❯ fitscube combine -h
usage: fitscube combine [-h] [-o] [--create-blanks] [--time-domain]
                        [--spec-file SPEC_FILE | --specs SPECS [SPECS ...] | --ignore-spec] [-v]
                        [--max-workers MAX_WORKERS] [--bounding-box] [--invalidate-zeros]
                        [--no-blank-zero-beams] [--floating {8,16,32,64}] [--compress]
                        [--compress-method {gzip,pgzip}]
                        file_list [file_list ...] out_cube

positional arguments:
  file_list             List of FITS files to combine (in frequency or time order)
  out_cube              Output FITS file

options:
  -h, --help            show this help message and exit
  -o, --overwrite       Overwrite output file if it exists
  --create-blanks       Try to create a blank cube with evenly spaced frequencies
  --time-domain         Flag for constructing a time-domain cube
  --spec-file SPEC_FILE
                        File containing frequencies in Hz or times in MJD s (if --time-domain ==
                        True)
  --specs SPECS [SPECS ...]
                        List of frequencies or times in Hz or MJD s respectively
  --ignore-spec         Ignore frequency or time information and just stack (probably not what you
                        want)
  -v, --verbosity       Increase output verbosity
  --max-workers MAX_WORKERS
                        Maximum number of workers to use for concurrent processing
  --bounding-box        Attempt to consider padded images when creating the cube. Requires an
                        extract read of the input data.
  --invalidate-zeros    Set pixels whose values are exactly zero to NaNs
  --no-blank-zero-beams
                        Keep images whose restoring beam is exactly zero. By default such images
                        are blanked with NaNs, as a zero beam means no PSF was fitted for that
                        plane (e.g. a wsclean model image planted by -fit-spectral-pol)
  --floating {8,16,32,64}
                        The number of floating point bits to use in the out cube. If None the
                        input data precision is used.
  --compress            Gzip-compress the output cube
  --compress-method {gzip,pgzip}
                        Compression backend for --compress. 'gzip' is stdlib and needs nothing
                        extra. 'pgzip' parallelizes across --max-workers threads but requires the
                        'pgzip' extra: pip install fitscube[pgzip]

❯ fitscube extract -h
usage: fitscube extract [-h] [--channel-index CHANNEL_INDEX] [--hdu-index HDU_INDEX] [-v] [--overwrite] [--output-path OUTPUT_PATH] fits_cube

positional arguments:
  fits_cube             The cube to extract a plane from

options:
  -h, --help            show this help message and exit
  --channel-index CHANNEL_INDEX
                        The channel to extract
  --hdu-index HDU_INDEX
                        The HDU index of the data card containing the cube data
  -v, --verbosity       Increase output verbosity
  --overwrite           overwrite the output file, if it exists.
  --output-path OUTPUT_PATH
                        The name of the new output file. If not specified it is generated from the fits cube name

Python:

from pathlib import Path

from fitscube import combine_fits

file_list = list(Path().glob("*.fits"))

frequencies = combine_fits(
    file_list
)

Zero-sized beams

wsclean writes BMAJ = BMIN = 0 into an image when no PSF was fitted for that plane. This happens with -fit-spectral-pol, where the channels that were not imaged directly are populated with the model image instead. Those planes look like data, but they are not comparable to the rest of the cube and will produce nasty effects downstream.

By default fitscube blanks any such plane with NaNs and logs a warning. Pass --no-blank-zero-beams (or blank_zero_beams=False in Python) to keep the image data as it is.

Either way, a zero beam is never written into the beam table as a literal zero - it gets the same np.finfo(np.float32).tiny sentinel that already marks a NaN PSF, which is what keeps CASA happy. This is applied per beam, so a zero BPA on an otherwise valid beam is left alone.

Convolving to a common resolution

See RACS-Tools.

License

MIT

Contributing

Contributions are welcome. Please open an issue or pull request.