Coverage for / home / casatest / venv / lib / python3.12 / site-packages / casatasks / private / imagerhelpers / _gclean.py: 61%
554 statements
« prev ^ index » next coverage.py v7.13.0, created at 2025-12-19 19:37 +0000
1# vim: syntax=python
2########################################################################
3#
4# _gclean.py
5#
6# Copyright (C) 2021,2022,2023
7# Associated Universities, Inc. Washington DC, USA.
8#
9# This script is free software; you can redistribute it and/or modify it
10# under the terms of the GNU Library General Public License as published by
11# the Free Software Foundation; either version 2 of the License, or (at your
12# option) any later version.
13#
14# This library is distributed in the hope that it will be useful, but WITHOUT
15# ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
16# FITNESS FOR A PARTICULAR PURPOSE. See the GNU Library General Public
17# License for more details.
18#
19# You should have received a copy of the GNU Library General Public License
20# along with this library; if not, write to the Free Software Foundation,
21# Inc., 675 Massachusetts Ave, Cambridge, MA 02139, USA.
22#
23# Correspondence concerning AIPS++ should be adressed as follows:
24# Internet email: casa-feedback@nrao.edu.
25# Postal address: AIPS++ Project Office
26# National Radio Astronomy Observatory
27# 520 Edgemont Road
28# Charlottesville, VA 22903-2475 USA
29#
31from __future__ import absolute_import
32from casatools.typecheck import CasaValidator as _val_ctor
33_pc = _val_ctor( )
34from casatools.coercetype import coerce as _coerce
35from casatools.errors import create_error_string
36from casatasks.private.task_logging import start_log as _start_log
37from casatasks.private.task_logging import end_log as _end_log
38from casatasks.private.task_logging import except_log as _except_log
40import os
41import asyncio
42from functools import reduce
43import copy
44import numpy as np
45import shutil
46import time
47import subprocess
49import sys
53from casatasks.private import cleanhelper
54from casatasks.private.imagerhelpers.imager_return_dict import ImagingDict, ConvergenceResult, StopCodes
55from casatasks.private.imagerhelpers.input_parameters import ImagerParameters, sanitize_tclean_inputs
56from casatasks import deconvolve, tclean, imstat
57from casatasks import casalog
59###
60### import check versions
61###
62_GCV001 = True
63_GCV002 = True
64_GCV003 = True
65_GCV004 = True
67# from casatasks.private.imagerhelpers._gclean import gclean
68class gclean:
69 '''
70 tclean ---- Radio Interferometric Image Reconstruction
72 Form images from visibilities and reconstruct a sky model.
73 This task handles continuum images and spectral line cubes,
74 supports outlier fields, contains standard clean based algorithms
75 along with algorithms for multi-scale and wideband image
76 reconstruction, widefield imaging correcting for the w-term,
77 full primary-beam imaging and joint mosaic imaging (with
78 heterogeneous array support for ALMA).
80 --------- parameter descriptions ---------------------------------------------
82 vis Name(s) of input visibility file(s)
83 default: none;
84 example: vis='ngc5921.ms'
85 vis=['ngc5921a.ms','ngc5921b.ms']; multiple MSs
86 selectdata Enable data selection parameters.
87 field to image or mosaic. Use field id(s) or name(s).
88 ['go listobs' to obtain the list id's or names]
89 default: ''= all fields.
90 If field string is a non-negative integer, it is assumed to
91 be a field index, otherwise it is assumed to be a
92 field name.
93 field='0~2'; field ids 0,1,2.
94 field='0,4,5~7'; field ids 0,4,5,6,7.
95 field='3C286,3C295'; field names 3C286 and 3C295.
96 field = '3,4C\*'; field id 3, all names starting with 4C.
97 For multiple MS input, a list of field strings can be used:
98 field = ['0~2','0~4']; field ids 0-2 for the first MS and 0-4
99 for the second.
100 field = '0~2'; field ids 0-2 for all input MSs.
101 spw l window/channels.
102 NOTE: channels not selected here will contain all zeros if
103 selected by other subparameters.
104 default: ''=all spectral windows and channels.
105 spw='0~2,4'; spectral windows 0,1,2,4 (all channels).
106 spw='0:5~61'; spw 0, channels 5 to 61.
107 spw='<2'; spectral windows less than 2 (i.e. 0,1).
108 spw='0,10,3:3~45'; spw 0,10 all channels, and spw 3
109 channels 3 to 45.
110 spw='0~2:2~6'; spw 0,1,2 with channels 2 through 6 in each.
111 For multiple MS input, a list of spw strings can be used:
112 spw=['0','0~3']; spw ids 0 for the first MS and 0-3 for the second.
113 spw='0~3' spw ids 0-3 for all input MS.
114 spw='3:10~20;50~60' for multiple channel ranges within spw id 3.
115 spw='3:10~20;50~60,4:0~30' for different channel ranges for spw ids 3 and 4.
116 spw='0:0~10,1:20~30,2:1;2;3'; spw 0, channels 0-10,
117 spw 1, channels 20-30, and spw 2, channels, 1,2 and 3.
118 spw='1~4;6:15~48' for channels 15 through 48 for spw ids 1,2,3,4 and 6.
119 timerange Range of time to select from data
121 default: '' (all); examples,
122 timerange = 'YYYY/MM/DD/hh:mm:ss~YYYY/MM/DD/hh:mm:ss'
123 Note: if YYYY/MM/DD is missing date defaults to first
124 day in data set.
125 timerange='09:14:0~09:54:0' picks 40 min on first day.
126 timerange='25:00:00~27:30:00' picks 1 hr to 3 hr
127 30min on NEXT day.
128 timerange='09:44:00' pick data within one integration
129 of time.
130 timerange='> 10:24:00' data after this time.
131 For multiple MS input, a list of timerange strings can be
132 used:
133 timerange=['09:14:0~09:54:0','> 10:24:00'].
134 timerange='09:14:0~09:54:0''; apply the same timerange for
135 all input MSs.
136 uvrange Select data within uvrange (default unit is meters)
137 default: '' (all); example:
138 uvrange='0~1000klambda'; uvrange from 0-1000 kilo-lambda.
139 uvrange='> 4klambda';uvranges greater than 4 kilo lambda.
140 For multiple MS input, a list of uvrange strings can be
141 used:
142 uvrange=['0~1000klambda','100~1000klamda'].
143 uvrange='0~1000klambda'; apply 0-1000 kilo-lambda for all
144 input MSs.
145 uvrange='0~1000'; apply 0-1000 meter for all input MSs.
146 antenna Select data based on antenna/baseline
148 default: '' (all)
149 If antenna string is a non-negative integer, it is
150 assumed to be an antenna index, otherwise, it is
151 considered an antenna name.
152 antenna='5\&6'; baseline between antenna index 5 and
153 index 6.
154 antenna='VA05\&VA06'; baseline between VLA antenna 5
155 and 6.
156 antenna='5\&6;7\&8'; baselines 5-6 and 7-8.
157 antenna='5'; all baselines with antenna index 5.
158 antenna='05'; all baselines with antenna number 05
159 (VLA old name).
160 antenna='5,6,9'; all baselines with antennas 5,6,9
161 index number.
162 For multiple MS input, a list of antenna strings can be
163 used:
164 antenna=['5','5\&6'];
165 antenna='5'; antenna index 5 for all input MSs.
166 antenna='!DV14'; use all antennas except DV14.
167 scan Scan number range
169 default: '' (all).
170 example: scan='1~5'.
171 For multiple MS input, a list of scan strings can be used:
172 scan=['0~100','10~200'].
173 scan='0~100; scan ids 0-100 for all input MSs.
174 observation Observation ID range
175 default: '' (all).
176 example: observation='1~5'.
177 intent Scan Intent(s)
179 default: '' (all).
180 example: intent='TARGET_SOURCE'.
181 example: intent='TARGET_SOURCE1,TARGET_SOURCE2'.
182 example: intent='TARGET_POINTING\*'.
183 datacolumn Data column to image (data or observed, corrected)
184 default:'corrected'
185 ( If 'corrected' does not exist, it will use 'data' instead )
186 imagename Pre-name of output images
188 example : imagename='try'
190 Output images will be (a subset of) :
192 try.psf - Point Spread Function (PSF).
193 try.residual - Residual image.
194 try.image - Restored image.
195 try.model - Model image (contains only flux components).
196 try.sumwt - Single pixel image containing sum-of-weights.
197 (for natural weighting, sensitivity=1/sqrt(sumwt)).
198 try.pb - Primary Beam (PB) model (values depend on the gridder used).
200 A-projection algorithms (gridder=mosaic,awproject, awp2) will
201 compute the following images too.
203 try.weight - FT of gridded weights or the
204 un-normalized sum of PB-square (for all pointings).
205 Here, PB = sqrt(weight) normalized to a maximum of 1.0.
207 For multi-term wideband imaging, all relevant images above will
208 have additional .tt0,.tt1, etc suffixes to indicate Taylor terms,
209 plus the following extra output images.
210 try.alpha - spectral index.
211 try.alpha.error - estimate of error on spectral index.
212 try.beta - spectral curvature (if nterms \> 2).
214 Tip : Include a directory name in 'imagename' for all
215 output images to be sent there instead of the
216 current working directory : imagename='mydir/try'.
218 Tip : Restarting an imaging run without changing 'imagename'
219 implies continuation from the existing model image on disk.
220 - If 'startmodel' was initially specified it needs to be set to ""
221 for the restart run (or tclean will exit with an error message).
222 - By default, the residual image and psf will be recomputed
223 but if no changes were made to relevant parameters between
224 the runs, set calcres=False, calcpsf=False to resume directly from
225 the minor cycle without the (unnecessary) first major cycle.
226 To automatically change 'imagename' with a numerical
227 increment, set restart=False (see tclean docs for 'restart').
229 Note : All imaging runs will by default produce restored images.
230 For a niter=0 run, this will be redundant and can optionally
231 be turned off via the 'restoration=T/F' parameter.
232 imsize Number of pixels
233 example:
235 imsize = [350,250].
236 imsize = 500 is equivalent to [500,500].
238 To take proper advantage of internal optimized FFT routines, the
239 number of pixels must be even and factorizable by 2,3,5 only.
240 To find the nearest optimal imsize to that desired by the user, please use the following tool method:
242 from casatools import synthesisutils
243 su = synthesisutils()
244 su.getOptimumSize(345)
245 Output : 360
246 cell Cell size
247 example: cell=['0.5arcsec,'0.5arcsec'] or
248 cell=['1arcmin', '1arcmin'].
249 cell = '1arcsec' is equivalent to ['1arcsec','1arcsec'].
250 phasecenter Phase center of the image (string or field id); if the phasecenter is the name known major solar system object ('MERCURY', 'VENUS', 'MARS', 'JUPITER', 'SATURN', 'URANUS', 'NEPTUNE', 'PLUTO', 'SUN', 'MOON') or is an ephemerides table then that source is tracked and the background sources get smeared. There is a special case, when phasecenter='TRACKFIELD', which will use the ephemerides or polynomial phasecenter in the FIELD table of the MS's as the source center to track.
252 Note : If unspecified, tclean will use the phase-center from the first data field of the MS (or list of MSs) selected for imaging.
254 example: phasecenter='6'.
255 phasecenter='J2000 19h30m00 -40d00m00'.
256 phasecenter='J2000 292.5deg -40.0deg'.
257 phasecenter='J2000 5.105rad -0.698rad'.
258 phasecenter='ICRS 13:05:27.2780 -049.28.04.458'.
259 phasecenter='myComet_ephem.tab'.
260 phasecenter='MOON'.
261 phasecenter='TRACKFIELD'.
262 stokes Stokes Planes to make
263 default='I'; example: stokes='IQUV';
264 Options: 'I','Q','U','V','IV','QU','IQ','UV','IQUV','RR','LL','XX','YY','RRLL','XXYY','pseudoI'
266 Note : Due to current internal code constraints, if any correlation pair
267 is flagged, by default, no data for that row in the MS will be used.
268 So, in an MS with XX,YY, if only YY is flagged, neither a
269 Stokes I image nor an XX image can be made from those data points.
270 In such a situation, please split out only the unflagged correlation into
271 a separate MS, or use the option 'pseudoI'.
273 Note : The 'pseudoI' option is a partial solution, allowing Stokes I imaging
274 when either of the parallel-hand correlations are unflagged.
276 The remaining constraints shall be removed (where logical) in a future release.
277 projection Coordinate projection
278 Examples : SIN, NCP.
279 A list of supported (but untested) projections can be found here :
280 http://casa.nrao.edu/active/docs/doxygen/html/classcasa_1_1Projection.html#a3d5f9ec787e4eabdce57ab5edaf7c0cd
281 startmodel Name of starting model image
283 The contents of the supplied starting model image will be
284 copied to the imagename.model before the run begins.
286 example : startmodel = 'singledish.im'.
288 For deconvolver='mtmfs', one image per Taylor term must be provided.
289 example : startmodel = ['try.model.tt0', 'try.model.tt1'].
290 startmodel = ['try.model.tt0'] will use a starting model only
291 for the zeroth order term.
292 startmodel = ['','try.model.tt1'] will use a starting model only
293 for the first order term.
295 This starting model can be of a different image shape and size from
296 what is currently being imaged. If so, an image regrid is first triggered
297 to resample the input image onto the target coordinate system.
299 A common usage is to set this parameter equal to a single dish image.
301 Negative components in the model image will be included as is.
303 Note : If an error occurs during image resampling/regridding,
304 please try using task imregrid to resample the starting model
305 image onto a CASA image with the target shape and
306 coordinate system before supplying it via startmodel.
307 specmode Spectral definition mode (mfs,cube,cubedata, cubesource, mvc)
309 specmode='mfs' : Continuum imaging with only one output image channel.
310 (mode='cont' can also be used here)
312 specmode='cube' : Spectral line imaging with one or more channels.
313 Parameters start, width,and nchan define the spectral
314 coordinate system and can be specified either in terms
315 of channel numbers, frequency or velocity in whatever
316 spectral frame is specified in 'outframe'.
317 All internal and output images are made with outframe as the
318 base spectral frame. However imaging code internally uses the fixed
319 spectral frame, LSRK for automatic internal software
320 Doppler correction, so that a spectral line observed over an
321 extended time range will line up appropriately.
322 Therefore the output images have additional spectral frame conversion
323 layer in LSRK on the top the base frame.
328 specmode='cubedata' : Spectral line imaging with one or more channels.
329 There is no internal software Doppler correction, so
330 a spectral line observed over an extended time range
331 may be smeared out in frequency. There is strictly
332 no valid spectral frame with which to associate with the
333 output images, thus the image spectral frame will
334 be labelled "Undefined".
337 specmode='cubesource': Spectral line imaging while
338 tracking moving source (near field or solar system
339 objects). The velocity of the source is accounted
340 and the frequency reported is in the source frame.
341 As there is no "SOURCE" frame defined in CASA,
342 the frame in the image will be labelled "REST" (but do note the
343 velocity of a given line reported may be different from the rest frame
344 velocity if the emission region is moving w.r.t the systemic
345 velocity frame of the source).
347 specmode='mvc' : Multiterm continuum imaging with cube major cycles.
348 This mode requires deconvolver='mtmfs' with nterms>1
349 and user-set choices of 'reffreq' and 'nchan'.
351 The output images and minor cycle are similar to specmode='mfs'
352 with deconvolver='mtmfs', but the major cycles are done in
353 cube mode (and require a setting of 'reffreq' and 'nchan').
354 By default, frequency-dependent primary beam correction is
355 applied to each channel, before being combined across frequency
356 to make the inputs to the 'mtmfs' deconvolver. This results in
357 implicit wideband pb-correction, with the deconvolver seeing only
358 the sky spectral structure.
360 Note : There is currently no option to turn off wideband pb correction
361 as part of the flat-sky normalization between the major and minor cycles.
362 Therefore, 'mvc' with the 'standard' and 'wproject' gridders will also apply
363 pblimits per channel, masking all regions outside of pblimit.
364 An option to retain sources outside the pblimit will be added in a future release.
366 Note : Below is some guidance for choosing 'nchan' and 'reffreq' :
368 The cube produced by the major cycle is used in a linear least square fits for Taylor
369 polynomials per pixel. Therefore, one only needs as many channels in the cube, as
370 required for an accurate polynomial fit for sources that have the strongest
371 spectral structure.
373 In general, 'nchan' needs to be greater than or equal to 'nterms', and the
374 frequency range selected by the data will be evenly split into nchan channels.
375 For a low-order polynomial fit, only a small number (around 10)
376 channels are typically needed (for VLA/ALMA bandwidth ratios).
377 'nchan=-1' applies a heuristic that results in a default of 10 cube channels
378 for a 2:1 bandwidth ratio.
380 nchan = MAX( bandwidth/(0.1*startfreq) , nterms+1 )
382 Note: When running in parallel, the nchan selected may limit the speedup if it
383 is smaller than the number of processes used.
385 The 'reffreq' is the reference frequency used for the Taylor polynomial expansion.
386 By default, in specmode='mvc', reffreq is set to the middle of the selected
387 frequency range.
388 reffreq Reference frequency of the output image coordinate system.
390 Example : reffreq='1.5GHz' as a string with units.
392 By default, it is calculated as the middle of the selected frequency range.
394 For deconvolver='mtmfs' the Taylor expansion is also done about
395 this specified reference frequency.
396 nchan Number of channels in the output image.
397 For default (=-1), the number of channels will be automatically determined
398 based on data selected by 'spw' with 'start' and 'width'.
399 It is often easiest to leave nchan at the default value.
400 example: nchan=100
401 start First channel (e.g. start=3,start=\'1.1GHz\',start=\'15343km/s\')
402 of output cube images specified by data channel number (integer),
403 velocity (string with a unit), or frequency (string with a unit).
404 Default:''; The first channel is automatically determined based on
405 the 'spw' channel selection and 'width'.
406 channels in 'spw'.
407 Since the integer number in 'start' represents the data channel number,
408 when the channel number is used along with the spectral window id selection
409 in 'spw', 'start' specified as an integer should be carefully set otherwise
410 it may result in the blank image channels if the 'start' channel (i.e. absolute
411 channel number) is outside of the channel range specified in 'spw'.
412 In such a case, 'start' can be left as a default (='') to ensure
413 matching with the data spectral channel selection.
414 For specmode='cube', when velocity or frequency is used it is
415 interpreted with the frame defined in outframe. [The parameters of
416 the desired output cube can be estimated by using the 'transform'
417 functionality of 'plotms'].
418 examples: start='5.0km/s'; 1st channel, 5.0km/s in outframe.
419 start='22.3GHz'; 1st channel, 22.3GHz in outframe.
420 width Channel width (e.g. width=2,width=\'0.1MHz\',width=\'10km/s\') of output cube images
421 specified by data channel number (integer), velocity (string with a unit), or
422 or frequency (string with a unit).
423 Default:''; data channel width.
424 The sign of width defines the direction of the channels to be incremented.
425 For width specified in velocity or frequency with '-' in front gives image channels in
426 decreasing velocity or frequency, respectively.
427 For specmode='cube', when velocity or frequency is used it is interpreted with
428 the reference frame defined in outframe.
429 examples: width='2.0km/s'; results in channels with increasing velocity.
430 width='-2.0km/s'; results in channels with decreasing velocity.
431 width='40kHz'; results in channels with increasing frequency.
432 width=-2; results in channels averaged of 2 data channels incremented from
433 high to low channel numbers.
434 outframe Spectral reference frame in which to interpret \'start\' and \'width\'
435 Options: '','LSRK','LSRD','BARY','GEO','TOPO','GALACTO','LGROUP','CMB'
436 example: outframe='bary' for Barycentric frame.
438 REST -- Rest frequency.
439 LSRD -- Local Standard of Rest (J2000).
440 -- as the dynamical definition (IAU, [9,12,7] km/s in galactic coordinates).
441 LSRK -- LSR as a kinematical (radio) definition.
442 -- 20.0 km/s in direction ra,dec = [270,+30] deg (B1900.0).
443 BARY -- Barycentric (J2000).
444 GEO --- Geocentric.
445 TOPO -- Topocentric.
446 GALACTO -- Galacto centric (with rotation of 220 km/s in direction l,b = [90,0] deg.
447 LGROUP -- Local group velocity -- 308km/s towards l,b = [105,-7] deg (F. Ghigo).
448 CMB -- CMB velocity -- 369.5km/s towards l,b = [264.4, 48.4] deg (F. Ghigo).
449 DEFAULT = LSRK.
450 veltype Velocity type (radio, z, ratio, beta, gamma, optical)
451 For 'start' and/or 'width' specified in velocity, specifies the velocity definition
452 Options: 'radio','optical','z','beta','gamma','optical'
453 NOTE: the viewer always defaults to displaying the 'radio' frame,
454 but that can be changed in the position tracking pull down.
456 The different types (with F = f/f0, the frequency ratio), are:
458 Z = (-1 + 1/F).
459 RATIO = (F) \*.
460 RADIO = (1 - F).
461 OPTICAL == Z.
462 BETA = ((1 - F^2)/(1 + F^2)).
463 GAMMA = ((1 + F^2)/2F) \*.
464 RELATIVISTIC == BETA (== v/c).
465 DEFAULT == RADIO.
466 Note that the ones with an '\*' have no real interpretation
467 (although the calculation will proceed) if given as a velocity.
468 restfreq List of rest frequencies or a rest frequency in a string.
469 Specify rest frequency to use for output image.
471 Currently it uses the first rest frequency in the list for translation of
472 velocities. The list will be stored in the output images.
473 Default: []; look for the rest frequency stored in the MS, if not available,
474 use center frequency of the selected channels.
475 examples: restfreq=['1.42GHz'].
476 restfreq='1.42GHz'.
477 interpolation Spectral interpolation (nearest,linear,cubic)
479 Interpolation rules to use when binning data channels onto image channels
480 and evaluating visibility values at the centers of image channels.
482 Note : 'linear' and 'cubic' interpolation requires data points on both sides of
483 each image frequency. Errors are therefore possible at edge channels, or near
484 flagged data channels. When image channel width is much larger than the data
485 channel width there is nothing much to be gained using linear or cubic thus
486 not worth the extra computation involved.
487 perchanweightdensity When calculating weight density for Briggs
488 style weighting in a cube, this parameter
489 determines whether to calculate the weight
490 density for each channel independently
491 (the default, True)
492 or a common weight density for all of the selected
493 data. This parameter has no
494 meaning for continuum (specmode='mfs') imaging
495 or for natural and radial weighting schemes.
496 For cube imaging
497 perchanweightdensity=True is a recommended
498 option that provides more uniform
499 sensitivity per channel for cubes, but with
500 generally larger psfs than the
501 perchanweightdensity=False (prior behavior)
502 option. When using Briggs style weight with
503 perchanweightdensity=True, the imaging weight
504 density calculations use only the weights of
505 data that contribute specifically to that
506 channel. On the other hand, when
507 perchanweightdensity=False, the imaging
508 weight density calculations sum all of the
509 weights from all of the data channels
510 selected whose (u,v) falls in a given uv cell
511 on the weight density grid. Since the
512 aggregated weights, in any given uv cell,
513 will change depending on the number of
514 channels included when imaging, the psf
515 calculated for a given frequency channel will
516 also necessarily change, resulting in
517 variability in the psf for a given frequency
518 channel when perchanweightdensity=False. In
519 general, perchanweightdensity=False results
520 in smaller psfs for the same value of
521 robustness compared to
522 perchanweightdensity=True, but the rms noise
523 as a function of channel varies and increases
524 toward the edge channels;
525 perchanweightdensity=True provides more
526 uniform sensitivity per channel for
527 cubes. This may make it harder to find
528 estimates of continuum when
529 perchanweightdensity=False. If you intend to
530 image a large cube in many smaller subcubes
531 and subsequently concatenate, it is advisable
532 to use perchanweightdensity=True to avoid
533 surprisingly varying sensitivity and psfs
534 across the concatenated cube.
535 gridder Gridding options (standard, wproject, widefield, mosaic, awproject, awp2, awphpg)
538 The following options choose different gridding convolution
539 functions for the process of convolutional resampling of the measured
540 visibilities onto a regular uv-grid prior to an inverse FFT.
541 Model prediction (degridding) also uses these same functions.
542 Several wide-field effects can be accounted for via careful choices of
543 convolution functions. Gridding (degridding) runtime will rise in
544 proportion to the support size of these convolution functions (in uv-pixels).
546 standard : Prolate Spheroid with 7x7 uv pixel support size.
548 [ This mode can also be invoked using 'ft' or 'gridft' ]
550 wproject : W-Projection algorithm to correct for the widefield
551 non-coplanar baseline effect. [Cornwell et.al 2008]
553 wprojplanes is the number of distinct w-values at
554 which to compute and use different gridding convolution
555 functions (see help for wprojplanes).
556 Convolution function support size can range
557 from 5x5 to few 100 x few 100.
559 [ This mode can also be invoked using 'wprojectft' ]
561 widefield : Facetted imaging with or without W-Projection per facet.
563 A set of facets x facets subregions of the specified image
564 are gridded separately using their respective phase centers
565 (to minimize max W). Deconvolution is done on the joint
566 full size image, using a PSF from the first subregion.
568 wprojplanes=1 : standard prolate spheroid gridder per facet.
569 wprojplanes > 1 : W-Projection gridder per facet.
570 nfacets=1, wprojplanes > 1 : Pure W-Projection and no facetting.
571 nfacets=1, wprojplanes=1 : Same as standard,ft,gridft.
573 A combination of facetting and W-Projection is relevant only for
574 very large fields of view. (In our current version of tclean, this
575 combination runs only with parallel=False.
577 mosaic : A-Projection with azimuthally symmetric beams without
578 sidelobes, beam rotation or squint correction.
579 Gridding convolution functions per visibility are computed
580 from FTs of PB models per antenna.
581 This gridder can be run on single fields as well as mosaics.
583 VLA : PB polynomial fit model (Napier and Rots, 1982).
584 EVLA : PB polynomial fit model (Perley, 2015).
585 ALMA : Airy disks for a 10.7m dish (for 12m dishes) and
586 6.25m dish (for 7m dishes) each with 0.75m
587 blockages (Hunter/Brogan 2011). Joint mosaic
588 imaging supports heterogeneous arrays for ALMA.
590 Typical gridding convolution function support sizes are
591 between 7 and 50 depending on the desired
592 accuracy (given by the uv cell size or image field of view).
594 [ This mode can also be invoked using 'mosaicft' or 'ftmosaic' ]
596 awproject : A-Projection with azimuthally asymmetric beams and
597 including beam rotation, squint correction,
598 conjugate frequency beams and W-projection.
599 [Bhatnagar et.al, 2008]
601 Gridding convolution functions are computed from
602 aperture illumination models per antenna and optionally
603 combined with W-Projection kernels and a prolate spheroid.
604 This gridder can be run on single fields as well as mosaics.
606 VLA : Uses ray traced model (VLA and EVLA) including feed
607 leg and subreflector shadows, off-axis feed location
608 (for beam squint and other polarization effects), and
609 a Gaussian fit for the feed beams (Brisken 2009)
610 ALMA : Similar ray-traced model as above (but the correctness
611 of its polarization properties remains un-verified).
613 Typical gridding convolution function support sizes are
614 between 7 and 50 depending on the desired
615 accuracy (given by the uv cell size or image field of view).
616 When combined with W-Projection they can be significantly larger.
618 [ This mode can also be invoked using 'awprojectft' ]
621 awp2 : A-Projection with azimuthally asymmetric beams and
622 including beam rotation, squint correction and W-projection.
623 [Bhatnagar et.al, 2008]
625 Gridding convolution functions are computed from
626 aperture illumination models (assuming similar antennas) and optionally
627 combined with W-Projection kernels.
628 This gridder can be run on single fields as well as mosaics.
629 The other sub-parameters that are of significance when using this gridder
630 are wprojplanes, computepastep, mosweight, usepointing, pblimit and normtype.
632 Only supports VLA : Uses ray traced model (VLA and EVLA) including feed
633 leg and subreflector shadows, off-axis feed location
634 (for beam squint and other polarization effects), and
635 a Gaussian fit for the feed beams (Ref: Brisken 2009)
637 For squint correction the value passed in computepastep has to be smaller than 180.
638 Anything larger awp2 will use an average of LL and RR beams. If computepastep=5, for
639 e.g., PB every 5 degrees, over the range of parallactic angle covered by the data, will be calculated
640 and the nearest beam to every integration will be used to correct for the squint between the L and R beams.
642 NOTE : For mtmfs with nterms >1 and using awp2 gridder, for accurate results always use specmode="mvc"
643 as awp2 with specmode="mfs" does not use conjugate beams to remove the spectral
644 index of the primary beam.
646 awphpg : Implementation of the high performance gridder (HPG; Pokorny, ngVLA Computing Memo #5).
647 For CASA 6.7.0 this mode is only available on the internal VLASS release of CASA.
648 It will be made available for general use in a future CASA release.
651 imagemosaic : (untested implementation).
653 Grid and iFT each pointing separately and combine the
654 images as a linear mosaic (weighted by a PB model) in
655 the image domain before a joint minor cycle.
657 VLA/ALMA PB models are same as for gridder='mosaicft'.
659 ------ Notes on PB models :
661 (1) Several different sources of PB models are used in the modes
662 listed above. This is partly for reasons of algorithmic flexibility
663 and partly due to the current lack of a common beam model
664 repository or consensus on what beam models are most appropriate.
666 (2) For ALMA and gridder='mosaic', ray-traced (TICRA) beams
667 are also available via the vpmanager tool.
668 For example, call the following before the tclean run.
669 vp.setpbimage(telescope="ALMA",
670 compleximage='/home/casa/data/trunk/alma/responses/ALMA_0_DV__0_0_360_0_45_90_348.5_373_373_GHz_ticra2007_VP.im',
671 antnames=['DV'+'%02d'%k for k in range(25)])
672 vp.saveastable('mypb.tab')
673 Then, supply vptable='mypb.tab' to tclean.
674 ( Currently this will work only for non-parallel runs )
677 ------ Note on PB masks :
679 In tclean, A-Projection gridders (mosaic, awproject, and awp2) produce a
680 .pb image and use the 'pblimit' subparameter to decide normalization
681 cutoffs and construct an internal T/F mask in the .pb and .image images.
682 However, this T/F mask cannot directly be used during deconvolution
683 (which needs a 1/0 mask). There are two options for making a pb based
684 deconvolution mask.
685 -- Run tclean with niter=0 to produce the .pb, construct a 1/0 image
686 with the desired threshold (using ia.open('newmask.im');
687 ia.calc('iif("xxx.pb">0.3,1.0,0.0)');ia.close() for example),
688 and supply it via the 'mask' parameter in a subsequent run
689 (with calcres=F and calcpsf=F to restart directly from the minor cycle).
690 -- Run tclean with usemask='pb' for it to automatically construct
691 a 1/0 mask from the internal T/F mask from .pb at a fixed 0.2 threshold.
694 ----- Making PBs for gridders other than mosaic, awproject, awp2
696 After the PSF generation, a PB is constructed using the same
697 models used in gridder='mosaic' but just evaluated in the image
698 domain without consideration to weights.
699 facets Number of facets on a side
701 A set of (facets x facets) subregions of the specified image
702 are gridded separately using their respective phase centers
703 (to minimize max W). Deconvolution is done on the joint
704 full size image, using a PSF from the first subregion/facet.
706 In our current version of tclean, facets>1 may be used only
707 with parallel=False.
708 psfphasecenter For mosaic use psf centered on this
709 optional direction. You may need to use
710 this if for example the mosaic does not
711 have any pointing in the center of the
712 image. Another reason; as the psf is
713 approximate for a mosaic, this may help
714 to deconvolve a non central bright source
715 well and quickly.
717 example:
719 psfphasecenter='6' #center psf on field 6.
720 psfphasecenter='J2000 19h30m00 -40d00m00'.
721 psfphasecenter='J2000 292.5deg -40.0deg'.
722 psfphasecenter='J2000 5.105rad -0.698rad'.
723 psfphasecenter='ICRS 13:05:27.2780 -049.28.04.458'.
724 wprojplanes Number of distinct w-values at which to compute and use different
725 gridding convolution functions for W-Projection
727 An appropriate value of wprojplanes depends on the presence/absence
728 of a bright source far from the phase center, the desired dynamic
729 range of an image in the presence of a bright far out source,
730 the maximum w-value in the measurements, and the desired trade off
731 between accuracy and computing cost.
733 As a (rough) guide, VLA L-Band D-config may require a
734 value of 128 for a source 30arcmin away from the phase
735 center. A-config may require 1024 or more. To converge to an
736 appropriate value, try starting with 128 and then increasing
737 it if artifacts persist. W-term artifacts (for the VLA) typically look
738 like arc-shaped smears in a synthesis image or a shift in source
739 position between images made at different times. These artifacts
740 are more pronounced the further the source is from the phase center.
742 There is no harm in simply always choosing a large value (say, 1024)
743 but there will be a significant performance cost to doing so, especially
744 for gridder='awproject' where it is combined with A-Projection.
746 wprojplanes=-1 is an option for gridder='widefield' or 'wproject'
747 in which the number of planes is automatically computed.
748 vptable vpmanager
750 vptable="" : Choose default beams for different telescopes.
751 ALMA : Airy disks.
752 EVLA : old VLA models.
754 Other primary beam models can be chosen via the vpmanager tool.
756 Step 1 : Set up the vpmanager tool and save its state in a table.
758 vp.setpbpoly(telescope='EVLA', coeff=[1.0, -1.529e-3, 8.69e-7, -1.88e-10])
759 vp.saveastable('myvp.tab')
761 Step 2 : Supply the name of that table in tclean.
763 tclean(....., vptable='myvp.tab',....)
765 Please see the documentation for the vpmanager for more details on how to
766 choose different beam models. Work is in progress to update the defaults
767 for EVLA and ALMA.
769 Note : AWProjection currently does not use this mechanism to choose
770 beam models. It instead uses ray-traced beams computed from
771 parameterized aperture illumination functions, which are not
772 available via the vpmanager. So, gridder='awproject' does not allow
773 the user to set this parameter.
774 mosweight When doing Brigg's style weighting (including uniform) to perform the weight density calculation for each field indepedently if True. If False the weight density is calculated from the average uv distribution of all the fields.
775 aterm Use aperture illumination functions during gridding.
777 This parameter turns on the A-term of the AW-Projection gridder.
778 Gridding convolution functions are constructed from aperture illumination
779 function models of each antenna.
780 psterm Include the Prolate Spheroidal (PS) funtion as the anti-aliasing
781 operator in the gridding convolution functions used for gridding.
783 Setting this parameter to true is necessary when aterm is set to
784 false. It can be set to false when aterm is set to true, though
785 with this setting effects of aliasing may be there in the image,
786 particularly near the edges.
788 When set to true, the .pb images will contain the fourier transform
789 of the of the PS funtion.
791 For more information on the functional
792 effects of the psterm, aterm and wprojplanes settings, see the
793 'Wide-field Imaging' pages in CASA Docs (https://casadocs.readthedocs.io).
794 wbawp Use frequency dependent A-terms.
795 Scale aperture illumination functions appropriately with frequency
796 when gridding and combining data from multiple channels.
797 conjbeams Use conjugate frequency for wideband A-terms.
799 While gridding data from one frequency channel, choose a convolution
800 function from a 'conjugate' frequency such that the resulting baseline
801 primary beam is approximately constant across frequency. For a system in
802 which the primary beam scales with frequency, this step will eliminate
803 instrumental spectral structure from the measured data and leave only the
804 sky spectrum for the minor cycle to model and reconstruct [Bhatnagar et al., ApJ, 2013].
806 As a rough guideline for when this is relevant, a source at the half power
807 point of the PB at the center frequency will see an artificial spectral
808 index of -1.4 due to the frequency dependence of the PB [Sault and Wieringa, 1994].
809 If left uncorrected during gridding, this spectral structure must be modeled
810 in the minor cycle (using the mtmfs algorithm) to avoid dynamic range limits
811 (of a few hundred for a 2:1 bandwidth).
812 This works for specmode='mfs' and its value is ignored for cubes.
813 cfcache Convolution function cache directory name.
815 Name of a directory in which to store gridding convolution functions.
816 This cache is filled at the beginning of an imaging run. This step can be time
817 consuming but the cache can be reused across multiple imaging runs that
818 use the same image parameters (cell size, image size , spectral data
819 selections, wprojplanes, wbawp, psterm, aterm). The effect of the wbawp,
820 psterm and aterm settings is frozen-in in the cfcache. Using an existing cfcache
821 made with a different setting of these parameters will not reflect the current
822 settings.
824 In a parallel execution, the construction of the cfcache is also parallelized
825 and the time to compute scales close to linearly with the number of compute
826 cores used. With the re-computation of Convolution Functions (CF) due to PA
827 rotation turned-off (the computepastep parameter), the total number of in the
828 cfcache can be computed as [No. of wprojplanes x No. of selected spectral windows x 4]
830 By default, cfcache = imagename + '.cf'
831 usepointing The usepointing flag informs the gridder that it should utilize the pointing table
832 to use the correct direction in which the antenna is pointing with respect to the pointing phasecenter.
833 computepastep Parallactic angle interval after the AIFs are recomputed (deg).
835 This parameter controls the accuracy of the aperture illumination function
836 used with AProjection for alt-az mount dishes where the AIF rotates on the
837 sky as the synthesis image is built up. Once the PA in the data changes by
838 the given interval, AIFs are re-computed at the new PA.
840 A value of 360.0 deg (the default) implies no re-computation due to PA rotation.
841 AIFs are computed for the PA value of the first valid data received and used for
842 all of the data.
844 For gridder=awp2 a value of 180.0 deg or larger implies no squint correction will be
845 attempted i.e an average beam of the left hand and right hand polarization will be calculated
846 rotatepastep Parallactic angle interval after which the nearest AIF is rotated (deg)
848 Instead of recomputing the AIF for every timestep's parallactic angle,
849 the nearest existing AIF is used and rotated
850 after the PA changed by rotatepastep value.
852 A value of 360.0 deg (the default) disables rotation of the AIF.
854 For example, computepastep=360.0 and rotatepastep=5.0 will compute
855 the AIFs at only the starting parallactic angle and all other timesteps will
856 use a rotated version of that AIF at the nearest 5.0 degree point.
857 pointingoffsetsigdev Corrections for heterogenous and time-dependent pointing
858 offsets via AWProjection are controlled by this parameter.
859 It is a vector of 2 ints or doubles each of which is interpreted
860 in units of arcsec. Based on the first threshold, a clustering
861 algorithm is applied to entries from the POINTING subtable
862 of the MS to determine how distinct antenna groups for which
863 the pointing offset must be computed separately. The second
864 number controls how much a pointing change across time can
865 be ignored and after which an antenna rebinning is required.
868 Note : The default value of this parameter is [], due a programmatic constraint.
869 If run with this value, it will internally pick [600,600] and exercise the
870 option of using large tolerances (10arcmin) on both axes. Please choose
871 a setting explicitly for runs that need to use this parameter.
873 Note : This option is available only for gridder='awproject' and usepointing=True and
874 and has been validated primarily with VLASS on-the-fly mosaic data
875 where POINTING subtables have been modified after the data are recorded.
878 Examples of parameter usage :
880 [100.0,100.0] : Pointing offsets of 100 arcsec or less are considered
881 small enough to be ignored. Using large values for both
882 indicates a homogeneous array.
885 [10.0, 100.0] : Based on entries in the POINTING subtable, antennas
886 are grouped into clusters based on a 10arcsec bin size.
887 All antennas in a bin are given a pointing offset calculated
888 as the average of the offsets of all antennas in the bin.
889 On the time axis, offset changes upto 100 arcsec will be ignored.
891 [10.0,10.0] : Calculate separate pointing offsets for each antenna group
892 (with a 10 arcsec bin size). As a function of time, recalculate
893 the antenna binning if the POINTING table entries change by
894 more than 10 arcsec w.r.to the previously computed binning.
896 [1.0, 1.0] : Tight tolerances will imply a fully heterogenous situation where
897 each antenna gets its own pointing offset. Also, time-dependent
898 offset changes greater than 1 arcsec will trigger recomputes of
899 the phase gradients. This is the most general situation and is also
900 the most expensive option as it constructs and uses separate
901 phase gradients for all baselines and timesteps.
903 For VLASS 1.1 data with two kinds of pointing offsets, the recommended
904 setting is [ 30.0, 30.0 ].
906 For VLASS 1.2 data with only the time-dependent pointing offsets, the
907 recommended setting is [ 300.0, 30.0 ] to turn off the antenna grouping
908 but to retain the time dependent corrections required from one timestep
909 to the next.
910 pblimit PB gain level at which to cut off normalizations.
912 Divisions by .pb during normalizations have a cut off at a .pb gain
913 level given by pblimit. Outside this limit, image values are set to zero.
914 Additionally, by default, an internal T/F mask is applied to the .pb, .image and
915 .residual images to mask out (T) all invalid pixels outside the pblimit area.
917 Note : This internal T/F mask cannot be used as a deconvolution mask.
918 To do so, please follow the steps listed above in the Notes for the
919 'gridder' parameter.
921 Note : To prevent the internal T/F mask from appearing in anything other
922 than the .pb and .image.pbcor images, 'pblimit' can be set to a
923 negative number.
924 The absolute value will still be used as a valid 'pblimit' for normalization
925 purposes. So, for example, pick pblimit=-0.1 (and not pblimit=-1).
926 A tclean restart using existing output images on disk that already
927 have this T/F mask in the .residual and .image but only pblimit set
928 to a negative value, will remove this mask after the next major cycle.
930 Note : An existing internal T/F mask may be removed from an image as
931 follows (without needing to re-run tclean itself).
932 ia.open('test.image');
933 ia.maskhandler(op='set', name='');
934 ia.done()
935 normtype Normalization type (flatnoise, flatsky, pbsquare).
937 Gridded (and FT'd) images represent the PB-weighted sky image.
938 Qualitatively it can be approximated as two instances of the PB
939 applied to the sky image (one naturally present in the data
940 and one introduced during gridding via the convolution functions).
942 xxx.weight : Weight image approximately equal to sum ( square ( pb ) )
943 xxx.pb : Primary beam calculated as sqrt ( xxx.weight )
945 normtype='flatnoise' : Divide the raw image by sqrt(.weight) so that
946 the input to the minor cycle represents the
947 product of the sky and PB. The noise is 'flat'
948 across the region covered by each PB.
950 normtype='flatsky' : Divide the raw image by .weight so that the input
951 to the minor cycle represents only the sky.
952 The noise is higher in the outer regions of the
953 primary beam where the sensitivity is low.
955 normtype='pbsquare' : No normalization after gridding and FFT.
956 The minor cycle sees the sky times pb square
957 deconvolver Name of minor cycle algorithm (hogbom,clark,multiscale,mtmfs,mem,clarkstokes,asp)
959 Each of the following algorithms operate on residual images and PSFs
960 from the gridder and produce output model and restored images.
961 Minor cycles stop and a major cycle is triggered when cyclethreshold
962 or cycleniter are reached. For all methods, components are picked from
963 the entire extent of the image or (if specified) within a mask.
965 hogbom : An adapted version of Hogbom Clean [Hogbom, 1974].
966 - Find the location of the peak residual.
967 - Add this delta function component to the model image.
968 - Subtract a scaled and shifted PSF of the same size as the image
969 from regions of the residual image where the two overlap.
970 - Repeat.
972 clark : An adapted version of Clark Clean [Clark, 1980].
973 - Find the location of max(I^2+Q^2+U^2+V^2).
974 - Add delta functions to each stokes plane of the model image.
975 - Subtract a scaled and shifted PSF within a small patch size
976 from regions of the residual image where the two overlap.
977 - After several iterations trigger a Clark major cycle to subtract
978 components from the visibility domain, but without de-gridding.
979 - Repeat.
981 ( Note : 'clark' maps to imagermode='' in the old clean task.
982 'clark_exp' is another implementation that maps to
983 imagermode='mosaic' or 'csclean' in the old clean task
984 but the behavior is not identical. For now, please
985 use deconvolver='hogbom' if you encounter problems. )
987 clarkstokes : Clark Clean operating separately per Stokes plane.
989 (Note : 'clarkstokes_exp' is an alternate version. See above.)
991 multiscale : MultiScale Clean [Cornwell, 2008].
992 - Smooth the residual image to multiple scale sizes.
993 - Find the location and scale at which the peak occurs.
994 - Add this multiscale component to the model image.
995 - Subtract a scaled,smoothed,shifted PSF (within a small
996 patch size per scale) from all residual images.
997 - Repeat from step 2.
999 mtmfs : Multi-term (Multi Scale) Multi-Frequency Synthesis [Rau and Cornwell, 2011].
1000 - Smooth each Taylor residual image to multiple scale sizes.
1001 - Solve a NTxNT system of equations per scale size to compute
1002 Taylor coefficients for components at all locations.
1003 - Compute gradient chi-square and pick the Taylor coefficients
1004 and scale size at the location with maximum reduction in
1005 chi-square.
1006 - Add multi-scale components to each Taylor-coefficient
1007 model image.
1008 - Subtract scaled,smoothed,shifted PSF (within a small patch size
1009 per scale) from all smoothed Taylor residual images.
1010 - Repeat from step 2.
1013 mem : Maximum Entropy Method [Cornwell and Evans, 1985].
1014 - Iteratively solve for values at all individual pixels via the
1015 MEM method. It minimizes an objective function of
1016 chi-square plus entropy (here, a measure of difference
1017 between the current model and a flat prior model).
1019 (Note : This MEM implementation is not very robust.
1020 Improvements will be made in the future.)
1022 asp : Adaptive Scale Pixel algorithm [Bhatnagar and Cornwell, 2004].
1023 - Define a set of initial scales defined as 0, W, 2W 4W and 8W.
1024 where W is a 2D Gaussian fitting width to the PSF.
1025 - Smooth the residual image by a Gaussian beam at initial scales.
1026 - Search for the global peak (F) among these smoothed residual images.
1027 - form an active Aspen set: amplitude(F), amplitude location(x,y).
1028 - Optimize the Aspen set by minimizing the objective function RI-Aspen*PSF,
1029 where RI is the residual image and * is the convulition operation.
1030 - Compute the model image and update the residual image
1031 - Repeat from step 2
1032 scales List of scale sizes (in pixels) for multi-scale and mtmfs algorithms.
1033 --> scales=[0,6,20]
1034 This set of scale sizes should represent the sizes
1035 (diameters in units of number of pixels)
1036 of dominant features in the image being reconstructed.
1038 The smallest scale size is recommended to be 0 (point source),
1039 the second being the size of the synthesized beam and the third being 3-5
1040 times the synthesized beam, etc. For example, if the synthesized
1041 beam is 10" FWHM and cell=2",try scales = [0,5,15].
1043 For numerical stability, the largest scale must be
1044 smaller than the image (or mask) size and smaller than or
1045 comparable to the scale corresponding to the lowest measured
1046 spatial frequency (as a scale size much larger than what the
1047 instrument is sensitive to is unconstrained by the data making
1048 it harder to recover from errors during the minor cycle).
1049 nterms Number of Taylor coefficients in the spectral model.
1051 - nterms=1 : Assume flat spectrum source.
1052 - nterms=2 : Spectrum is a straight line with a slope.
1053 - nterms=N : A polynomial of order N-1.
1055 From a Taylor expansion of the expression of a power law, the
1056 spectral index is derived as alpha = taylorcoeff_1 / taylorcoeff_0.
1058 Spectral curvature is similarly derived when possible.
1060 The optimal number of Taylor terms depends on the available
1061 signal to noise ratio, bandwidth ratio, and spectral shape of the
1062 source as seen by the telescope (sky spectrum x PB spectrum).
1064 nterms=2 is a good starting point for wideband EVLA imaging
1065 and the lower frequency bands of ALMA (when fractional bandwidth
1066 is greater than 10%) and if there is at least one bright source for
1067 which a dynamic range of greater than few 100 is desired.
1069 Spectral artifacts for the VLA often look like spokes radiating out from
1070 a bright source (i.e. in the image made with standard mfs imaging).
1071 If increasing the number of terms does not eliminate these artifacts,
1072 check the data for inadequate bandpass calibration. If the source is away
1073 from the pointing center, consider including wide-field corrections too.
1075 (Note : In addition to output Taylor coefficient images .tt0,.tt1,etc
1076 images of spectral index (.alpha), an estimate of error on
1077 spectral index (.alpha.error) and spectral curvature (.beta,
1078 if nterms is greater than 2) are produced.
1079 - These alpha, alpha.error and beta images contain
1080 internal T/F masks based on a threshold computed
1081 as peakresidual/10. Additional masking based on
1082 .alpha/.alpha.error may be desirable.
1083 - .alpha.error is a purely empirical estimate derived
1084 from the propagation of error during the division of
1085 two noisy numbers (alpha = xx.tt1/xx.tt0) where the
1086 'error' on tt1 and tt0 are simply the values picked from
1087 the corresponding residual images. The absolute value
1088 of the error is not always accurate and it is best to interpret
1089 the errors across the image only in a relative sense.
1090 smallscalebias A numerical control to bias the scales when using multi-scale or mtmfs algorithms.
1091 The peak from each scale's smoothed residual is
1092 multiplied by ( 1 - smallscalebias \* scale/maxscale )
1093 to increase or decrease the amplitude relative to other scales,
1094 before the scale with the largest peak is chosen.
1095 Smallscalebias can be varied between -1.0 and 1.0.
1096 A score of 0.0 gives all scales equal weight (default).
1097 A score larger than 0.0 will bias the solution towards smaller scales.
1098 A score smaller than 0.0 will bias the solution towards larger scales.
1099 The effect of smallscalebias is more pronounced when using multi-scale relative to mtmfs.
1100 fusedthreshold ring Hogbom Clean (number in units of Jy).
1102 fusedthreshold = 0.0001 : 0.1 mJy.
1104 This is a subparameter of the Asp Clean deconvolver. When peak residual
1105 is lower than the threshold, Asp Clean is "switched to Hogbom Clean" (i.e. only use the 0 scale for cleaning) for
1106 the following number of iterations until it switches back to Asp Clean.
1108 NumberIterationsInHogbom = 50 + 2 * (exp(0.05 * NthHogbom) - 1)
1110 , where NthHogbom is the number of times Hogbom Clean has been triggered.
1112 When the Asp Clean detects it is approaching convergence, it uses only the 0 scale for the following number of iterations for better computational efficiency.
1114 NumberIterationsInHogbom = 500 + 2 * (exp(0.05 * NthHogbom) - 1)
1116 Set 'fusedthreshold = -1' to make the Asp Clean deconvolver never "switch" to Hogbom Clean.
1117 largestscale xels) allowed for the initial guess for the Asp Clean deconvolver.
1119 largestscale = 100
1121 The default initial scale sizes used by Asp Clean is [0, w, 2w, 4w, 8w],
1122 where `w` is the PSF width. The default `largestscale` is -1 which indicates
1123 users accept these initial scales. If `largestscale` is set, the initial scales
1124 would be [0, w, ... up to the `largestscale`]. This is only an initial guess,
1125 and actual fitted scale sizes may evolve from these initial values.
1127 It is recommended not to set `largestscale` unless Asp Clean picks a large
1128 scale that has no constraints from the data (the UV hole issue).
1129 restoration e.
1131 Construct a restored image : imagename.image by convolving the model
1132 image with a clean beam and adding the residual image to the result.
1133 If a restoringbeam is specified, the residual image is also
1134 smoothed to that target resolution before adding it in.
1136 If a .model does not exist, it will make an empty one and create
1137 the restored image from the residuals ( with additional smoothing if needed ).
1138 With algorithm='mtmfs', this will construct Taylor coefficient maps from
1139 the residuals and compute .alpha and .alpha.error.
1140 restoringbeam ize to use.
1142 - restoringbeam='' or [''].
1143 A Gaussian fitted to the PSF main lobe (separately per image plane).
1145 - restoringbeam='10.0arcsec'.
1146 Use a circular Gaussian of this width for all planes.
1148 - restoringbeam=['8.0arcsec','10.0arcsec','45deg'].
1149 Use this elliptical Gaussian for all planes.
1151 - restoringbeam='common'.
1152 Automatically estimate a common beam shape/size appropriate for
1153 all planes. This option can be used when the beam shape is different as a function of frequency, and will smooth all planes to a single beam, defined by the largest beam in the cube.
1155 Note : For any restoring beam different from the native resolution
1156 the model image is convolved with the beam and added to
1157 residuals that have been convolved to the same target resolution.
1158 pbcor the output restored image.
1160 A new image with extension .image.pbcor will be created from
1161 the evaluation of .image / .pb for all pixels above the specified pblimit.
1163 Note : Stand-alone PB-correction can be triggered by re-running
1164 tclean with the appropriate imagename and with
1165 niter=0, calcpsf=False, calcres=False, pbcor=True, vptable='vp.tab'
1166 ( where vp.tab is the name of the vpmanager file;
1167 see the inline help for the 'vptable' parameter ). Alternatively, task impbcor can be used for primary beam correction using the .image and .pb files.
1169 Note : For deconvolver='mtmfs', pbcor will divide each Taylor term image by the .tt0 average PB.
1170 For all gridders, this calculation is accurate for small fractional bandwidths.
1172 For large fractional bandwidths, please use one of the following options.
1174 (a) For single pointings, run the tclean task with specmode='mfs', deconvolver='mtmfs',
1175 and gridder='standard' with pbcor=True or False.
1176 If a PB-corrected spectral index is required,
1177 please use the widebandpbcor task to apply multi-tern PB-correction.
1179 (b) For mosaics, run tclean task with specmode='mfs', deconvolver='mtmfs',
1180 and gridder='awproject' , wbawp=True, conjbeams=True, with pbcor=True.
1181 This option applies wideband PB correction as part of the gridding step and
1182 pbcor=True will be accurate because the spectral index map will already
1183 be PB-corrected.
1185 (c) For mosaics, run tclean with specmode='mvc', deconvolver='mtmfs',
1186 and gridder='mosaic' or 'awp2' with pbcor=True.
1187 This option applies wideband PB-correction to channelized residual images
1188 prior to the minor cycle and pbcor=True will be accurate because the spectral
1189 index map will already be PB-corrected.
1191 Note : Frequency-dependent PB corrections are typically required for full-band imaging with the VLA.
1192 Wideband PB corrections are required when the amplitude of the
1193 brightest source is known accurately enough to be sensitive
1194 to the difference in the PB gain between the upper and lower
1195 end of the band at its location. As a guideline, the artificial spectral
1196 index due to the PB is -1.4 at the 0.5 gain level and less than -0.2
1197 at the 0.9 gain level at the middle frequency )
1198 outlierfile Name of outlier-field image definitions.
1200 A text file containing sets of parameter=value pairs,
1201 one set per outlier field.
1203 Example : outlierfile='outs.txt'
1205 Contents of outs.txt :
1207 imagename=tst1
1208 nchan=1
1209 imsize=[80,80]
1210 cell=[8.0arcsec,8.0arcsec]
1211 phasecenter=J2000 19:58:40.895 +40.55.58.543
1212 mask=circle[[40pix,40pix],10pix]
1214 imagename=tst2
1215 nchan=1
1216 imsize=[100,100]
1217 cell=[8.0arcsec,8.0arcsec]
1218 phasecenter=J2000 19:58:40.895 +40.56.00.000
1219 mask=circle[[60pix,60pix],20pix]
1221 The following parameters are currently allowed to be different between
1222 the main field and the outlier fields (i.e. they will be recognized if found
1223 in the outlier text file). If a parameter is not listed, the value is picked from
1224 what is defined in the main task input.
1226 imagename, imsize, cell, phasecenter, startmodel, mask
1227 specmode, nchan, start, width, nterms, reffreq,
1228 gridder, deconvolver, wprojplanes.
1230 Note : 'specmode' is an option, so combinations of mfs and cube
1231 for different image fields, for example, are supported.
1232 'deconvolver' and 'gridder' are also options that allow different
1233 imaging or deconvolution algorithm per image field.
1235 For example, multiscale with wprojection and 16 w-term planes
1236 on the main field and mtmfs with nterms=3 and wprojection
1237 with 64 planes on a bright outlier source for which the frequency
1238 dependence of the primary beam produces a strong effect that
1239 must be modeled. The traditional alternative to this approach is
1240 to first image the outlier, subtract it out of the data (uvsub) and
1241 then image the main field.
1242 weighting Weighting scheme (natural,uniform,briggs,superuniform,radial, briggsabs, briggsbwtaper).
1244 During gridding of the dirty or residual image, each visibility value is
1245 multiplied by a weight before it is accumulated on the uv-grid.
1246 The PSF's uv-grid is generated by gridding only the weights (weightgrid).
1248 weighting='natural' : Gridding weights are identical to the data weights
1249 from the MS. For visibilities with similar data weights,
1250 the weightgrid will follow the sample density
1251 pattern on the uv-plane. This weighting scheme
1252 provides the maximum imaging sensitivity at the
1253 expense of a PSF with possibly wider main lobes and high sidelobes.
1254 It is most appropriate for detection experiments
1255 where sensitivity is most important.
1257 weighting='uniform' : Gridding weights per visibility data point are the
1258 original data weights divided by the total weight of
1259 all data points that map to the same uv grid cell :
1260 ' data_weight / total_wt_per_cell '.
1262 The weightgrid is as close to flat as possible resulting
1263 in a PSF with a narrow main lobe and suppressed
1264 sidelobes. However, since heavily sampled areas of
1265 the uv-plane get down-weighted, the imaging
1266 sensitivity is not as high as with natural weighting.
1267 It is most appropriate for imaging experiments where
1268 a well behaved PSF can help the reconstruction.
1270 weighting='briggs' : Gridding weights per visibility data point are given by
1271 'data_weight / ( A \* total_wt_per_cell + B ) ' where
1272 A and B vary according to the 'robust' parameter.
1274 robust = -2.0 maps to A=1,B=0 or uniform weighting.
1275 robust = +2.0 maps to natural weighting.
1276 (robust=0.5 is equivalent to robust=0.0 in AIPS IMAGR.)
1278 Robust/Briggs weighting generates a PSF that can
1279 vary smoothly between 'natural' and 'uniform' and
1280 allow customized trade-offs between PSF shape and
1281 imaging sensitivity.
1282 weighting='briggsabs' : Experimental option.
1283 Same as Briggs except the formula is different A=
1284 robust\*robust and B is dependent on the
1285 noise per visibility estimated. Giving noise='0Jy'
1286 is a not a reasonable option.
1287 In this mode (or formula) robust values
1288 from -2.0 to 0.0 only make sense (2.0 and
1289 -2.0 will get the same weighting)
1291 weighting='superuniform' : This is similar to uniform weighting except that
1292 the total_wt_per_cell is replaced by the
1293 total_wt_within_NxN_cells around the uv cell of
1294 interest. N=7 is the default (when the
1295 parameter 'npixels' is set to 0 with 'superuniform')
1297 This method tends to give a PSF with inner
1298 sidelobes that are suppressed as in uniform
1299 weighting but with far-out sidelobes closer to
1300 natural weighting. The peak sensitivity is also
1301 closer to natural weighting.
1303 weighting='radial' : Gridding weights are given by ' data_weight \* uvdistance '
1304 This method approximately minimizes rms sidelobes
1305 for an east-west synthesis array.
1307 weighting='briggsbwtaper' : A modified version of Briggs weighting for cubes where an inverse uv taper,
1308 which is proportional to the fractional bandwidth of the entire cube,
1309 is applied per channel. The objective is to modify cube (perchanweightdensity = True)
1310 imaging weights to have a similar density to that of the continuum imaging weights.
1311 This is currently an experimental weighting scheme being developed for ALMA.
1313 For more details on weighting please see Chapter3
1314 of Dan Briggs' thesis (http://www.aoc.nrao.edu/dissertations/dbriggs)
1315 robust Robustness parameter for Briggs weighting.
1317 robust = -2.0 maps to uniform weighting.
1318 robust = +2.0 maps to natural weighting.
1319 (robust=0.5 is equivalent to robust=0.0 in AIPS IMAGR.)
1320 noise noise parameter for briggs abs mode weighting
1321 npixels Number of pixels to determine uv-cell size for super-uniform weighting
1322 (0 defaults to -/+ 3 pixels).
1324 npixels -- uv-box used for weight calculation
1325 a box going from -npixel/2 to +npixel/2 on each side
1326 around a point is used to calculate weight density.
1328 npixels=2 goes from -1 to +1 and covers 3 pixels on a side.
1330 npixels=0 implies a single pixel, which does not make sense for
1331 superuniform weighting. Therefore, for 'superuniform'
1332 weighting, if npixels=0 it will be forced to 6 (or a box
1333 of -3pixels to +3pixels) to cover 7 pixels on a side.
1334 uvtaper uv-taper on outer baselines in uv-plane.
1336 Apply a Gaussian taper in addition to the weighting scheme specified
1337 via the 'weighting' parameter. Higher spatial frequencies are weighted
1338 down relative to lower spatial frequencies to suppress artifacts
1339 arising from poorly sampled areas of the uv-plane. It is equivalent to
1340 smoothing the PSF obtained by other weighting schemes and can be
1341 specified either as the HWHM of a Gaussian in uv-space (eg. units of lambda)
1342 or as the FWHM of a Gaussian in the image domain (eg. angular units like arcsec).
1344 uvtaper = [bmaj, bmin, bpa].
1346 Note : FWHM_uv_lambda = (4 log2) / ( pi * FWHM_lm_radians ).
1348 A FWHM_lm of 100.000 arcsec maps to a HWHM_uv of 910.18 lambda.
1349 A FWHM_lm of 1 arcsec maps to a HWHM_uv of 91 klambda.
1351 default: uvtaper=[]; no Gaussian taper applied.
1352 example: uvtaper=['5klambda'] circular taper of HWHM=5 kilo-lambda.
1353 uvtaper=['5klambda','3klambda','45.0deg'] uv-domain HWHM.
1354 uvtaper=['50arcsec','30arcsec','30.0deg'] : image domain FWHM.
1355 uvtaper=['10arcsec'] : image domain FWHM.
1356 uvtaper=['300.0'] default units are lambda in aperture plane.
1357 niter Maximum number of iterations.
1359 A stopping criterion based on total iteration count.
1360 Currently the parameter type is defined as an integer therefore the integer value
1361 larger than 2147483647 will not be set properly as it causes an overflow.
1363 Iterations are typically defined as the selecting one flux component
1364 and partially subtracting it out from the residual image.
1366 niter=0 : Do only the initial major cycle (make dirty image, psf, pb, etc).
1368 niter larger than zero : Run major and minor cycles.
1370 Note : Global stopping criteria vs major-cycle triggers.
1372 In addition to global stopping criteria, the following rules are
1373 used to determine when to terminate a set of minor cycle iterations
1374 and trigger major cycles [derived from Cotton-Schwab Clean, 1984].
1376 'cycleniter' : controls the maximum number of iterations per image
1377 plane before triggering a major cycle.
1378 'cyclethreshold' : Automatically computed threshold related to the
1379 max sidelobe level of the PSF and peak residual.
1380 Divergence, detected as an increase of 10% in peak residual from the
1381 minimum so far (during minor cycle iterations).
1383 The first criterion to be satisfied takes precedence.
1385 Note : Iteration counts for cubes or multi-field images :
1386 For images with multiple planes (or image fields) on which the
1387 deconvolver operates in sequence, iterations are counted across
1388 all planes (or image fields). The iteration count is compared with
1389 'niter' only after all channels/planes/fields have completed their
1390 minor cycles and exited either due to 'cycleniter' or 'cyclethreshold'.
1391 Therefore, the actual number of iterations reported in the logger
1392 can sometimes be larger than the user specified value in 'niter'.
1393 For example, with niter=100, cycleniter=20,nchan=10,threshold=0,
1394 a total of 200 iterations will be done in the first set of minor cycles
1395 before the total is compared with niter=100 and it exits.
1397 Note : Additional global stopping criteria include:
1398 - no change in peak residual across two major cycles.
1399 - a 50% or more increase in peak residual across one major cycle.
1400 gain Loop gain.
1402 Fraction of the source flux to subtract out of the residual image
1403 for the CLEAN algorithm and its variants.
1405 A low value (0.2 or less) is recommended when the sky brightness
1406 distribution is not well represented by the basis functions used by
1407 the chosen deconvolution algorithm. A higher value can be tried when
1408 there is a good match between the true sky brightness structure and
1409 the basis function shapes. For example, for extended emission,
1410 multiscale clean with an appropriate set of scale sizes will tolerate
1411 a higher loop gain than Clark clean.
1412 threshold Stopping threshold (number in units of Jy, or string).
1414 A global stopping threshold that the peak residual (within clean mask)
1415 across all image planes is compared to.
1417 threshold = 0.005 : 5mJy
1418 threshold = '5.0mJy'
1420 Note : A 'cyclethreshold' is internally computed and used as a major cycle
1421 trigger. It is related to what fraction of the PSF can be reliably
1422 used during minor cycle updates of the residual image. By default
1423 the minor cycle iterations terminate once the peak residual reaches
1424 the first sidelobe level of the brightest source.
1426 'cyclethreshold' is computed as follows using the settings in
1427 parameters 'cyclefactor','minpsffraction','maxpsffraction','threshold' :
1429 psf_fraction = max_psf_sidelobe_level \* 'cyclefactor'
1430 psf_fraction = max(psf_fraction, 'minpsffraction');
1431 psf_fraction = min(psf_fraction, 'maxpsffraction');
1432 cyclethreshold = peak_residual \* psf_fraction
1433 cyclethreshold = max( cyclethreshold, 'threshold' )
1435 If nsigma is set (>0.0), the N-sigma threshold is calculated (see
1436 the description under nsigma), then cyclethreshold is further modified as,
1438 cyclethreshold = max( cyclethreshold, nsgima_threshold ).
1441 'cyclethreshold' is made visible and editable only in the
1442 interactive GUI when tclean is run with interactive=True.
1443 nsigma Multiplicative factor for rms-based threshold stopping.
1445 N-sigma threshold is calculated as nsigma \* rms value per image plane determined
1446 from a robust statistics. For nsigma > 0.0, in a minor cycle, a maximum of the two values,
1447 the N-sigma threshold and cyclethreshold, is used to trigger a major cycle
1448 (see also the descreption under 'threshold').
1449 Set nsigma=0.0 to preserve the previous tclean behavior without this feature.
1450 The top level parameter, fastnoise is relevant for the rms noise calculation which is used
1451 to determine the threshold.
1453 The parameter 'nsigma' may be an int, float, or a double.
1454 cycleniter Maximum number of minor-cycle iterations (per plane) before triggering
1455 a major cycle.
1457 For example, for a single plane image, if niter=100 and cycleniter=20,
1458 there will be 5 major cycles after the initial one (assuming there is no
1459 threshold based stopping criterion). At each major cycle boundary, if
1460 the number of iterations left over (to reach niter) is less than cycleniter,
1461 it is set to the difference.
1463 Note : cycleniter applies per image plane, even if cycleniter x nplanes
1464 gives a total number of iterations greater than 'niter'. This is to
1465 preserve consistency across image planes within one set of minor
1466 cycle iterations.
1467 cyclefactor Scaling on PSF sidelobe level to compute the minor-cycle stopping threshold.
1469 Please refer to the Note under the documentation for 'threshold' that
1470 discussed the calculation of 'cyclethreshold'.
1472 cyclefactor=1.0 results in a cyclethreshold at the first sidelobe level of
1473 the brightest source in the residual image before the minor cycle starts.
1475 cyclefactor=0.5 allows the minor cycle to go deeper.
1476 cyclefactor=2.0 triggers a major cycle sooner.
1477 minpsffraction PSF fraction that marks the max depth of cleaning in the minor cycle.
1479 Please refer to the Note under the documentation for 'threshold' that
1480 discussed the calculation of 'cyclethreshold'.
1482 For example, minpsffraction=0.5 will stop cleaning at half the height of
1483 the peak residual and trigger a major cycle earlier.
1484 maxpsffraction PSF fraction that marks the minimum depth of cleaning in the minor cycle.
1486 Please refer to the Note under the documentation for 'threshold' that
1487 discussed the calculation of 'cyclethreshold'.
1489 For example, maxpsffraction=0.8 will ensure that at least the top 20
1490 percent of the source will be subtracted out in the minor cycle even if
1491 the first PSF sidelobe is at the 0.9 level (an extreme example), or if the
1492 cyclefactor is set too high for anything to get cleaned.
1493 nmajor The nmajor parameter limits the number of minor and major cycle sets
1494 that tclean executes. It is defined as the number of major cycles after the
1495 initial set of minor cycle iterations. In other words, the count of nmajor does
1496 not include the initial residual calculation that occurs when calcres=True.
1498 A setting of nmajor=-1 implies no limit (default -1).
1499 A setting of nmajor=0 implies nothing other than the initial residual calculation
1500 A setting of nmajor>0 imples that nmajor sets of minor and major cycles will
1501 be done in addition to the initial residual calculation.
1503 If the major cycle limit is reached, stopcode 9 will be returned. Other stopping
1504 criteria (such as threshold) could cause tclean to stop in fewer than this
1505 number of major cycles. If tclean reaches another stopping criteria, first
1506 or at the same time as nmajor, then that stopcode will be returned instead.
1508 Note however that major cycle ids in the log messages as well as in the return
1509 dictionary do begin with 1 for the initial residual calculation, when it exists.
1511 Example 1 : A tclean run with 'nmajor=5' and 'calcres=True' will iterate for
1512 5 major cycles (not counting the initial residual calculation). But, the return
1513 dictionary will show 'nmajordone:6'. If 'calcres=False', then the return
1514 dictionary will show 'nmajordone:5'.
1516 Example 2 : For both the following cases, there will be a printout in the logs
1517 "Running Major Cycle 1" and the return value will include "nmajordone: 1",
1518 however there is a difference in the purpose of the major cycle and the
1519 number of minor cycles executed:
1520 Case 1; nmajor=0, calcres=True: The major cycle done is for the creation
1521 of the residual, and no minor cycles are executed.
1522 Case 2; nmajor=1, calcres=False: The major cycle is done as part of the
1523 major/minor cycle loop, and 1 minor cycle will be executed.
1524 usemask Type of mask(s) to be used for deconvolution.
1526 user: (default) mask image(s) or user specified region file(s) or string CRTF expression(s).
1527 subparameters: mask, pbmask.
1528 pb: primary beam mask.
1529 subparameter: pbmask.
1531 Example: usemask="pb", pbmask=0.2.
1532 Construct a mask at the 0.2 pb gain level.
1533 (Currently, this option will work only with
1535 gridders that produce .pb (i.e. mosaic, awp2 and awproject)
1536 or if an externally produced .pb image exists on disk)
1539 auto-multithresh : auto-masking by multiple thresholds for deconvolution.
1540 subparameters : sidelobethreshold, noisethreshold, lownoisethreshold, negativethrehsold, smoothfactor,
1541 minbeamfrac, cutthreshold, pbmask, growiterations, dogrowprune, minpercentchange, verbose.
1542 Additional top level parameter relevant to auto-multithresh: fastnoise.
1544 if pbmask is >0.0, the region outside the specified pb gain level is excluded from
1545 image statistics in determination of the threshold.
1550 Note: By default the intermediate mask generated by automask at each deconvolution cycle
1551 is over-written in the next cycle but one can save them by setting
1552 the environment variable, SAVE_ALL_AUTOMASKS="true".
1553 (e.g. in the CASA prompt, os.environ['SAVE_ALL_AUTOMASKS']="true" )
1554 The saved CASA mask image name will be imagename.mask.autothresh#, where
1555 # is the iteration cycle number.
1556 mask Mask (a list of image name(s) or region file(s) or region string(s).
1559 The name of a CASA image or region file or region string that specifies
1560 a 1/0 mask to be used for deconvolution. Only locations with value 1 will
1561 be considered for the centers of flux components in the minor cycle.
1562 If regions specified fall completely outside of the image, tclean will throw an error.
1564 Manual mask options/examples :
1566 mask='xxx.mask' : Use this CASA image named xxx.mask and containing
1567 ones and zeros as the mask.
1568 If the mask is only different in spatial coordinates from what is being made
1569 it will be resampled to the target coordinate system before being used.
1570 The mask has to have the same shape in velocity and Stokes planes
1571 as the output image. Exceptions are single velocity and/or single
1572 Stokes plane masks. They will be expanded to cover all velocity and/or
1573 Stokes planes of the output cube.
1575 [ Note : If an error occurs during image resampling or
1576 if the expected mask does not appear, please try
1577 using tasks 'imregrid' or 'makemask' to resample
1578 the mask image onto a CASA image with the target
1579 shape and coordinates and supply it via the 'mask'
1580 parameter. ]
1583 mask='xxx.crtf' : A text file with region strings and the following on the first line
1584 ( #CRTFv0 CASA Region Text Format version 0 )
1585 This is the format of a file created via the viewer's region
1586 tool when saved in CASA region file format.
1588 mask='circle[[40pix,40pix],10pix]' : A CASA region string.
1590 mask=['xxx.mask','xxx.crtf', 'circle[[40pix,40pix],10pix]'] : a list of masks.
1596 Note : Mask images for deconvolution must contain 1 or 0 in each pixel.
1597 Such a mask is different from an internal T/F mask that can be
1598 held within each CASA image. These two types of masks are not
1599 automatically interchangeable, so please use the makemask task
1600 to copy between them if you need to construct a 1/0 based mask
1601 from a T/F one.
1603 Note : Work is in progress to generate more flexible masking options and
1604 enable more controls.
1605 pbmask Sub-parameter for usemask: primary beam mask.
1607 Examples : pbmask=0.0 (default, no pb mask).
1608 pbmask=0.2 (construct a mask at the 0.2 pb gain level).
1609 sidelobethreshold Sub-parameter for "auto-multithresh": mask threshold based on sidelobe levels: sidelobethreshold \* max_sidelobe_level \* peak residual.
1610 noisethreshold Sub-parameter for "auto-multithresh": mask threshold based on the noise level: noisethreshold \* rms + location (=median).
1612 The rms is calculated from the median absolute deviation (MAD), with rms = 1.4826\*MAD.
1613 lownoisethreshold Sub-parameter for "auto-multithresh": mask threshold to grow previously masked regions via binary dilation: lownoisethreshold \* rms in residual image + location (=median).
1615 The rms is calculated from the median absolute deviation (MAD), with rms = 1.4826\*MAD.
1616 negativethreshold Sub-parameter for "auto-multithresh": mask threshold for negative features: -1.0* negativethreshold \* rms + location(=median).
1618 The rms is calculated from the median absolute deviation (MAD), with rms = 1.4826\*MAD.
1619 smoothfactor Sub-parameter for "auto-multithresh": smoothing factor in a unit of the beam.
1620 minbeamfrac Sub-parameter for "auto-multithresh": minimum beam fraction in size to prune masks smaller than mimbeamfrac \* beam
1621 <=0.0 : No pruning
1622 cutthreshold Sub-parameter for "auto-multithresh": threshold to cut the smoothed mask to create a final mask: cutthreshold \* peak of the smoothed mask.
1623 growiterations Sub-parameter for "auto-multithresh": Maximum number of iterations to perform using binary dilation for growing the mask.
1624 dogrowprune Experimental sub-parameter for "auto-multithresh": Do pruning on the grow mask.
1625 minpercentchange If the change in the mask size in a particular channel is less than minpercentchange, stop masking that channel in subsequent cycles. This check is only applied when noise based threshold is used and when the previous clean major cycle had a cyclethreshold value equal to the clean threshold. Values equal to -1.0 (or any value less than 0.0) will turn off this check (the default). Automask will still stop masking if the current channel mask is an empty mask and the noise threshold was used to determine the mask.
1626 verbose he summary of automasking at the end of each automasking process
1627 is printed in the logger. Following information per channel will be listed in the summary.
1629 chan: channel number.
1630 masking?: F - stop updating automask for the subsequent iteration cycles.
1631 RMS: robust rms noise.
1632 peak: peak in residual image.
1633 thresh_type: type of threshold used (noise or sidelobe).
1634 thresh_value: the value of threshold used.
1635 N_reg: number of the automask regions.
1636 N_pruned: number of the automask regions removed by pruning.
1637 N_grow: number of the grow mask regions.
1638 N_grow_pruned: number of the grow mask regions removed by pruning.
1639 N_neg_pix: number of pixels for negative mask regions.
1641 Note that for a large cube, extra logging may slow down the process.
1642 fastnoise Only relevant when automask (user='multi-autothresh') and/or n-sigma stopping threshold (nsigma>0.0) are/is used. If it is set to True, a simpler but faster noise calucation is used.
1643 In this case, the threshold values are determined based on classic statistics (using all
1644 unmasked pixels for the calculations).
1646 If it is set to False, the new noise calculation
1647 method is used based on pre-existing mask.
1649 Case 1: no exiting mask.
1650 Calculate image statistics using Chauvenet algorithm.
1652 Case 2: there is an existing mask.
1653 Calculate image statistics by classical method on the region
1654 outside the mask and inside the primary beam mask.
1656 In all cases above RMS noise is calculated from the median absolute deviation (MAD).
1657 restart images (and start from an existing model image)
1658 or automatically increment the image name and make a new image set.
1660 True : Re-use existing images. If imagename.model exists the subsequent
1661 run will start from this model (i.e. predicting it using current gridder
1662 settings and starting from the residual image). Care must be taken
1663 when combining this option with startmodel. Currently, only one or
1664 the other can be used.
1666 startmodel='', imagename.model exists :
1667 - Start from imagename.model.
1668 startmodel='xxx', imagename.model does not exist :
1669 - Start from startmodel.
1670 startmodel='xxx', imagename.model exists :
1671 - Exit with an error message requesting the user to pick
1672 only one model. This situation can arise when doing one
1673 run with startmodel='xxx' to produce an output
1674 imagename.model that includes the content of startmodel,
1675 and wanting to restart a second run to continue deconvolution.
1676 Startmodel should be set to '' before continuing.
1678 If any change in the shape or coordinate system of the image is
1679 desired during the restart, please change the image name and
1680 use the startmodel (and mask) parameter(s) so that the old model
1681 (and mask) can be regridded to the new coordinate system before starting.
1683 False : A convenience feature to increment imagename with '_1', '_2',
1684 etc as suffixes so that all runs of tclean are fresh starts (without
1685 having to change the imagename parameter or delete images).
1687 This mode will search the current directory for all existing
1688 imagename extensions, pick the maximum, and adds 1.
1689 For imagename='try' it will make try.psf, try_2.psf, try_3.psf, etc.
1691 This also works if you specify a directory name in the path :
1692 imagename='outdir/try'. If './outdir' does not exist, it will create it.
1693 Then it will search for existing filenames inside that directory.
1695 If outlier fields are specified, the incrementing happens for each
1696 of them (since each has its own 'imagename'). The counters are
1697 synchronized across imagefields, to make it easier to match up sets
1698 of output images. It adds 1 to the 'max id' from all outlier names
1699 on disk. So, if you do two runs with only the main field
1700 (imagename='try'), and in the third run you add an outlier with
1701 imagename='outtry', you will get the following image names
1702 for the third run : 'try_3' and 'outtry_3' even though
1703 'outry' and 'outtry_2' have not been used.
1704 savemodel Options to save model visibilities (none, virtual, modelcolumn).
1706 Often, model visibilities must be created and saved in the MS
1707 to be later used for self-calibration (or to just plot and view them).
1709 none : Do not save any model visibilities in the MS. The MS is opened
1710 in readonly mode.
1712 Model visibilities can be predicted in a separate step by
1713 restarting tclean with niter=0,savemodel=virtual or modelcolumn
1714 and not changing any image names so that it finds the .model on
1715 disk (or by changing imagename and setting startmodel to the
1716 original imagename).
1718 virtual : In the last major cycle, save the image model and state of the
1719 gridder used during imaging within the SOURCE subtable of the
1720 MS. Images required for de-gridding will also be stored internally.
1721 All future references to model visibilities will activate the
1722 (de)gridder to compute them on-the-fly. This mode is useful
1723 when the dataset is large enough that an additional model data
1724 column on disk may be too much extra disk I/O, when the
1725 gridder is simple enough that on-the-fly recomputing of the
1726 model visibilities is quicker than disk I/O.
1727 For e.g. that gridder='awproject' and 'awp2' does not support virtual model.
1729 modelcolumn : In the last major cycle, save predicted model visibilities
1730 in the MODEL_DATA column of the MS. This mode is useful when
1731 the de-gridding cost to produce the model visibilities is higher
1732 than the I/O required to read the model visibilities from disk.
1733 This mode is currently required for gridder='awproject' and 'awp2'.
1734 This mode is also required for the ability to later pull out
1735 model visibilities from the MS into a python array for custom
1736 processing.
1738 Note 1 : The imagename.model image on disk will always be constructed
1739 if the minor cycle runs. This savemodel parameter applies only to
1740 model visibilities created by de-gridding the model image.
1742 Note 2 : It is possible for an MS to have both a virtual model
1743 as well as a model_data column, but under normal operation,
1744 the last used mode will get triggered. Use the delmod task to
1745 clear out existing models from an MS if confusion arises.
1746 Note 3: when parallel=True, use savemodel='none'; Other options are not yet ready
1747 for use in parallel. If model visibilities need to be saved (virtual or modelcolumn):
1748 please run tclean in serial mode with niter=0; after the parallel run
1749 calcres Calculate initial residual image.
1751 This parameter controls what the first major cycle does.
1753 calcres=False with niter greater than 0 will assume that
1754 a .residual image already exists and that the minor cycle can
1755 begin without recomputing it.
1757 calcres=False with niter=0 implies that only the PSF will be made
1758 and no data will be gridded.
1760 calcres=True requires that calcpsf=True or that the .psf and .sumwt
1761 images already exist on disk (for normalization purposes).
1763 Usage example : For large runs (or a pipeline scripts) it may be
1764 useful to first run tclean with niter=0 to create
1765 an initial .residual to look at and perhaps make
1766 a custom mask for. Imaging can be resumed
1767 without recomputing it.
1768 calcpsf Calculate PSF
1770 This parameter controls what the first major cycle does.
1772 calcpsf=False will assume that a .psf image already exists
1773 and that the minor cycle can begin without recomputing it.
1774 psfcutoff When the .psf image is created a 2 dimensional Gaussian is fit to the main lobe of the PSF.
1775 Which pixels in the PSF are fitted is determined by psfcutoff.
1776 The default value of psfcutoff is 0.35 and can varied from 0.01 to 0.99.
1777 Fitting algorithm:
1778 - A region of 41 x 41 pixels around the peak of the PSF is compared against the psfcutoff.
1779 Sidelobes are ignored by radially searching from the PSF peak.
1780 - Calculate the bottom left corner (blc) and top right corner (trc) from the points. Expand blc and trc with a number of pixels (5).
1781 - Create a new sub-matrix from blc and trc.
1782 - Interpolate matrix to a target number of points (3001) using CUBIC spline.
1783 - All the non-sidelobe points, in the interpolated matrix, that are above the psfcutoff are used to fit a Gaussian.
1784 A Levenberg-Marquardt algorithm is used.
1785 - If the fitting fails the algorithm is repeated with the psfcutoff decreased (psfcutoff=psfcutoff/1.5).
1786 A message in the log will apear if the fitting fails along with the new value of psfcutoff.
1787 This will be done up to 50 times if fitting fails.
1788 This Gaussian beam is defined by a major axis, minor axis, and position angle.
1789 During the restoration process, this Gaussian beam is used as the Clean beam.
1790 Varying psfcutoff might be useful for producing a better fit for highly non-Gaussian PSFs, however, the resulting fits should be carefully checked.
1791 This parameter should rarely be changed.
1793 (This is not the support size for clark clean.)
1794 parallel Run major cycles in parallel.
1796 Parallel tclean will run only if casa has already been started using mpirun.
1797 Please refer to external resources on high performance computing for details on how to start this on your system.
1799 Example : mpirun -n 3 -xterm 0 `which casa`
1801 Continuum Imaging :
1802 - Data are partitioned (in time) into NProc pieces.
1803 - Gridding/iFT is done separately per partition.
1804 - Images (and weights) are gathered and then normalized.
1805 - One non-parallel minor cycle is run.
1806 - Model image is scattered to all processes.
1807 - Major cycle is done in parallel per partition.
1809 Cube Imaging :
1810 - Data and Image coordinates are partitioned (in freq) into NProc pieces.
1811 - Each partition is processed independently (major and minor cycles).
1812 - All processes are synchronized at major cycle boundaries for convergence checks.
1813 - At the end, cubes from all partitions are concatenated along the spectral axis.
1815 Note 1 : Iteration control for cube imaging is independent per partition.
1816 - There is currently no communication between them to synchronize
1817 information such as peak residual and cyclethreshold. Therefore,
1818 different chunks may trigger major cycles at different levels.
1819 (Proper synchronization of iteration control is work in progress.)
1820 [1;42mRETURNS[1;m void
1822 --------- examples -----------------------------------------------------------
1826 For more information, see the task pages of tclean in CASA Docs:
1828 https://casadocs.readthedocs.io
1834 '''
1836 _info_group = """imaging"""
1837 _info_desc = """Radio Interferometric Image Reconstruction"""
1839 def _tclean( self, *args, **kwargs ):
1840 """ Calls tclean records the arguments in the local history of tclean calls.
1841casatasks/src/private/imagerhelpers/_gclean.py
1842 The full tclean history for this instance can be retrieved via the cmds() method."""
1843 arg_s = ', '.join( map( lambda a: self._history_filter(len(self._exe_cmds), None, repr(a)), args ) )
1844 kw_s = ', '.join( map( lambda kv: self._history_filter(len(self._exe_cmds), kv[0], "%s=%s" % (kv[0],repr(kv[1]))), kwargs.items()) )
1845 if len(arg_s) > 0 and len(kw_s) > 0:
1846 parameters = arg_s + ", " + kw_s
1847 else:
1848 parameters = arg_s + kw_s
1849 self._exe_cmds.append( "tclean( %s )" % parameters )
1850 self._exe_cmds_per_iter[-1] += 1
1851 return tclean( *args, **kwargs )
1853 def _deconvolve( self, *args, **kwargs ):
1854 """ Calls deconvolve records the arguments in the local history of deconvolve calls.
1856 The full deconvolve history for this instance can be retrieved via the cmds() method."""
1857 arg_s = ', '.join( map( lambda a: self._history_filter(len(self._exe_cmds), None, repr(a)), args ) )
1858 kw_s = ', '.join( map( lambda kv: self._history_filter(len(self._exe_cmds), kv[0], "%s=%s" % (kv[0],repr(kv[1]))), kwargs.items()) )
1859 if len(arg_s) > 0 and len(kw_s) > 0:
1860 parameters = arg_s + ", " + kw_s
1861 else:
1862 parameters = arg_s + kw_s
1863 self._exe_cmds.append( "deconvolve( %s )" % parameters )
1864 self._exe_cmds_per_iter[-1] += 1
1865 return deconvolve( *args, **kwargs )
1868 def _deconv_over_fields(self, niter, mask=None, threshold=None, restoration=False, hasit_per_field=None):
1869 """
1870 Loop over each field, deconvolve, merge the return dicts, all that good
1871 stuff. The only configurable option here is self._niter, and
1872 restoreation, which is different depending on whether we want to
1873 deconvolve (niter > 0) or update the mask (niter = 0), or restore
1874 (niter=0, restoration=True).
1876 The niter passed in is tracked by self._niter for when we need to do deconvolution,
1877 and is set to 0 otherwise.
1879 The other parameters that come in from the GUI are used directly within this function,
1880 such as self._cyclethreshold, self._nsigma, self._gain, etc.
1881 """
1883 initialize_merge_dict = False
1884 counter = 0
1886 # Loop over each image (outliers included) and deconvolve.
1887 # This is not necessary for tclean, since we want a joint gridding of all the images
1888 for fidx, fieldid in enumerate(self._paramlist.allimpars.keys()):
1889 # Check if convergence info exists
1890 # If it exists, skip field if it has converged
1891 if hasit_per_field is not None and restoration == False:
1892 if hasit_per_field[fidx]:
1893 # If convergence is hit, skip this field
1894 continue
1895 elif restoration == True:
1896 # If restoration is True, we always want to run the deconvolution
1897 # even if convergence is hit, so we can restore the image.
1898 # Force niter = 0 regardless of input
1899 niter = 0
1900 pass
1902 if threshold is None:
1903 # User input threshold over all fields from initial iclean call
1904 this_threshold = self._threshold
1905 else:
1906 # Calculated threshold per field
1907 this_threshold = threshold[fidx]
1909 tmpiterpars = self._paramlist.iterpars
1910 tmpdecpars = self._paramlist.alldecpars[fieldid]
1912 thismask = tmpdecpars['mask'] if mask is None else mask
1914 deconv_ret = self._deconvolve(imagename=tmpdecpars['imagename'], startmodel=tmpdecpars['startmodel'],
1915 deconvolver=tmpdecpars['deconvolver'], scales=tmpdecpars['scales'], nterms=tmpdecpars['nterms'],
1916 smallscalebias=tmpdecpars['scalebias'], restoration=restoration, restoringbeam=tmpdecpars['restoringbeam'],
1917 niter = niter, gain=self._gain, threshold=this_threshold, nsigma=self._nsigma,
1918 interactive = False, fullsummary=True, fastnoise=tmpdecpars['fastnoise'], usemask=tmpdecpars['usemask'],
1919 mask = thismask, pbmask=tmpdecpars['pbmask'], sidelobethreshold=self._sidelobethreshold,
1920 noisethreshold=self._noisethreshold, lownoisethreshold=self._lownoisethreshold,
1921 negativethreshold=self._negativethreshold, smoothfactor=tmpdecpars['smoothfactor'],
1922 minbeamfrac=self._minbeamfrac, cutthreshold=tmpdecpars['cutthreshold'],
1923 growiterations=tmpdecpars['growiterations'], dogrowprune=self._dogrowprune,
1924 verbose=tmpdecpars['verbose'])
1927 # Since we are looping through field_id externally, deconvolve
1928 # always thinks it's imaging field 0 Fix this manually so when we
1929 # do the concat below we are not over-writing records
1930 if fidx > 0 and self._nfields > 1:
1931 deconv_ret['summaryminor'][fidx] = deconv_ret['summaryminor'].pop(0)
1932 initialize_merge_dict = True
1934 if self._nfields > 1 and counter == 0:
1935 deconv_ret_merge = copy.deepcopy(deconv_ret)
1937 # Only start merging after the first non-converged field has been processed
1938 if counter > 0:
1939 deconv_ret_merge = ImagingDict.concat_multifield_deconv(deconv_ret_merge, deconv_ret)
1941 counter += 1
1943 # Move the merged dict back to the original var name
1944 if self._nfields > 1 and initialize_merge_dict:
1945 deconv_ret = copy.deepcopy(deconv_ret_merge)
1947 return deconv_ret
1951 def _remove_tree( self, directory ):
1952 if os.path.isdir(directory):
1953 shutil.rmtree(directory)
1954 self._exe_cmds.append( f'''shutil.rmtree( {repr(directory)} )''' )
1955 self._exe_cmds_per_iter[-1] += 1
1957 def _log( self, history=False ):
1958 """ Returns the history of all tclean calls for this instance. If ``history``
1959 is set to True then the full history will be returned, otherwise the commands
1960 executed for generating the latest result are returned.
1961 """
1963 if history:
1964 return self._exe_cmds
1965 else:
1966 if self._exe_cmds_per_iter[-1] > 0:
1967 # Return the last N commands
1968 return self._exe_cmds[-self._exe_cmds_per_iter[-1]:]
1969 else:
1970 # If convergence is hit, no commands were run so return nothing
1971 return []
1974 def update( self, msg ):
1975 """ Interactive clean parameters update.
1977 Args:
1978 msg: dict with possible keys 'niter', 'cycleniter', 'nmajor', 'threshold', 'cyclefactor'
1980 Returns:
1981 stopcode : Stop code in case of error (-1 on error, 0 if no error), int
1982 stopdesc : Exception error message, str
1983 """
1984 if 'niter' in msg:
1985 try:
1986 self._niter = int(msg['niter'])
1987 if self._niter < -1:
1988 return -1, f"niter must be >= -1"
1989 except ValueError as err:
1990 return -1, "niter must be an integer"
1992 if 'cycleniter' in msg:
1993 try:
1994 self._cycleniter = int(msg['cycleniter'])
1995 if self._cycleniter < -1:
1996 return -1, f"cycleniter must be >= -1"
1997 except ValueError:
1998 return -1, "cycleniter must be an integer"
2000 if 'nmajor' in msg:
2001 try:
2002 self._nmajor = int(msg['nmajor'])
2003 if self._nmajor < -1:
2004 return -1, f"nmajor must be >= -1"
2005 except ValueError as e:
2006 return -1, "nmajor must be an integer"
2008 if 'threshold' in msg:
2009 try:
2010 self._threshold = float(msg['threshold'])
2011 if self._threshold < 0:
2012 return -1, f"threshold must be >= 0"
2013 except ValueError:
2014 if isinstance(msg['threshold'], str) and "jy" in msg['threshold'].lower():
2015 self._threshold_to_float(msg['threshold']) # Convert str to float
2016 else:
2017 return -1, f"threshold must be a number, or a number with units (Jy/mJy/uJy)"
2020 if 'nsigma' in msg:
2021 try:
2022 self._nsigma = float(msg['nsigma'])
2023 if self._nsigma < 0:
2024 return -1, f"nsigma must be >= 0"
2025 except ValueError:
2026 return -1, "nsigma must be a number"
2029 if 'gain' in msg:
2030 try:
2031 self._gain = float(msg['gain'])
2032 if self._gain <= 0:
2033 return -1, f"gain must be > 0"
2034 except ValueError:
2035 return -1, "gain must be a number"
2038 if 'cyclefactor' in msg:
2039 try:
2040 self._cyclefactor = float(msg['cyclefactor'])
2041 if self._cyclefactor <= 0:
2042 return -1, f"cyclefactor must be > 0"
2043 except ValueError:
2044 return -1, "cyclefactor must be a number"
2046 ### Automasking parameters
2048 if 'dogrowprune' in msg:
2049 try:
2050 self._dogrowprune = bool(msg['dogrowprune'])
2051 except ValueError:
2052 return -1, "dogrowprune must be a boolean"
2054 if 'noisethreshold' in msg:
2055 try:
2056 self._noisethreshold = float(msg['noisethreshold'])
2057 if self._noisethreshold < 0:
2058 return -1, f"noisethreshold must be >= 0"
2059 except ValueError:
2060 return -1, "noisethreshold must be a number"
2062 if 'sidelobethreshold' in msg:
2063 try:
2064 self._sidelobethreshold = float(msg['sidelobethreshold'])
2065 if self._sidelobethreshold < 0:
2066 return -1, f"sidelobethreshold must be >= 0"
2067 except ValueError:
2068 return -1, "sidelobethreshold must be a number"
2070 if 'lownoisethreshold' in msg:
2071 try:
2072 self._lownoisethreshold = float(msg['lownoisethreshold'])
2073 if self._lownoisethreshold < 0:
2074 return -1, f"lownoisethreshold must be >= 0"
2075 except ValueError:
2076 return -1, "lownoisethreshold must be a number"
2078 if 'minbeamfrac' in msg:
2079 try:
2080 self._minbeamfrac = float(msg['minbeamfrac'])
2081 if self._minbeamfrac < 0:
2082 return -1, f"minbeamfrac must be >= 0"
2083 except ValueError:
2084 return -1, "minbeamfrac must be a number"
2086 if 'negativethreshold' in msg:
2087 try:
2088 self._negativethreshold = float(msg['negativethreshold'])
2089 if self._negativethreshold < 0:
2090 return -1, f"negativethreshold must be >= 0"
2091 except ValueError:
2092 return -1, "negativethreshold must be a number"
2095 return 0, ""
2098 def _threshold_to_float(self, msg=None):
2099 # Convert threshold from string to float if necessary
2100 if msg is not None:
2101 if isinstance(msg, str):
2102 if "mJy" in msg:
2103 self._threshold = float(msg.replace("mJy", "")) / 1e3
2104 elif "uJy" in msg:
2105 self._threshold = float(msg.replace("uJy", "")) / 1e6
2106 elif "Jy" in msg:
2107 self._threshold = float(msg.replace("Jy", ""))
2108 else:
2109 if isinstance(self._threshold, str):
2110 if "mJy" in self._threshold:
2111 self._threshold = float(self._threshold.replace("mJy", "")) / 1e3
2112 elif "uJy" in self._threshold:
2113 self._threshold = float(self._threshold.replace("uJy", "")) / 1e6
2114 elif "Jy" in self._threshold:
2115 self._threshold = float(self._threshold.replace("Jy", ""))
2118 @staticmethod
2119 def _validate(doc, schema):
2120 """ Validate the input parameters against the schema.
2122 Args:
2123 doc: dict, input parameters
2124 schema: dict, schema for the input parameters
2126 Returns:
2127 bool: True if the input parameters are valid, False otherwise
2128 """
2129 return _pc.validate(doc, schema)
2132 def _initialize_convergence_result(self):
2133 """
2134 Initialize the convergence result with the correct number of outlier fields specified.
2135 The structure of the convergence result is a list-of-lists, each element of the outer
2136 list is a field.
2138 The inner list is structured as
2140 self.CR.convergence_result = [None,None,None,None,None,{'field_name':{'chan': None, 'major': None},}]
2141 ^^^^ ^^^^ ^^^^ ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^----->>> convergence info - one dict of dicts per field
2142 | | | | +----->>> Number of global iterations remaining for current run (niterleft)
2143 | | | +---------->>> Number of major cycles remaining for current run (nmajorleft)
2144 | | +---------------->>> major cycles done for current run (nmajordone)
2145 | +------------------>>> tclean stopcode
2146 """
2148 #for nn in range(self._nfields):
2149 # self.CR.convergence_result.append([None, None, None, None, None, { 'chan': None, 'major': None }])
2150 self.CR.convergence_result.append(None, None, None, None, None)
2152 field_dict = {}
2153 for nn in range(self._nfields):
2154 field_dict[self._fieldnames[nn]] = {'chan':None, 'major':None}
2156 self.CR.convergence_result.append(field_dict)
2160 def __init__(self, vis='', selectdata=True, field='', spw='', timerange='', uvrange='', antenna='', scan='', observation='', intent='', datacolumn='corrected', imagename='', imsize=[ int(100) ], cell=[ ], phasecenter='', stokes='I', projection='SIN', startmodel='', specmode='mfs', reffreq='', nchan=int(-1), start='', width='', outframe='LSRK', veltype='radio', restfreq=[ ], interpolation='linear', perchanweightdensity=True, gridder='standard', facets=int(1), psfphasecenter='', wprojplanes=int(1), vptable='', mosweight=True, aterm=True, psterm=False, wbawp=True, conjbeams=False, cfcache='', usepointing=False, computepastep=float(360.0), rotatepastep=float(360.0), pointingoffsetsigdev=[ ], pblimit=float(0.2), normtype='flatnoise', deconvolver='hogbom', scales=[ ], nterms=int(2), smallscalebias=float(0.0), fusedthreshold=float(0.0), largestscale=int(-1), restoration=True, restoringbeam=[ ], pbcor=False, outlierfile='', weighting='natural', robust=float(0.5), noise='1.0Jy', npixels=int(0), uvtaper=[ '' ], niter=int(0), gain=float(0.1), threshold=float(0.0), nsigma=float(0.0), cycleniter=int(-1), cyclefactor=float(1.0), minpsffraction=float(0.05), maxpsffraction=float(0.8), nmajor=int(-1), usemask='user', mask='', pbmask=float(0.0), sidelobethreshold=float(3.0), noisethreshold=float(5.0), lownoisethreshold=float(1.5), negativethreshold=float(0.0), smoothfactor=float(1.0), minbeamfrac=float(0.3), cutthreshold=float(0.01), growiterations=int(75), dogrowprune=True, minpercentchange=float(-1.0), verbose=False, fastnoise=True, restart=True, savemodel='none', calcres=True, calcpsf=True, psfcutoff=float(0.35), parallel=False, history_filter=lambda index, arg, history_value: history_value ):
2162 schema = { 'vis': {'anyof': [{'type': 'cReqPath', 'coerce': _coerce.expand_path}, {'type': 'cReqPathVec', 'coerce': [_coerce.to_list,_coerce.expand_pathvec]}]}, 'selectdata': {'type': 'cBool'}, 'field': {'anyof': [{'type': 'cStr', 'coerce': _coerce.to_str}, {'type': 'cStrVec', 'coerce': [_coerce.to_list,_coerce.to_strvec]}]}, 'spw': {'anyof': [{'type': 'cStr', 'coerce': _coerce.to_str}, {'type': 'cStrVec', 'coerce': [_coerce.to_list,_coerce.to_strvec]}]}, 'timerange': {'anyof': [{'type': 'cStr', 'coerce': _coerce.to_str}, {'type': 'cStrVec', 'coerce': [_coerce.to_list,_coerce.to_strvec]}]}, 'uvrange': {'anyof': [{'type': 'cStr', 'coerce': _coerce.to_str}, {'type': 'cStrVec', 'coerce': [_coerce.to_list,_coerce.to_strvec]}]}, 'antenna': {'anyof': [{'type': 'cStr', 'coerce': _coerce.to_str}, {'type': 'cStrVec', 'coerce': [_coerce.to_list,_coerce.to_strvec]}]}, 'scan': {'anyof': [{'type': 'cStr', 'coerce': _coerce.to_str}, {'type': 'cStrVec', 'coerce': [_coerce.to_list,_coerce.to_strvec]}]}, 'observation': {'anyof': [{'type': 'cStr', 'coerce': _coerce.to_str}, {'type': 'cInt'}]}, 'intent': {'anyof': [{'type': 'cStr', 'coerce': _coerce.to_str}, {'type': 'cStrVec', 'coerce': [_coerce.to_list,_coerce.to_strvec]}]}, 'datacolumn': {'type': 'cStr', 'coerce': _coerce.to_str}, 'imagename': {'anyof': [{'type': 'cInt'}, {'type': 'cStr', 'coerce': _coerce.to_str}, {'type': 'cStrVec', 'coerce': [_coerce.to_list,_coerce.to_strvec]}]}, 'imsize': {'anyof': [{'type': 'cInt'}, {'type': 'cIntVec', 'coerce': [_coerce.to_list,_coerce.to_intvec]}]}, 'cell': {'anyof': [{'type': 'cIntVec', 'coerce': [_coerce.to_list,_coerce.to_intvec]}, {'type': 'cStr', 'coerce': _coerce.to_str}, {'type': 'cFloat', 'coerce': _coerce.to_float}, {'type': 'cStrVec', 'coerce': [_coerce.to_list,_coerce.to_strvec]}, {'type': 'cInt'}, {'type': 'cFloatVec', 'coerce': [_coerce.to_list,_coerce.to_floatvec]}]}, 'phasecenter': {'anyof': [{'type': 'cInt'}, {'type': 'cStr', 'coerce': _coerce.to_str}]}, 'stokes': {'type': 'cStr', 'coerce': _coerce.to_str, 'allowed': [ 'I', 'IQUV', 'UV', 'RRLL', 'IQ', 'V', 'pseudoI', 'QU', 'YY', 'RR', 'Q', 'U', 'IV', 'XX', 'XXYY', 'LL' ]}, 'projection': {'type': 'cStr', 'coerce': _coerce.to_str}, 'startmodel': {'type': 'cVariant', 'coerce': [_coerce.to_variant]}, 'specmode': {'type': 'cStr', 'coerce': _coerce.to_str, 'allowed': [ 'cont', 'cubedata', 'cube', 'cubesource', 'mfs', 'mvc' ]}, 'reffreq': {'type': 'cVariant', 'coerce': [_coerce.to_variant]}, 'nchan': {'type': 'cInt'}, 'start': {'type': 'cVariant', 'coerce': [_coerce.to_variant]}, 'width': {'type': 'cVariant', 'coerce': [_coerce.to_variant]}, 'outframe': {'type': 'cStr', 'coerce': _coerce.to_str}, 'veltype': {'type': 'cStr', 'coerce': _coerce.to_str}, 'restfreq': {'type': 'cVariant', 'coerce': [_coerce.to_variant]}, 'interpolation': {'type': 'cStr', 'coerce': _coerce.to_str, 'allowed': [ 'nearest', 'linear', 'cubic' ]}, 'perchanweightdensity': {'type': 'cBool'}, 'gridder': {'type': 'cStr', 'coerce': _coerce.to_str, 'allowed': [ 'widefield', 'wproject', 'awphpg', 'imagemosaic', 'standard', 'awproject', 'wprojectft', 'mosaicft', 'ft', 'ftmosaic', 'mosaic', 'awprojectft', 'gridft', 'awp2' ]}, 'facets': {'type': 'cInt'}, 'psfphasecenter': {'anyof': [{'type': 'cInt'}, {'type': 'cStr', 'coerce': _coerce.to_str}]}, 'wprojplanes': {'type': 'cInt'}, 'vptable': {'type': 'cStr', 'coerce': _coerce.to_str}, 'mosweight': {'type': 'cBool'}, 'aterm': {'type': 'cBool'}, 'psterm': {'type': 'cBool'}, 'wbawp': {'type': 'cBool'}, 'conjbeams': {'type': 'cBool'}, 'cfcache': {'type': 'cStr', 'coerce': _coerce.to_str}, 'usepointing': {'type': 'cBool'}, 'computepastep': {'type': 'cFloat', 'coerce': _coerce.to_float}, 'rotatepastep': {'type': 'cFloat', 'coerce': _coerce.to_float}, 'pointingoffsetsigdev': {'anyof': [{'type': 'cIntVec', 'coerce': [_coerce.to_list,_coerce.to_intvec]}, {'type': 'cFloatVec', 'coerce': [_coerce.to_list,_coerce.to_floatvec]}]}, 'pblimit': {'type': 'cFloat', 'coerce': _coerce.to_float}, 'normtype': {'type': 'cStr', 'coerce': _coerce.to_str}, 'deconvolver': {'type': 'cStr', 'coerce': _coerce.to_str, 'allowed': [ 'clarkstokes_exp', 'mtmfs', 'mem', 'clarkstokes', 'hogbom', 'clark_exp', 'clark', 'asp', 'multiscale' ]}, 'scales': {'anyof': [{'type': 'cIntVec', 'coerce': [_coerce.to_list,_coerce.to_intvec]}, {'type': 'cFloatVec', 'coerce': [_coerce.to_list,_coerce.to_floatvec]}]}, 'nterms': {'type': 'cInt'}, 'smallscalebias': {'type': 'cFloat', 'coerce': _coerce.to_float}, 'fusedthreshold': {'type': 'cFloat', 'coerce': _coerce.to_float}, 'largestscale': {'type': 'cInt'}, 'restoration': {'type': 'cBool'}, 'restoringbeam': {'anyof': [{'type': 'cStr', 'coerce': _coerce.to_str}, {'type': 'cStrVec', 'coerce': [_coerce.to_list,_coerce.to_strvec]}]}, 'pbcor': {'type': 'cBool'}, 'outlierfile': {'type': 'cStr', 'coerce': _coerce.to_str}, 'weighting': {'type': 'cStr', 'coerce': _coerce.to_str, 'allowed': [ 'briggsabs', 'briggs', 'briggsbwtaper', 'natural', 'radial', 'superuniform', 'uniform' ]}, 'robust': {'type': 'cFloat', 'coerce': _coerce.to_float, 'min': -2.0, 'max': 2.0}, 'noise': {'type': 'cVariant', 'coerce': [_coerce.to_variant]}, 'npixels': {'type': 'cInt'}, 'uvtaper': {'type': 'cStrVec', 'coerce': [_coerce.to_list,_coerce.to_strvec]}, 'niter': {'type': 'cInt'}, 'gain': {'type': 'cFloat', 'coerce': _coerce.to_float}, 'threshold': {'type': 'cVariant', 'coerce': [_coerce.to_variant]}, 'nsigma': {'type': 'cFloat', 'coerce': _coerce.to_float}, 'cycleniter': {'type': 'cInt'}, 'cyclefactor': {'type': 'cFloat', 'coerce': _coerce.to_float}, 'minpsffraction': {'type': 'cFloat', 'coerce': _coerce.to_float}, 'maxpsffraction': {'type': 'cFloat', 'coerce': _coerce.to_float}, 'nmajor': {'type': 'cInt'}, 'usemask': {'type': 'cStr', 'coerce': _coerce.to_str, 'allowed': [ 'user', 'pb', 'auto-multithresh' ]}, 'mask': {'anyof': [{'type': 'cStr', 'coerce': _coerce.to_str}, {'type': 'cStrVec', 'coerce': [_coerce.to_list,_coerce.to_strvec]}]}, 'pbmask': {'type': 'cFloat', 'coerce': _coerce.to_float}, 'sidelobethreshold': {'type': 'cFloat', 'coerce': _coerce.to_float}, 'noisethreshold': {'type': 'cFloat', 'coerce': _coerce.to_float}, 'lownoisethreshold': {'type': 'cFloat', 'coerce': _coerce.to_float}, 'negativethreshold': {'type': 'cFloat', 'coerce': _coerce.to_float}, 'smoothfactor': {'type': 'cFloat', 'coerce': _coerce.to_float}, 'minbeamfrac': {'type': 'cFloat', 'coerce': _coerce.to_float}, 'cutthreshold': {'type': 'cFloat', 'coerce': _coerce.to_float}, 'growiterations': {'type': 'cInt'}, 'dogrowprune': {'type': 'cBool'}, 'minpercentchange': {'type': 'cFloat', 'coerce': _coerce.to_float}, 'verbose': {'type': 'cBool'}, 'fastnoise': {'type': 'cBool'}, 'restart': {'type': 'cBool'}, 'savemodel': {'type': 'cStr', 'coerce': _coerce.to_str, 'allowed': [ 'none', 'virtual', 'modelcolumn' ]}, 'calcres': {'type': 'cBool'}, 'calcpsf': {'type': 'cBool'}, 'psfcutoff': {'type': 'cFloat', 'coerce': _coerce.to_float}, 'parallel': {'type': 'cBool'}, }
2163 doc = { 'vis': vis, 'selectdata': selectdata, 'field': field, 'spw': spw, 'timerange': timerange, 'uvrange': uvrange, 'antenna': antenna, 'scan': scan, 'observation': observation, 'intent': intent, 'datacolumn': datacolumn, 'imagename': imagename, 'imsize': imsize, 'cell': cell, 'phasecenter': phasecenter, 'stokes': stokes, 'projection': projection, 'startmodel': startmodel, 'specmode': specmode, 'reffreq': reffreq, 'nchan': nchan, 'start': start, 'width': width, 'outframe': outframe, 'veltype': veltype, 'restfreq': restfreq, 'interpolation': interpolation, 'perchanweightdensity': perchanweightdensity, 'gridder': gridder, 'facets': facets, 'psfphasecenter': psfphasecenter, 'wprojplanes': wprojplanes, 'vptable': vptable, 'mosweight': mosweight, 'aterm': aterm, 'psterm': psterm, 'wbawp': wbawp, 'conjbeams': conjbeams, 'cfcache': cfcache, 'usepointing': usepointing, 'computepastep': computepastep, 'rotatepastep': rotatepastep, 'pointingoffsetsigdev': pointingoffsetsigdev, 'pblimit': pblimit, 'normtype': normtype, 'deconvolver': deconvolver, 'scales': scales, 'nterms': nterms, 'smallscalebias': smallscalebias, 'fusedthreshold': fusedthreshold, 'largestscale': largestscale, 'restoration': restoration, 'restoringbeam': restoringbeam, 'pbcor': pbcor, 'outlierfile': outlierfile, 'weighting': weighting, 'robust': robust, 'noise': noise, 'npixels': npixels, 'uvtaper': uvtaper, 'niter': niter, 'gain': gain, 'threshold': threshold, 'nsigma': nsigma, 'cycleniter': cycleniter, 'cyclefactor': cyclefactor, 'minpsffraction': minpsffraction, 'maxpsffraction': maxpsffraction, 'nmajor': nmajor, 'usemask': usemask, 'mask': mask, 'pbmask': pbmask, 'sidelobethreshold': sidelobethreshold, 'noisethreshold': noisethreshold, 'lownoisethreshold': lownoisethreshold, 'negativethreshold': negativethreshold, 'smoothfactor': smoothfactor, 'minbeamfrac': minbeamfrac, 'cutthreshold': cutthreshold, 'growiterations': growiterations, 'dogrowprune': dogrowprune, 'minpercentchange': minpercentchange, 'verbose': verbose, 'fastnoise': fastnoise, 'restart': restart, 'savemodel': savemodel, 'calcres': calcres, 'calcpsf': calcpsf, 'psfcutoff': psfcutoff, 'parallel': parallel, }
2164 assert self._validate(doc, schema), create_error_string(_pc.errors)
2166 if len(imagename) == 0:
2167 raise Exception("imagename must be specified")
2169 self._vis = vis
2170 self._imagename = imagename
2171 self._selectdata = selectdata
2172 self._imsize = imsize
2173 self._cell = cell
2174 self._phasecenter = phasecenter
2175 self._projection = projection
2176 self._stokes = stokes
2177 self._startmodel = startmodel
2178 self._specmode = specmode
2179 self._reffreq = reffreq
2180 self._nchan = nchan
2181 self._start = start
2182 self._width = width
2183 self._outframe = outframe
2184 self._veltype = veltype
2185 self._restfreq = restfreq
2186 self._interpolation = interpolation
2187 self._perchanweightdensity = perchanweightdensity
2188 self._gridder = gridder
2189 self._facets = facets
2190 self._psfphasecenter = psfphasecenter
2191 self._wprojplanes = wprojplanes
2192 self._mosweight = mosweight
2193 self._aterm = aterm
2194 self._psterm = psterm
2195 self._wbawp = wbawp
2196 self._conjbeams = conjbeams
2197 self._usepointing = usepointing
2198 self._cfcache = cfcache
2199 self._pointingoffsetsigdev = pointingoffsetsigdev
2200 self._vpable = vptable
2201 self._computepastep = computepastep
2202 self._rotatepastep = rotatepastep
2203 self._pblimit = pblimit
2204 self._normtype = normtype
2205 self._deconvolver = deconvolver
2206 self._smallscalebias = smallscalebias
2207 self._fusedthreshold = fusedthreshold
2208 self._largestscale = largestscale
2209 self._niter = niter
2210 self._threshold = threshold
2211 self._cycleniter = cycleniter
2212 self._minpsffraction = minpsffraction
2213 self._maxpsffraction = maxpsffraction
2214 self._nsigma = nsigma
2215 self._nmajor = nmajor
2216 self._cyclefactor = cyclefactor
2217 self._scales = scales
2218 self._restoringbeam = restoringbeam
2219 self._pbcor = pbcor
2220 self._outlierfile = outlierfile
2221 self._nterms = nterms
2222 self._exe_cmds = [ ]
2223 self._exe_cmds_per_iter = [ ]
2224 self._history_filter = history_filter
2225 self._finalized = False
2226 self._field = field
2227 self._spw = spw
2228 self._timerange = timerange
2229 self._uvrange = uvrange
2230 self._antenna = antenna
2231 self._scan = scan
2232 self._observation = observation
2233 self._intent = intent
2234 self._datacolumn = datacolumn
2235 self._weighting = weighting
2236 self._robust = robust
2237 self._noise = noise
2238 self._uvtaper = uvtaper
2239 self._npixels = npixels
2240 self._gain = gain
2241 self._pbmask = pbmask
2242 self._sidelobethreshold = sidelobethreshold
2243 self._noisethreshold = noisethreshold
2244 self._lownoisethreshold = lownoisethreshold
2245 self._negativethreshold = negativethreshold
2246 self._smoothfactor = smoothfactor
2247 self._minbeamfrac = minbeamfrac
2248 self._cutthreshold = cutthreshold
2249 self._growiterations = growiterations
2250 self._dogrowprune = dogrowprune
2251 self._minpercentchange = minpercentchange
2252 self._verbose = verbose
2253 self._fastnoise = fastnoise
2254 self._savemodel = savemodel
2255 self._parallel = parallel
2256 self._usemask = usemask
2257 self._mask = mask
2258 self._restoration = restoration
2259 self._restart = restart
2260 self._calcres = calcres
2261 self._calcpsf = calcpsf
2262 self._psfcutoff = psfcutoff
2263 self.global_imdict = ImagingDict()
2264 self.current_imdict = ImagingDict()
2265 self._major_done = 0
2266 self.hasit = StopCodes(major=0,minor=0) # Convergence flag
2267 self._has_restored = False
2268 self._nfields = 1
2269 self._fieldnames = ['', ] # Default is a list of len 1
2271 self.stopdescription = '' # Convergence flag
2272 self._all_input_params = locals()
2273 del self._all_input_params['self']
2275 self._initial_mask_exists = []
2276 self._paramlist = None
2278 # Convert threshold from string to float, interpreting units.
2279 # XXX : We should ideally use quantities, but we are trying to
2280 # stick to "public API" funtions inside _gclean
2281 self._threshold_to_float()
2283 # Read in params and calculate nfields
2284 self._get_imager_parameters(locals())
2285 # Initialize convergence result for each field being imaged
2286 # this has been replaced by the ConvergenceResult class
2287 #self._initialize_convergence_result()
2289 self.hasit_per_field = [False for _ in range(self._nfields)]
2290 self.CR = ConvergenceResult(self._nfields, self._fieldnames)
2293 def _get_imager_parameters(self, inpparams):
2294 """
2295 Initialize ImagerParameters object and set the parameters
2296 """
2298 # Drop unnecessary params
2299 del inpparams['self']
2300 del inpparams['schema']
2301 del inpparams['history_filter']
2303 # deal with parameters that are not the same name
2304 inpparams, defparm = sanitize_tclean_inputs(inpparams)
2306 # Stolen from task_tclean.py
2307 bparm = {k: inpparams[k] if k in inpparams else defparm[k] for k in defparm.keys()}
2308 self._paramlist = ImagerParameters(**bparm)
2310 self._nfields = len(self._paramlist.allimpars.keys())
2311 self._fieldnames = [val['imagename'] for key, val in self._paramlist.allimpars.items()]
2313 return self._paramlist
2316 def __update_convergence(self, field=0):
2317 """
2318 Accumulates the per-channel/stokes summaryminor keys across all major cycle calls so far.
2320 The "iterDone" key will be replaced with "iterations", and for the "iterations" key,
2321 the value in the returned cummulative record will be a rolling sum of iterations done
2322 for tclean calls so far, one value per minor cycle.
2323 For example, if there have been two clean calls, and in the first call channel 0 had
2324 [1] iteration in 1 minor cycle, and for the second call channel 0 had [6, 10, 9, 1]
2325 iterations in 4 minor cycles), then the resultant "iterations" key for channel 0 would be:
2326 [1, 7, 17, 26, 27]
2328 Inputs:
2329 field: int, field index to update convergence for. Default is 0, which is the main field.
2330 If there are multiple fields, this will be the index of the field to update.
2331 Returns:
2332 outrec: dict, a nested dictionary with the following structure:
2333 {
2334 channel_id: {
2335 stokes_id: {
2336 summary_key: [values, one per minor cycle]
2337 }
2338 }
2339 }
2340 """
2342 keys = ['modelFlux', 'iterDone', 'peakRes', 'stopCode', 'cycleThresh']
2344 # Grab tuples of keys of interest
2345 outrec = {}
2346 for nn in range(self.global_imdict.nchan):
2347 outrec[nn] = {}
2348 for ss in range(self.global_imdict.nstokes):
2349 outrec[nn][ss] = {}
2350 for key in keys:
2351 # Replace iterDone with iterations
2352 if key == 'iterDone':
2353 # Maintain cumulative sum of iterations per entry
2354 outrec[nn][ss]['iterations'] = np.cumsum(self.global_imdict.get_key(key, field=field, stokes=ss, chan=nn))
2355 # Replace iterDone with iterations
2356 #outrec[nn][ss]['iterations'] = self.global_imdict.get_key(key, stokes=ss, chan=nn)
2357 else:
2358 outrec[nn][ss][key] = self.global_imdict.get_key(key, field=field, stokes=ss, chan=nn)
2360 return outrec
2364 def __add_per_major_items( self, tclean_ret, major_ret):
2365 '''Add meta-data about the whole major cycle, including 'cyclethreshold'
2366 '''
2367 # Original if code :
2368 #rdict = dict( major=dict( cyclethreshold=[tclean_ret['cyclethreshold']] if major_ret is None else (major_ret['cyclethreshold'] + [tclean_ret['cyclethreshold']]) ),
2369 # chan=chan_ret )
2370 # ---
2371 # original else code
2372 #rdict = dict( major=dict( cyclethreshold=major_ret['cyclethreshold'].append(tclean_ret['cyclethreshold']) ),
2373 # chan=chan_ret )
2375 rdict = {}
2377 for fidx, fname in enumerate(self._fieldnames):
2378 if 'cyclethreshold' in tclean_ret:
2379 if all([mj['major'] is None for mj in major_ret.values()]):
2380 rdict[fname] = dict(major=dict(cyclethreshold=[tclean_ret['cyclethreshold']]), chan=self.__update_convergence(field=fidx))
2381 else:
2382 rdict[fname] = dict(major=dict(cyclethreshold=major_ret[fname]['major']['cyclethreshold'] + [tclean_ret['cyclethreshold']]), chan=self.__update_convergence(field=fidx))
2383 else:
2384 # It really shouldn't get to this point...
2385 rdict = dict(major=dict([mj['major']['cyclethreshold'] for mj in major_ret.values()] + [tclean_ret['cyclethreshold']]), chan=self.__update_convergence(field=fidx))
2387 return rdict
2391 def _calc_deconv_controls(self, imdict, field = 0, niterleft=0, threshold=0, cycleniter=-1, hasit_per_field=None):
2392 """
2393 Calculate cycleniter and cyclethreshold for deconvolution.
2394 """
2396 use_cycleniter = niterleft #niter - imdict.returndict['iterdone']
2398 if cycleniter > -1 : # User is forcing this number
2399 use_cycleniter = min(cycleniter, use_cycleniter)
2401 psffrac = imdict.returndict['maxpsfsidelobe'] * self._cyclefactor
2402 psffrac = max(psffrac, self._minpsffraction)
2403 psffrac = min(psffrac, self._maxpsffraction)
2405 # Each field needs it's own cyclethreshold
2406 cyclethresh_per_field = []
2407 for fidx, fieldid in enumerate(self._paramlist.allimpars.keys()):
2408 if hasit_per_field is not None:
2409 if hasit_per_field[fidx]:
2410 # If convergence is hit, skip this field
2411 cyclethresh_per_field.append(None)
2412 continue
2414 if fidx not in imdict.returndict['summaryminor']:
2415 # If the field is not in the summaryminor, then we cannot calculate the cyclethreshold
2416 # So set to the default for now
2417 # This can happen if the field has not been imaged yet.
2418 cyclethresh_per_field.append(self._threshold)
2419 continue
2421 cyclethreshold = psffrac * imdict.get_peakres(field=fidx)
2422 cyclethreshold = max(cyclethreshold, threshold)
2423 cyclethresh_per_field.append(cyclethreshold)
2425 return int(use_cycleniter), cyclethresh_per_field
2428 def _check_initial_mask(self):
2429 """
2430 Check if a mask from a previous run exists on disk or not.
2431 """
2433 image_products = self.image_products()
2435 for product in image_products:
2436 if self._usemask == 'user' and self._mask == '':
2437 maskname = product['maskname']
2439 if os.path.exists(maskname):
2440 self._initial_mask_exists.append(True)
2441 else:
2442 self._initial_mask_exists.append(False)
2445 def _fix_initial_mask(self):
2446 """
2447 If on start up, no user mask is provided, then flip the initial mask to
2448 be all zeros for interactive use.
2449 """
2451 from casatools import image
2452 ia = image()
2454 image_products = self.image_products()
2456 if self._usemask == 'user' and self._mask == '':
2457 for idx, product in enumerate(image_products):
2458 maskname = product['maskname']
2460 # This means the mask was newly created by deconvolve, so flip it
2461 if os.path.exists(maskname) and self._initial_mask_exists[idx] is False:
2462 ia.open(maskname)
2463 ia.set(0.0)
2464 ia.close()
2467 def _update_peakres_masksum_multifield(self):
2468 """
2469 Update peak residual value in the return dict for each image on disk -
2470 main field and outliers
2471 """
2473 image_products = self.image_products()
2475 imstats = []
2477 for idx, product in enumerate(image_products):
2478 residname = product['residualname']
2479 maskname = product['maskname']
2481 statsdict = {}
2483 if not os.path.exists(maskname):
2484 maskname = ''
2486 peakres = imstat(imagename=residname, mask=f'''"{maskname}"''')['max']
2487 if len(maskname) > 0:
2488 masksum = imstat(imagename=maskname)['sum']
2489 else:
2490 masksum = []
2492 if len(peakres) > 0:
2493 peakres = peakres[0]
2494 else:
2495 peakres = None
2497 if len(masksum) > 0:
2498 masksum = masksum[0]
2499 else:
2500 masksum = None
2502 statsdict['idx'] = idx
2503 statsdict['imagename'] = product['imagename']
2504 statsdict['peakres'] = peakres
2505 statsdict['masksum'] = masksum
2507 imstats.append(statsdict)
2509 #return peakres, masksum
2510 return imstats
2513 def image_products(self, imageidx = None):
2514 """
2515 Function that returns the image names for all fields.
2517 Optionally if imageidx is specified, return the image names
2518 only for that specific index.
2520 By default it returns a list of dicts, where each element in the list is a
2521 field. The dicts have the following keys:
2523 - imageidx: index of the field (0 for the main field, 1 for the first outlier etc.)
2524 - imagename: name of the image
2525 - maskname: name of the mask
2526 - residualname: name of the residual
2527 - imagepath: path to the image
2528 """
2530 # Use pysynthimager to get the names of the outlierfields
2531 # Populate paths and names, generate list
2532 # Use this function instead of self._imagename wherever possible.
2534 allimpars = self._paramlist.allimpars
2535 imlist = []
2537 for key, value in allimpars.items():
2538 if imageidx is not None:
2539 if imageidx != int(key):
2540 continue
2542 imdict = {}
2544 fullpath = os.path.abspath(value['imagename'])
2545 imname = os.path.basename(fullpath)
2546 deconvolver = value['deconvolver']
2548 if deconvolver == 'mtmfs':
2549 suffix = '.tt0'
2550 else:
2551 suffix = ''
2553 imdict['name'] = imname
2554 imdict['imagename'] = imname + f'.image{suffix}'
2555 imdict['imagepath'] = os.path.dirname(fullpath)
2557 imdict['maskname'] = imname + '.mask'
2558 imdict['residualname'] = imname + f'.residual{suffix}'
2560 # Append to list
2561 imlist.append(imdict)
2563 return imlist
2565 def _generate_convergence_log(self, hasit_per_field, stopcode_min_per_field, stopdesc_per_field):
2566 """
2567 Generate a string that summarizes the convergence information for all fields, that
2568 can be passed in to the logger function.
2570 Inputs:
2571 hasit_per_field: list of bools, whether convergence was hit for each field
2572 stopcode_min_per_field: list of ints, the stop code for each field
2573 stopdesc_per_field: list of str, the stop description for each field
2575 Returns:
2576 logstr : str, a string that summarizes the convergence information for all fields
2577 """
2579 logstr = ''
2580 logstr += f"#################################### Convergence Log ###################################<br>"
2581 logstr += f"# {'Field Name':^20s} | {'Description':^46s} | Code (maj, min)<br>"
2582 logstr += "#----------------------------------------------------------------------------------------<br>"
2583 for fidx, fieldname in enumerate(self._fieldnames):
2584 if hasit_per_field[fidx]:
2585 logstr += f"# {fieldname:^20s} | {stopdesc_per_field[fidx]:^46s} | ({hasit_per_field[fidx]}, {stopcode_min_per_field[fidx]})<br>"
2586 else:
2587 logstr += f"# {fieldname:^20s} | {'Not converged':^46s} | ({hasit_per_field[fidx]}, {stopcode_min_per_field[fidx]})<br>"
2588 logstr += "########################################################################################<br>"
2589 logstr = logstr.replace(" ", " ") # Replace spaces with non-breaking spaces for HTML rendering
2590 return logstr
2593 def __next__( self ):
2594 """ Runs tclean and returns the (stopcode, convergence result) when executed with the python builtin next() function.
2596 The returned convergence result is a nested dictionary:
2597 {
2598 channel id: {
2599 stokes id: {
2600 summary key: [values, one per minor cycle]
2601 },
2602 },
2603 }
2605 See also: gclean.__update_convergence(...)
2606 """
2608 tclean_ret = {}
2609 deconv_ret = {}
2610 stopcode_min_per_field = [0 for _ in range(self._nfields)]
2611 stopdesc_per_field = ['' for _ in range(self._nfields)]
2613 # vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv------------>>> is done to produce an initial dirty image
2614 if self._niter < 1 and self.CR.convergence_result[2] is not None:
2615 self.CR.convergence_result = ( f'nothing to run, niter == {self._niter}',
2616 self.CR.convergence_result[1],
2617 self._major_done,
2618 self._nmajor,
2619 self._niter,
2620 self.CR.convergence_result[5] )
2621 return self.CR.convergence_result
2622 else:
2623 ### CALL SEQUENCE:
2624 ### tclean(niter=0),deconvolve(niter=0),tclean(niter=100),deconvolve(niter=0),tclean(niter=100),tclean(niter=0,restoration=True)
2625 self._exe_cmds_per_iter.append(0)
2626 try:
2627 if self.CR.convergence_result[1] is None:
2628 # initial call to tclean(...) creates the initial dirty image with niter=0
2630 # If calcres and calcpsf are False, no need to run initial tclean - assume image products already on disk.
2631 if not (self._calcres == False and self._calcpsf == False):
2632 casalog.post('Running initial major cycle to create first residual image.', 'INFO')
2633 print('Running initial major cycle to create first residual image.')
2634 tclean_ret = self._tclean( vis=self._vis, mask=self._mask, imagename=self._imagename, imsize=self._imsize, cell=self._cell, selectdata=self._selectdata, phasecenter=self._phasecenter, stokes=self._stokes,
2635 startmodel=self._startmodel, specmode=self._specmode, projection=self._projection, reffreq=self._reffreq, gridder=self._gridder, wprojplanes=self._wprojplanes, facets=self._facets,
2636 mosweight=self._mosweight, psfphasecenter = self._psfphasecenter, aterm=self._aterm, psterm=self._psterm, wbawp=self._wbawp, conjbeams=self._conjbeams, cfcache=self._cfcache,
2637 usepointing=self._usepointing, computepastep=self._computepastep, rotatepastep=self._rotatepastep, normtype=self._normtype, fusedthreshold=self._fusedthreshold,
2638 largestscale=self._largestscale, noise = self._noise, uvtaper=self._uvtaper, psfcutoff=self._psfcutoff, interpolation=self._interpolation,
2639 perchanweightdensity=self._perchanweightdensity, nchan=self._nchan, start=self._start, width=self._width, veltype=self._veltype, restfreq=self._restfreq, outframe=self._outframe,
2640 pointingoffsetsigdev=self._pointingoffsetsigdev, pblimit=self._pblimit, deconvolver=self._deconvolver, smallscalebias=self._smallscalebias, cyclefactor=self._cyclefactor,
2641 scales=self._scales, restoringbeam=self._restoringbeam, pbcor=self._pbcor, nterms=self._nterms, field=self._field, spw=self._spw, timerange=self._timerange, uvrange=self._uvrange,
2642 antenna=self._antenna, scan=self._scan, observation=self._observation, intent=self._intent, datacolumn=self._datacolumn, weighting=self._weighting, robust=self._robust,
2643 npixels=self._npixels, interactive=False, niter=0, gain=self._gain, calcres=self._calcres, calcpsf=self._calcpsf, restoration=False, parallel=self._parallel, fullsummary=True,
2644 outlierfile = self._outlierfile)
2646 # Check if a mask from a previous run exists on disk
2647 # This loops over all fields
2648 self._check_initial_mask()
2650 # deconv args are parsed correctly inside this function
2651 deconv_ret = self._deconv_over_fields(niter=0)
2652 # If no mask from a previous run exists, over-write the ones with zeros for the default mask
2653 self._fix_initial_mask()
2655 if len(tclean_ret) > 0 and len(deconv_ret) > 0:
2656 self.current_imdict.returndict = self.current_imdict.merge(tclean_ret, deconv_ret)
2657 elif len(tclean_ret) > 0 and len(deconv_ret) == 0:
2658 self.current_imdict.returndict = copy.deepcopy(tclean_ret)
2659 elif len(tclean_ret) == 0 and len(deconv_ret) > 0:
2660 self.current_imdict.returndict = copy.deepcopy(deconv_ret)
2661 else: # Both dicts are empty, this should never happen
2662 raise ValueError("Both tclean and deconvolve return dicts are empty. This should never happen.")
2665 self.global_imdict.returndict = self.current_imdict.returndict
2667 self.current_imdict.returndict['stopcode'] = self.hasit
2668 self.current_imdict.returndict['stopDescription'] = self.stopdescription
2669 self._major_done = 0
2670 else:
2671 # Reset convergence every time, since we return control to the GUI after a single major cycle
2672 self.current_imdict.returndict['iterdone'] = 0.
2674 # Mask can be updated here...
2675 # Check for mask update - peakres + masksum
2676 # imstats is a list of dicts
2677 imstats = self._update_peakres_masksum_multifield()
2679 self.hasit, self.stopdescription, self.hasit_per_field, __, __ = self.global_imdict.check_convergence(self._niter, self._threshold, self._nmajor, imstats=imstats)
2681 # Has not, i.e., not converged
2682 if not self.hasit.major:
2683 use_cycleniter, cyclethresh_per_field = self._calc_deconv_controls(self.current_imdict, niterleft=self._niter, threshold=self._threshold, cycleniter=self._cycleniter)
2685 # deconv args are parsed correctly inside this function
2686 # This will skip fields that have already converged
2687 deconv_ret = self._deconv_over_fields(niter=use_cycleniter, mask='', threshold=cyclethresh_per_field, hasit_per_field=self.hasit_per_field)
2689 # Joint major cycle over all fields
2690 # Run the major cycle
2691 tclean_ret = self._tclean( vis=self._vis, imagename=self._imagename, imsize=self._imsize, cell=self._cell,
2692 phasecenter=self._phasecenter, stokes=self._stokes, specmode=self._specmode, reffreq=self._reffreq,
2693 gridder=self._gridder, wprojplanes=self._wprojplanes, mosweight=self._mosweight, psterm=self._psterm,
2694 wbawp=self._wbawp, conjbeams=self._conjbeams, usepointing=self._usepointing, interpolation=self._interpolation,
2695 perchanweightdensity=self._perchanweightdensity, nchan=self._nchan, start=self._start,
2696 width=self._width, veltype=self._veltype, restfreq=self._restfreq, outframe=self._outframe,
2697 pointingoffsetsigdev=self._pointingoffsetsigdev, pblimit=self._pblimit, deconvolver=self._deconvolver,
2698 smallscalebias=self._smallscalebias, cyclefactor=self._cyclefactor, scales=self._scales,
2699 restoringbeam=self._restoringbeam, pbcor=self._pbcor, nterms=self._nterms, field=self._field,
2700 spw=self._spw, timerange=self._timerange, uvrange=self._uvrange, antenna=self._antenna,
2701 scan=self._scan, observation=self._observation, intent=self._intent, datacolumn=self._datacolumn,
2702 weighting=self._weighting, robust=self._robust, npixels=self._npixels, interactive=False,
2703 niter=0, restart=True, calcpsf=False, calcres=True, restoration=False, threshold=self._threshold,
2704 nsigma=self._nsigma, cycleniter=self._cycleniter, nmajor=1, gain=self._gain,
2705 sidelobethreshold=self._sidelobethreshold, noisethreshold=self._noisethreshold,
2706 lownoisethreshold=self._lownoisethreshold, negativethreshold=self._negativethreshold,
2707 minbeamfrac=self._minbeamfrac, growiterations=self._growiterations, dogrowprune=self._dogrowprune,
2708 minpercentchange=self._minpercentchange, fastnoise=self._fastnoise, savemodel='none',
2709 maxpsffraction=self._maxpsffraction,
2710 minpsffraction=self._minpsffraction, parallel=self._parallel, fullsummary=True, outlierfile=self._outlierfile)
2712 # Replace return dict with new return dict
2713 # The order of the dicts into merge is important.
2714 self.current_imdict.returndict = self.current_imdict.merge(tclean_ret, deconv_ret)
2716 # Append new return dict to global return dict
2717 self.global_imdict.returndict = self.global_imdict.concat(self.global_imdict.returndict, self.current_imdict.returndict)
2718 self._major_done = self.current_imdict.returndict['nmajordone']
2720 ## Decrement count for the major cycle just done...
2721 self.__decrement_counts()
2723 cycleniterleft = self._cycleniter - self.current_imdict.returndict['iterdone']
2725 # Use global imdict for convergence check
2726 if deconv_ret['stopcode'] == 7: ## Tell the convergence checker that the mask is zero and iterations were skipped
2727 self.hasit, self.stopdescription, self.hasit_per_field, stopcode_min_per_field, stopdesc_per_field = self.global_imdict.check_convergence(self._niter, self._threshold, self._nmajor, cycleniter=cycleniterleft)
2728 else:
2729 imstats = self._update_peakres_masksum_multifield()
2730 self.hasit, self.stopdescription, self.hasit_per_field, stopcode_min_per_field, stopdesc_per_field = self.global_imdict.check_convergence(self._niter, self._threshold, self._nmajor, cycleniter=cycleniterleft, imstats=imstats)
2732 # XXX : This is breaking, fix this to get the logging working correctly.
2733 logstr = self._generate_convergence_log(self.hasit_per_field, stopcode_min_per_field, stopdesc_per_field)
2734 self._exe_cmds.append(logstr)
2735 self._exe_cmds_per_iter[-1] += 1
2737 self.global_imdict.returndict['stopcode'] = self.hasit
2738 self.global_imdict.returndict['stopDescription'] = self.stopdescription
2741 if not self.hasit.major and self._usemask == 'auto-multithresh':
2742 # If we haven't converged, run deconvolve to update the mask
2743 # Note : This is only necessary if using auto-multithresh. If using the interactive viewer to draw, the mask gets updated as the regions
2744 # are added or removed in the interactive viewer.
2746 # deconv args are parsed correctly inside this function
2747 deconv_ret = self._deconv_over_fields(niter=0, mask='', hasit_per_field=self.hasit_per_field)
2749 if len(self.global_imdict.returndict) > 0 and 'summaryminor' in self.global_imdict.returndict and sum(map(len,self.global_imdict.returndict['summaryminor'].values())) > 0:
2750 self.CR.convergence_result = ( self.global_imdict.returndict['stopDescription'] if 'stopDescription' in self.global_imdict.returndict else '',
2751 self.global_imdict.returndict['stopcode'] if 'stopcode' in self.global_imdict.returndict else 0,
2752 self._major_done,
2753 self._nmajor,
2754 self._niter,
2755 self.__add_per_major_items( self.global_imdict.returndict,
2756 #self.CR.convergence_result[5]['major'],
2757 self.CR.convergence_result[5]))
2758 #self.__update_convergence( ) ) )
2759 else:
2760 self.CR.convergence_result = ( f'tclean returned an empty result',
2761 self.CR.convergence_result[1],
2762 self._major_done,
2763 self._nmajor,
2764 self._niter,
2765 self.CR.convergence_result[5] )
2766 except Exception as e:
2767 self.CR.convergence_result = ( str(e),
2768 -1,
2769 self._major_done,
2770 self._nmajor,
2771 self._niter,
2772 self.CR.convergence_result[5] )
2773 return self.CR.convergence_result
2775 return self.CR.convergence_result
2777 def __decrement_counts( self ):
2778 ## Update niterleft and nmajorleft now.
2779 if not any (self.hasit): ##If not yet converged.
2780 if self._nmajor != -1: ## If -1, don't touch it.
2781 self._nmajor = self._nmajor - 1
2782 if self._nmajor < 0: ## Force a floor
2783 self._nmajor = 0
2784 self._niter = self._niter - self.current_imdict.get_key('iterdone')
2785 if self._niter < 0: ## This can happen when we're counting niter across channels in a single minor cycle set, and it crosses the total.
2786 self._niter = 0 ## Force a floor
2787 else:
2788 return ##If convergence has been reached, don't try to decrement further.
2791 def __reflect_stop( self ):
2792 ## if python wasn't hacky, you would be able to try/except/raise in lambda
2793 #time.sleep(1)
2794 try:
2795 return self.__next__( )
2796 except StopIteration:
2797 raise StopAsyncIteration
2799 async def __anext__( self ):
2800 ### asyncio.run cannot be used here because this is called
2801 ### within an asyncio loop...
2802 loop = asyncio.get_event_loop( )
2803 result = await loop.run_in_executor( None, self.__reflect_stop )
2804 return result
2806 def __iter__( self ):
2807 return self
2809 def __aiter__( self ):
2810 return self
2812 def __split_filename( self, path ):
2813 return os.path.splitext(os.path.basename(path))
2815 def __default_mask_name( self ):
2816 imgparts = self.__split_filename( self._imagename )
2817 return f'{imgparts[0]}.mask'
2819 def __del__(self):
2820 pass
2822 def mask(self):
2823 #return self.__default_mask_name() if self._mask == '' else self._mask
2824 return f'{self._imagename}.mask' if self._mask == '' else self._mask
2826 def reset(self):
2827 #if not self._finalized:
2828 # raise RuntimeError('attempt to reset a gclean run that has not been finalized')
2829 self._finalized = False
2830 self.CR.convergence_result = ( None,
2831 self.CR.convergence_result[1],
2832 self._major_done,
2833 self._nmajor,
2834 self._niter,
2835 self.CR.convergence_result[5] )
2838 def restore(self):
2839 """
2840 Restores the final image, and returns a path to the restored image.
2841 If savemodel is specified, first predict the model column with tclean,
2842 then run restore with deconvolve.
2843 """
2845 casalog.post('Restoring final image (iterating over all fields)', 'INFO')
2846 print("Restoring the final image (iterating over all fields)")
2848 if self._savemodel != 'none':
2849 # Run the major cycle to predict the model column
2850 tclean_ret = self._tclean( vis=self._vis, imagename=self._imagename, imsize=self._imsize, cell=self._cell,
2851 phasecenter=self._phasecenter, stokes=self._stokes, specmode=self._specmode, reffreq=self._reffreq,
2852 gridder=self._gridder, wprojplanes=self._wprojplanes, mosweight=self._mosweight, psterm=self._psterm,
2853 wbawp=self._wbawp, conjbeams=self._conjbeams, usepointing=self._usepointing, interpolation=self._interpolation,
2854 perchanweightdensity=self._perchanweightdensity, nchan=self._nchan, start=self._start,
2855 width=self._width, veltype=self._veltype, restfreq=self._restfreq, outframe=self._outframe,
2856 pointingoffsetsigdev=self._pointingoffsetsigdev, pblimit=self._pblimit, deconvolver=self._deconvolver,
2857 smallscalebias=self._smallscalebias, cyclefactor=self._cyclefactor, scales=self._scales,
2858 restoringbeam=self._restoringbeam, pbcor=self._pbcor, nterms=self._nterms, field=self._field,
2859 spw=self._spw, timerange=self._timerange, uvrange=self._uvrange, antenna=self._antenna,
2860 scan=self._scan, observation=self._observation, intent=self._intent, datacolumn=self._datacolumn,
2861 weighting=self._weighting, robust=self._robust, npixels=self._npixels, interactive=False,
2862 niter=0, restart=True, calcpsf=False, calcres=False, restoration=False, threshold=self._threshold,
2863 nsigma=self._nsigma, cycleniter=self._cycleniter, nmajor=1, gain=self._gain,
2864 sidelobethreshold=self._sidelobethreshold, noisethreshold=self._noisethreshold,
2865 lownoisethreshold=self._lownoisethreshold, negativethreshold=self._negativethreshold,
2866 minbeamfrac=self._minbeamfrac, growiterations=self._growiterations, dogrowprune=self._dogrowprune,
2867 minpercentchange=self._minpercentchange, fastnoise=self._fastnoise, savemodel=self._savemodel,
2868 maxpsffraction=self._maxpsffraction, outlierfile=self._outlierfile,
2869 minpsffraction=self._minpsffraction, parallel=self._parallel, fullsummary=True )
2872 deconv_ret = self._deconv_over_fields(niter=0, mask='', restoration=True)
2874 self._has_restored = True
2876 # Write history into relevant images for all fields
2877 for fidx, field in enumerate(self._fieldnames):
2878 #params = cleanhelper.get_func_params(self.__init__, self._all_input_params)
2879 cleanhelper.write_tclean_history(self._paramlist.allimpars[f"{fidx}"]['imagename'], "InteractiveClean", tuple(self._all_input_params['doc'].items()), casalog)
2881 return self.global_imdict.returndict
2883 # return { "image": f"{self._imagename}.image" }
2885 def has_next(self):
2886 return not self._finalized