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:
If no input file is provided, a point source at the phase centre, with brightness of 1 Jy is used as the component.
- If either an HDF file or a TXT file is provided, the components are read from the file.
if
--num_bright_sources
is left asNone
, all of the components are used for the SkyModelif
--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