rascil_imager

rascil_imager is a command line app written using RASCIL. It supports three ways of making an image:

  • invert: Inverse Fourier Transform of the visibilities to make a dirty image (or point spread function)

  • cip: The SKA Continuum Imaging Pipeline.

  • ical: The SKA Iterative Calibration Pipeline (ICAL)

Notable features:

  • Reads a CASA MeasurementSet and writes FITS files

  • Image size can be a composite of 2, 3, 5

  • Distribute processing across processors using Dask

  • Multi Frequency Synthesis Multiscale CLEAN available, also with distribution of CLEAN over facets

  • Distribution of restoration over facets

  • Wide field imaging using the fast and accurate nifty gridder

  • Modelling of bright sources by fitting with sub-pixel locations

  • Selfcalibration available for atmosphere (T), complex gains (G), and bandpass (B)

  • Selection of data by uv range and r range (where r is the distance of station/dish from array centre

CLI arguments are grouped:

  • --mode prefixed parameters controls which algorithm is run.

  • --imaging prefixed parameters control the details of the imaging such as number of pixels, cellsize

  • --clean prefixed parameters control the clean deconvolutions (active only for modes cip and ical)

  • --calibration prefixed parameters control the calibration in the ICAL pipeline. (active only for mode ical)

  • --dask prefixed parameters control the use of Dask/rsexecute for distributing the processing

MeasurementSet ingest

Although a CASA MeasurementSet can hold heterogeneous observations, identified by data descriptors. rascil-imager can only process identical data descriptors from a MS. The number of channels and polarisation must be the same.

Each selected data descriptor is optionally split into a number of channels optionally averaged and placed into one Visibility.

For example, using the arguments:

--ingest_msname SNR_G55_10s.calib.ms --ingest_dd 0 1 2 3 --ingest_vis_nchan 64 \
--ingest_chan_per_vis 8 --ingest_average_vis True

will read data descriptors 0, 1, 2, 3, each of which has 64 channels. Each set of 64 channels are split into blocks of 8 and averaged. We thus end up with 32 separate datasets in RASCIL, each of which is a Visibility and has 1 channel, for a total of 32 channels. If the argument --ingest_average_vis is set to False, each Visibility has eight channels, for a total of 256 channels.

Selection

rascil_imager supports selection of data by uv range --imaging_uvmin --imaging_uvmax, and by dish/station based on distance from the array centre --imaging_rmin --imaging_rmax

Imaging

To make an image from visibilities or to predict visibilities from a model, it is necessary to use a gridder. Nifty gridder (https://gitlab.mpcdf.mpg.de/ift/nifty_gridder) is currently the best gridder to use in RASCIL. It is written in c and uses OpenMP to distribute the processing across multiple threads. The Nifty Gridder uses an improved wstacking algorithm uses many fewer w-planes than w stacking or w projection. It is not necessary to explicitly set the number of w-planes.

The gridder is set by the --imaging_context argument. The default, --imaging_context ng is the Nifty Gridder.

CLEAN

rascil-imager supports Hogbom CLEAN, MultiScale CLEAN, and Multi-Frequency Synthesis MultiScale Clean (also known as MMCLEAN). The first two work independently on different frequency channels, while MMClean works jointly cross all channels using a Taylor Series expansion in frequency for the emission.

The clean methods support a number of processing speed enhancements:

  • The multi-frequency-synthesis CLEAN works by fitting a Taylor series in frequency. The --ingest_chan_per_vis argument controls the aggregation of channels in the MeasurementSet to form image planes for the CLEAN. Within a Visibility the different channels are gridded together to form one image. Each image is then used in the mmclean algorithm. For example, a data set may have 256 channels spread over 4 data descriptors. We can split these into 32 BlockVisibilities and then run the mmclean over these 32 channels.

  • Only a limited central region of the PSF will be subtracted during the minor cycles.

  • The cleaning may be partitioned into overlapping facets, each of which is cleaned independently, and then merged with neighbours using a taper function. This works well for fields of compact sources but is likely to not perform well for extended emission.

  • The restoration may be distributed via subimages. This requires that the subimages have significant overlap such that the clean beam can fit within the overlap area.

Bright compact sources can optionally be represented by discrete components instead of pixels.

  • --clean_component_threshold 0.5 All sources > 0.5 Jy to be fitted

  • --clean_component_method fit non-linear last squares algorithm to find source parameters

The skymodel written at the end of processing will include both the image model and the skycomponents.

Polarisation

The polarisation processing behaviour is controlled by --image_pol.

  • --image_pol stokesI will image only the I Stokes parameter

  • --image_pol stokesIQUV will image all Stokes parameters I, Q, U, V

Note that the combination of MM CLEAN and stokesIQUV imaging is not likely to be meaningful.

Self-calibration

rascil-imager supports self-calibration as part of the imaging. At the end of each major cycle a calibration solution and application may optionally be performed.

Calibration uses the Hamaker Bregman Sault formalism with the following Jones matrices supported: T (Atmospheric phase), G (Electronics gain), B - (Bandpass).

An example consider the arguments:

calibration_T_first_selfcal = 2
calibration_T_phase_only = True
calibration_T_timeslice = None
calibration_G_first_selfcal = 5
calibration_G_phase_only = False
calibration_G_timeslice = 1200.0
calibration_B_first_selfcal = 8
calibration_B_phase_only = False
calibration_B_timeslice = 1.0e5
calibration_global_solution = True
calibration_calibration_context = "TGB"

These will perform a phase only solution of the T term after the second major cycle for every integration, solution of G after 5 major cycles with timescale of 1200s, and solution of B after 8 major cycles, integrating across all frequencies where appropriate. Note, that T and G terms are averages across frequency.

SkyModel in ICAL

When running rascil_imager in mode ical, optionally, an initial SkyModel can be used. To do this, set --use_initial_skymodel to True. The SkyModel is made up of model images (created based on input BlockVisibilities), and SkyComponents. The kind of SkyComponent(s) to use in the initial SkyModel is controlled by the --input_skycomponent_file and --num_bright_sources arguments:

  1. If no input file is provided, a point source at the phase centre, with brightness of 1 Jy is used as the component.

  2. If either an HDF file or a TXT file is provided, the components are read from the file.

    1. if --num_bright_sources is left as None, all of the components are used for the SkyModel

    2. if --num_bright_sources is an integer n (n>0), then n number of the brightest components are used for the SkyModel

This SkyModel is then overwritten during the remaining cycles of the run.

By default, --use_initial_skymodel is set to False, and hence no initial SkyModel is used.

In addition, you can decide whether to reset the initial skymodel after first calibration, or not, by setting the --calibration_reset_skymodel either to True or False.

Dask

Dask is used to distribute processing across multiple cores or nodes. The setup and execution of a set of workers is controlled by a scheduler. By default, rascil uses the process scheduler which sets up a number of processes each with a number of threads. If the host has 16 cores, the set up will be 4 processes each with 4 threads for a total of 16 Dask workers.

For distribution across a cluster, the Dask distributed processor is required. See RASCIL and DASK for more details.

Example script

The following runs the cip on a data set from the CASA examples:

#!/bin/bash
# Run this in the directory containing SNR_G55_10s.calib.ms
# (The dataset can be downloaded at
# http://casa.nrao.edu/Data/EVLA/SNRG55/SNR_G55_10s.calib.tar.gz)
python $RASCIL/rascil/apps/rascil_imager.py --mode cip \
--ingest_msname SNR_G55_10s.calib.ms --ingest_dd 0 1 2 3 --ingest_vis_nchan 64 \
--ingest_chan_per_vis 8 --ingest_average_vis True \
--imaging_npixel 1280 --imaging_cellsize 3.878509448876288e-05 \
--imaging_weighting robust --imaging_robustness -0.5 \
--clean_nmajor 5 --clean_algorithm mmclean --clean_scales 0 6 10 30 60 \
--clean_fractional_threshold 0.3 --clean_threshold 0.12e-3 --clean_nmoment 5 \
--clean_psf_support 640 --clean_restored_output integrated

Command line arguments

RASCIL continuum imager

usage: rascil_imager [-h] [--mode MODE] [--logfile LOGFILE]
                     [--performance_file PERFORMANCE_FILE]
                     [--ingest_msname INGEST_MSNAME]
                     [--ingest_dd [INGEST_DD ...]]
                     [--ingest_vis_nchan INGEST_VIS_NCHAN]
                     [--ingest_chan_per_vis INGEST_CHAN_PER_VIS]
                     [--ingest_average_vis INGEST_AVERAGE_VIS]
                     [--imaging_phasecentre IMAGING_PHASECENTRE]
                     [--imaging_pol IMAGING_POL]
                     [--imaging_nchan IMAGING_NCHAN]
                     [--imaging_context IMAGING_CONTEXT]
                     [--imaging_ng_threads IMAGING_NG_THREADS]
                     [--imaging_w_stacking IMAGING_W_STACKING]
                     [--imaging_flat_sky IMAGING_FLAT_SKY]
                     [--imaging_npixel IMAGING_NPIXEL]
                     [--imaging_cellsize IMAGING_CELLSIZE]
                     [--imaging_weighting IMAGING_WEIGHTING]
                     [--imaging_robustness IMAGING_ROBUSTNESS]
                     [--imaging_gaussian_taper IMAGING_GAUSSIAN_TAPER]
                     [--imaging_dopsf IMAGING_DOPSF]
                     [--imaging_dft_kernel IMAGING_DFT_KERNEL]
                     [--imaging_uvmax IMAGING_UVMAX]
                     [--imaging_uvmin IMAGING_UVMIN]
                     [--imaging_rmax IMAGING_RMAX]
                     [--imaging_rmin IMAGING_RMIN]
                     [--perform_flagging PERFORM_FLAGGING]
                     [--flagging_strategy_name FLAGGING_STRATEGY_NAME]
                     [--calibration_reset_skymodel CALIBRATION_RESET_SKYMODEL]
                     [--calibration_T_first_selfcal CALIBRATION_T_FIRST_SELFCAL]
                     [--calibration_T_phase_only CALIBRATION_T_PHASE_ONLY]
                     [--calibration_T_timeslice CALIBRATION_T_TIMESLICE]
                     [--calibration_G_first_selfcal CALIBRATION_G_FIRST_SELFCAL]
                     [--calibration_G_phase_only CALIBRATION_G_PHASE_ONLY]
                     [--calibration_G_timeslice CALIBRATION_G_TIMESLICE]
                     [--calibration_B_first_selfcal CALIBRATION_B_FIRST_SELFCAL]
                     [--calibration_B_phase_only CALIBRATION_B_PHASE_ONLY]
                     [--calibration_B_timeslice CALIBRATION_B_TIMESLICE]
                     [--calibration_global_solution CALIBRATION_GLOBAL_SOLUTION]
                     [--calibration_context CALIBRATION_CONTEXT]
                     [--use_initial_skymodel USE_INITIAL_SKYMODEL]
                     [--input_skycomponent_file INPUT_SKYCOMPONENT_FILE]
                     [--num_bright_sources NUM_BRIGHT_SOURCES]
                     [--calibrate_with_dp3 CALIBRATE_WITH_DP3]
                     [--input_dp3_skymodel INPUT_DP3_SKYMODEL]
                     [--clean_algorithm CLEAN_ALGORITHM]
                     [--clean_use_radler CLEAN_USE_RADLER]
                     [--clean_beam CLEAN_BEAM CLEAN_BEAM CLEAN_BEAM]
                     [--clean_scales [CLEAN_SCALES ...]]
                     [--clean_nmoment CLEAN_NMOMENT]
                     [--clean_nmajor CLEAN_NMAJOR] [--clean_niter CLEAN_NITER]
                     [--clean_psf_support CLEAN_PSF_SUPPORT]
                     [--clean_gain CLEAN_GAIN]
                     [--clean_threshold CLEAN_THRESHOLD]
                     [--clean_component_threshold CLEAN_COMPONENT_THRESHOLD]
                     [--clean_component_method CLEAN_COMPONENT_METHOD]
                     [--clean_fractional_threshold CLEAN_FRACTIONAL_THRESHOLD]
                     [--clean_facets CLEAN_FACETS]
                     [--clean_overlap CLEAN_OVERLAP]
                     [--clean_taper CLEAN_TAPER]
                     [--clean_restore_facets CLEAN_RESTORE_FACETS]
                     [--clean_restore_overlap CLEAN_RESTORE_OVERLAP]
                     [--clean_restore_taper CLEAN_RESTORE_TAPER]
                     [--clean_restored_output CLEAN_RESTORED_OUTPUT]
                     [--use_dask USE_DASK] [--dask_nthreads DASK_NTHREADS]
                     [--dask_memory DASK_MEMORY]
                     [--dask_memory_usage_file DASK_MEMORY_USAGE_FILE]
                     [--dask-nodes [DASK_NODES ...]]
                     [--dask_nworkers DASK_NWORKERS]
                     [--dask_scheduler DASK_SCHEDULER]
                     [--dask_scheduler_file DASK_SCHEDULER_FILE]
                     [--dask_tcp_timeout DASK_TCP_TIMEOUT]
                     [--dask_connect_timeout DASK_CONNECT_TIMEOUT]
                     [--dask_malloc_trim_threshold DASK_MALLOC_TRIM_THRESHOLD]

Named Arguments

--mode

Processing cip | ical | invert | load

Default: “cip”

--logfile

Name of logfile (default is to construct one from msname)

--performance_file

Name of json file to contain performance information

--ingest_msname

MeasurementSet to be read

--ingest_dd

Data descriptors in MS to read (all must have the same number of channels)

Default: [0]

--ingest_vis_nchan

Number of channels in a single data descriptor in the MS

--ingest_chan_per_vis

Number of channels per vis (before any average)

Default: 1

--ingest_average_vis

Average all channels in vis?

Default: “False”

--imaging_phasecentre

Phase centre (in SkyCoord string format)

--imaging_pol

RASCIL polarisation frame for image

Default: “stokesI”

--imaging_nchan

Number of channels per image

Default: 1

--imaging_context

Imaging context i.e. the gridder used 2d | ng

Default: “ng”

--imaging_ng_threads

Number of Nifty Gridder threads to use (4 is a good choice)

Default: 4

--imaging_w_stacking

Use the improved w stacking method in Nifty Gridder?

Default: True

--imaging_flat_sky

If using a primary beam, normalise to flat sky?

Default: False

--imaging_npixel

Number of pixels in ra, dec: Should be a composite of 2, 3, 5

--imaging_cellsize

Cellsize (radians). Default is to calculate.

--imaging_weighting

Type of weighting uniform or robust or natural)

Default: “uniform”

--imaging_robustness

Robustness for robust weighting

Default: 0.0

--imaging_gaussian_taper

Size of Gaussian smoothing, implemented as taper in weights (rad)

--imaging_dopsf

Make the PSF instead of the dirty image?

Default: “False”

--imaging_dft_kernel

DFT kernel: cpu_looped | gpu_raw

--imaging_uvmax

Maximum uv (wavelengths)

--imaging_uvmin

Minimum uv (wavelengths)

--imaging_rmax

Maximum distance of dish/station from array center (wavelengths)

--imaging_rmin

Minimum distance of dish/station from array center (wavelengths)

--perform_flagging

If enabled, runs AOFlagger flagging strategy

Default: “False”

--flagging_strategy_name

Contains the name of the flagging strategy to use when perform_flagging is True. There are strategies available for different telescopes: AARTFAAC, ARECIBO, ARECIBO 305M, BIGHORNS, EVLA, JVLA, LOFAR, MWA, PARKES, PKS, ATPKSMB, WSRT. If the desired telescope is not listed here, you can use one of the strategies defined in the AOFlagger repository (https://gitlab.com/aroffringa/aoflagger/-/tree/master/data/strategies) or define a new strategy interactively using the AOFlagger rfigui (https://aoflagger.readthedocs.io/en/latest/using_rfigui.html)

Default: “generic”

--calibration_reset_skymodel

Reset the initial skymodel after initial calibration?

Default: “True”

--calibration_T_first_selfcal

First selfcal for T (complex gain). T is common to both receptors

Default: 1

--calibration_T_phase_only

Phase only solution

Default: “True”

--calibration_T_timeslice

Solution length (s) 0 means minimum

--calibration_G_first_selfcal

First selfcal for G (complex gain). G is different for the two receptors

Default: 3

--calibration_G_phase_only

Phase only solution?

Default: “False”

--calibration_G_timeslice

Solution length (s) 0 means minimum

--calibration_B_first_selfcal

First selfcal for B (bandpass complex gain). B is complex gain per frequency.

Default: 4

--calibration_B_phase_only

Phase only solution

Default: “False”

--calibration_B_timeslice

Solution length (s)

--calibration_global_solution

Solve across frequency

Default: “True”

--calibration_context

Terms to solve (in order e.g. TGB)

Default: “T”

--use_initial_skymodel

Whether to use an initial SkyModel in ICAL or not

Default: False

--input_skycomponent_file

Input name of skycomponents file (in hdf or txt format) for initial SkyModel in ICAL

--num_bright_sources

Number of brightest sources to select for initial SkyModel (if None, use all sources from input file)

--calibrate_with_dp3

Enables calibration using DP3 Gaincal step (https://dp3.readthedocs.io/en/latest/steps/GainCal.html)

Default: False

--input_dp3_skymodel

Path to a .skymodel file as expected by DP3

--clean_algorithm

Type of deconvolution algorithm (hogbom or msclean or mmclean)

Default: “mmclean”

--clean_use_radler

If enabled, RADLER is used for deconvolution

Default: “False”

--clean_beam

Clean beam: major axis, minor axis, position angle (deg)

--clean_scales

Scales for multiscale clean (pixels) e.g. [0, 6, 10]

Default: [0]

--clean_nmoment

Number of frequency moments in mmclean (1 is a constant, 2 is linear, etc.)

Default: 4

--clean_nmajor

Number of major cycles in cip or ical

Default: 5

--clean_niter

Number of minor cycles in CLEAN (i.e. clean iterations)

Default: 1000

--clean_psf_support

Half-width of psf used in cleaning (pixels)

Default: 256

--clean_gain

Clean loop gain

Default: 0.1

--clean_threshold

Clean stopping threshold (Jy/beam)

Default: 0.0001

--clean_component_threshold

Sources with absolute flux > this level (Jy) are fit or extracted using skycomponents

--clean_component_method

Method to convert sources in image to skycomponents: ‘fit’ in frequency or ‘extract’ actual values

Default: “fit”

--clean_fractional_threshold

Fractional stopping threshold for major cycle

Default: 0.3

--clean_facets

Number of overlapping facets in faceted clean (along each axis)

Default: 1

--clean_overlap

Overlap of facets in clean (pixels)

Default: 32

--clean_taper

Type of interpolation between facets in deconvolution (none or linear or tukey)

Default: “tukey”

--clean_restore_facets

Number of overlapping facets in restore step (along each axis)

Default: 1

--clean_restore_overlap

Overlap of facets in restore step (pixels)

Default: 32

--clean_restore_taper

Type of interpolation between facets in restore step (none or linear or tukey)

Default: “tukey”

--clean_restored_output

Type of restored image output: taylor, list, or integrated

Default: “list”

--use_dask

Use Dask processing? False means that graphs are executed as they are constructed.

Default: “True”

--dask_nthreads

Number of threads in each Dask worker (None means Dask will choose)

--dask_memory

Memory per Dask worker (GB), e.g. 5GB (None means Dask will choose)

--dask_memory_usage_file

File in which to track Dask memory use (using dask-memusage)

--dask-nodes

Node names for SSHCluster

--dask_nworkers

Number of workers (None means Dask will choose)

--dask_scheduler

Externally defined Dask scheduler e.g. 127.0.0.1:8786 or ssh for SSHCluster or existing for current scheduler

--dask_scheduler_file

Externally defined Dask scheduler file to setup dask cluster

--dask_tcp_timeout

Dask TCP timeout

--dask_connect_timeout

Dask connect timeout

--dask_malloc_trim_threshold

Threshold for trimming memory on release (0 is aggressive)

Default: 0