scriptForImaging_template.py

#>>> =====================================================================================#
#>>>                        TEMPLATE IMAGING SCRIPT                                       #
#>>> =====================================================================================#
#>>>
#>>> Updated: Wed  11:49 EDT 2021

#>>> Lines beginning with '#>>>' are instructions to the data imager
#>>> and will be removed from the script delivered to the PI. If you
#>>> would like to include a comment that will be passed to the PI,
#>>> begin the line with a single '#', i.e., standard python comment
#>>> syntax.  Helpful tip: Use the commands %cpaste or %paste to copy
#>>> and paste indented sections of code into the casa command line.

#>>> The commands below serve as a guide to best practices for imaging
#>>> ALMA data. It does not replace careful thought on your part while
#>>> imaging the data. You can remove or modify sections as necessary
#>>> depending on your particular imaging case (e.g., no
#>>> self-calibration, continuum only.) Please read the NA Imaging
#>>> Guide
#>>> (https://staff.nrao.edu/wiki/bin/view/NAASC/NAImagingScripts) for
#>>> more information.  Before imaging, you should use the commands
#>>> the first section of this script to prep the data for imaging.
#>>> The commands in both sections should be able to be run as as
#>>> standard Python script. However, the cleaning in this script is
#>>> done interactively making the final product somewhat dependent on
#>>> the individual doing the clean -- please clean conservatively
#>>> (i.e., don't box every possible source). The final data products
#>>> are the primary beam corrected images (*.pbcor), and the primary
#>>> beams (*.pb). These images should be converted to fits at the end
#>>> of the script (see example at the end of this file).  This script
#>>> (and the associated guide) are under active development. Please
#>>> contact Amanda Kepley (akepley@nrao.edu) if you have any
#>>> suggested changes or find any bugs that are almost certainly
#>>> there.

# This script has been tested for CASA 6.2.1.7.


########################################
# Check CASA version

import re

try:
    import casalith
except:
    print("Script requires CASA 6.0 or greater")


if casalith.compare_version("<",[6,2,1,7]):
    print("Please use CASA version greater than or equal to 6.2.1.7 with this script")


##################################################
# Create an Averaged Continuum MS

#>>> Continuum images can be sped up considerably by averaging the data
#>>> together to reduce overall volume. Since the sensitivity of a
#>>> continuum image depends on its bandwidth, continuum images are
#>>> typically made by including as much bandwidth as possible in the
#>>> data while excluding any line emission. The following plotms command
#>>> pages through the spectral windows in a project allowing you to
#>>> identify channel ranges within spectral windows that do not include
#>>> *strong* line emission. You will form a continuum image by averaging
#>>> the line-free spws and/or channel ranges within spws. In most cases,
#>>> you will not need to create an image to select line channels,
#>>> although you can suggest this to the PI as a possible path for
#>>> future exploration in the README file for cases where there is
#>>> wide-spread line emission.
#>>>
#>>> For a project with continuum target sensitivities, it is worth
#>>> checking the OT to see what continuum bandwidth the PI was
#>>> anticipating. In many cases, the continuum-only windows will be
#>>> specified in the OT, in general these have the broadest bandwidths
#>>> (~2GHz) with a small number of channels (128).  However, other
#>>> windows may be combined with these designated continuum windows to
#>>> increase the continuum sensitivity. In general, it is not necessary
#>>> to include narrow spectral windows (<250MHz) in the continuum image.

finalvis='calibrated_final.ms' # This is your output ms from the data
                               # preparation script.

# Use plotms to identify line and continuum spectral windows.
#>>> If you have a project with multiple fields, you will want to run
#>>> the following plotms command separately for each source. If the
#>>> spectra for each field are significantly different from each other,
#>>> it may be necessary to make separate average continuum  and
#>>> continuum-subtracted measurement sets for each field.
plotms(vis=finalvis, xaxis='channel', yaxis='amplitude',
       ydatacolumn='data',
       avgtime='1e8', avgscan=True, avgchannel='1', 
       iteraxis='spw' )


#>>> Note that when you average channels in plotms, it displays
#>>> the "bin" number rather than the average channel number of each
#>>> bin. 

#>>> If you don't see any obvious lines in the above plot, you may to try
#>>> to set avgbaseline=True with uvrange (e.g., <100m). Limiting the
#>>> uvrange to the short baselines greatly improves the visibility of
#>>> lines with extended emission.

#>>> If you have multiple sources, the line channel ranges may be
#>>> different for different sources. Thus you would need to repeat the
#>>> process below for each source.


# Set spws to be used to form continuum
contspws = '0,1,2,3'

# If you have complex line emission and no dedicated continuum
# windows, you will need to flag the line channels prior to averaging.
flagmanager(vis=finalvis,mode='save',
            versionname='before_cont_flags')

initweights(vis=finalvis,wtmode='weight',dowtsp=True)

# Flag the "line channels"
flagchannels='2:1201~2199,3:1201~2199' # In this example , spws 2&3 have a line between channels 1201 and 2199 and spectral windows 0 and 1 are line-free.


flagdata(vis=finalvis,mode='manual',
          spw=flagchannels,flagbackup=False)

# check that flags are as expected, NOTE must check reload on plotms
# gui if its still open.
plotms(vis=finalvis,yaxis='amp',xaxis='channel',
       avgchannel='1',avgtime='1e8',avgscan=True,iteraxis='spw') 

# Average the channels within spws
contvis='calibrated_final_cont.ms'
rmtables(contvis)
os.system('rm -rf ' + contvis + '.flagversions')

#>>> Note that to mitigate bandwidth smearing, please keep the width
#>>> of averaged channels less than 125MHz in Band 3, 4, and 6, and 250MHz
#>>> in Band 7 for both TDM and FDM modes. For example, for a 2GHz TDM window
#>>> with 15.625 MHz channels, this means that the maximum width parameter
#>>> should be 8 channels for Bands 3, 4, and 6 and 16 channels for Band 7.
#>>> This is especially important for any long baseline data. These limits
#>>> have been designed to have minimize the reduction of the peak flux to
#>>> 95%. See the "for continuum" header for more information on the imaging
#>>> wiki for more infomration.

#>>> Note that in CASA 5.1, split2 is now split. Previously split2 was
#>>> needed to deal correctly with channelized weights.
split(vis=finalvis,
     spw=contspws,      
     outputvis=contvis,
      width=[256,8,8,8], # number of channels to average together. The final channel width should be less than 125MHz in Bands 3, 4, and 6 and 250MHz in Band 7.
     datacolumn='data')


# Check the weights. You will need to change antenna and field to
# appropriate values
plotms(vis=contvis, yaxis='wtsp',xaxis='freq',spw='',antenna='DA42',field='0')

# If you flagged any line channels, restore the previous flags
flagmanager(vis=finalvis,mode='restore',
            versionname='before_cont_flags')

# Inspect continuum for any problems
plotms(vis=contvis,xaxis='uvdist',yaxis='amp',coloraxis='spw')

# #############################################
# Image Parameters

#>>> You're now ready to image. Review the science goals in the OT and
#>>> set the relevant imaging parameters below. 

# source parameters
# ------------------

field='0' # science field(s). For a mosaic, select all mosaic fields. DO NOT LEAVE BLANK ('') OR YOU WILL POTENTIALLY TRIGGER A BUG IN CLEAN THAT WILL PUT THE WRONG COORDINATE SYSTEM ON YOUR FINAL IMAGE.
#gridder='standard' # uncomment if single field 
# gridder='mosaic' # uncomment if mosaic or if combining one 7m and one 12m pointing.
# phasecenter=3 # uncomment and set to field number for phase
                # center. Note lack of ''.  Use the weblog to
                # determine which pointing to use. Remember that the
                # field ids for each pointing will be re-numbered
                # after your initial split. You can also specify the
                # phase center using coordinates, e.g.,
                # phasecenter='J2000 19h30m00 -40d00m00'.
# phasecenter = 'TRACKFIELD' # If imaging an ephemeris object (planet, etc), the phasecenter needs to be TRACKFIELD, not a field number as above.


# image parameters.
# ----------------

#>>> Generally, you want 5-8 cells (i.e., pixels) across the narrowest
#>>> part of the beam. You can estimate the beam size using the following
#>>> equation: 206265.0/(longest baseline in wavelengths).  To determine
#>>> the longest baseline, use plotms with xaxis='uvwave' and
#>>> yaxis='amp'. Divide the estimated beam size by five to eight to get
#>>> your cell size. It's better to error on the side of slightly too
#>>> many cells per beam than too few. Once you have made an image,
#>>> please re-assess the cell size based on the beam of the image. You can
#>>> check your cell size using  au.pickCellSize('calibrated_final.ms'). Note
#>>> however, that this routine does not take into account the projection of
#>>> the baseline onto the sky, so the plotms method described above is overall
#>>> more accurate.

#>>> To determine the image size (i.e., the imsize parameter), first you
#>>> need to figure out whether the ms is a mosaic by either looking out
#>>> the output from listobs/vishead or checking the spatial setup in the OT. For
#>>> single fields, an imsize equal to the size of the primary beam is
#>>> usually sufficient. The ALMA 12m primary beam in arcsec scales as
#>>> 6300 / nu[GHz] and the ALMA 7m primary beam in arcsec scales as
#>>> 10608 / nu[GHz], where nu[GHz] is the sky frequency. However, if
#>>> there is significant point source and/or extended emission beyond
#>>> the edges of your initial images, you should increase the imsize to
#>>> incorporate more emission. For mosaics, you can get the imsize from
#>>> the spatial tab of the OT. The parameters "p length" and "q length"
#>>> specify the dimensions of the mosaic. If you're imaging a mosaic,
#>>> pad the imsize substantially to avoid artifacts. Note that the imsize
#>>> parameter is in PIXELS, not arcsec, so you will need to divide the image size
#>>> in arcsec by the pixel size to determine a value for imsize.

#>>> Note that you can check your image size using 
#>>> au.pickCellSize('calibrated_final.ms', imsize=True). This task
#>>> now works both mosaics and single fields, but has not been tested
#>>> extensively on mosaics. Please report any large issues to
#>>> Todd Hunter. Note that au.pickCellSize does not take
#>>> into account the projection of the baselines, so the plotms
#>>> method is more accurate.

cell='1arcsec' # cell size for imaging.
imsize = [128,128] # size of image in pixels.

# velocity parameters
# -------------------

outframe='lsrk' # velocity reference frame. 
veltype='radio' # velocity type. 

#>>> Note on veltype: For quality assurance purposes, we recommend keeping veltype
#>>> set to radio, regardless of the velocity frame listed the object in the OT.
#>>> If the sensitivity in the OT is defined using a velocity width, then the 'radio'
#>>> definition of the velocity frame is used regardless of the velocity
#>>> definition in the "source parameters" tab of the OT.

# imaging control
# ----------------

# The cleaning below is done interactively, so niter and threshold can
# be controlled within clean. 

weighting = 'briggs'
robust=0.5
niter=1000
threshold = '0.0mJy'

#>>> Guidelines for setting robust:

#>>> Robust < 0.0 is not recommended for mosaics with poor-uv
#>>> coverage. Using values of robust less than or equal to 0.0 will
#>>> lead to major artifacts in the images including uneven noise
#>>> across the image.

#>>> If you are uv-tapering the data, you should set robust=2 (natural
#>>> weighting) to avoid upweighting points that are going to be
#>>> downweighted by uv-taper.

#############################################
# Imaging the Continuuum

# Set the ms and continuum image name.
contvis = 'calibrated_final_cont.ms'         

#>>> Generate the relevant image name and copy and paste imagename for PI

#>>> The spwmap parameter renumbers the images to match the archival
#>>> spectral windows, which aren't renumbered, unlike the manual
#>>> imaging spectral windows, which are renumbered by split in the
#>>> imaging prep script. It is defined similarly to other spwmap
#>>> parameters. For example, to map spw (e.g., 0,1,2,3) in
#>>> *.split.cal to their original values (e.g., 17, 19, 21, 23), you
#>>> would set spwmap=[17,19,21,23]. The sciencespws variable was
#>>> created in the scriptForImagingPrep.py in the first step. You can
#>>> also do a listobs on the original ms to get the spectral windows,
#>>> e.g.,
#>>> listobs(vis='uid___A002_Xc3412f_X53ff.ms',intent='OBSERVE_TARGET*',spw='*FULL_RES*')

#>>>  aU.genImageName(vis=contvis,spw=list(map(int,contspws.split(','))),field=int(field.split('~')[0]),imtype='mfs',targettype='sci',stokes='I',mous='',modtext='manual',spwmap=list(map(int,sciencespws.split(','))))
contimagename = '' 

# If necessary, run the following commands to get rid of older clean
# data.

#clearcal(vis=contvis)
#delmod(vis=contvis)

for ext in ['.image','.mask','.model','.image.pbcor','.psf','.residual','.pb','.sumwt','.weight']:
    rmtables(contimagename+ext)

#>>> If you're going be be imaging with nterms>1, then you also need
#>>> to removed the *.tt0, and *.tt1 images in additional to those
#>>> listed above.

#>>> If the fractional bandwidth for the aggregate continuum is
#>>> greater than 10%, set deconvolver='mtmfs' to use multi-term,
#>>> multi-frequency synthesis. This algorithm takes into account the
#>>> spatial spectral index variations in an image.  Note that only
#>>> ALMA Band 3 and the lower end of Band 4 can have fractional
#>>> bandwidths of greater than 10% and only when both sidebands are
#>>> employed.

tclean(vis=contvis,
       imagename=contimagename,
       field=field,
       #  phasecenter=phasecenter, # uncomment if mosaic or imaging an ephemeris object     
       # mosweight = True, # uncomment if mosaic
       specmode='mfs',
       deconvolver='hogbom', 
       # Uncomment the below to image with nterms>1 when the fractional bandwidth is greater than 10%.
       #deconvolver='mtmfs',
       #nterms=2,
       imsize = imsize, 
       cell= cell, 
       weighting = weighting,
       robust = robust,
       niter = niter, 
       threshold = threshold,
       interactive = True,
       gridder = gridder,
       pbcor = True,
       usepointing=False)
       
#>>> If interactively cleaning (interactive=True), then note number of
#>>> iterations at which you stop for the PI. This number will help
#>>> the PI replicate the delivered images. Do not clean empty
#>>> images. Just click the red X to stop the interactive and note the
#>>> RMS.

#>>> Note RMS for PI. 

# If you'd like to redo your clean, but don't want to make a new mask
# use the following commands to save your original mask. This is an optional step.
#contmaskname = 'cont.mask'
##rmtables(contmaskname) # if you want to delete the old mask
#os.system('cp -ir ' + contimagename + '.mask ' + contmaskname)

##############################################
# Self-calibration on the continuum [OPTIONAL]

#>>> Self-calibration solutions can be determined when a source exhibits
#>>> a strong continuum emission, preferably near the phase
#>>> center. Self-calibration can be attempted in cases where the
#>>> expected rms is not reached on the continuum (for projects where the
#>>> sensitivity is defined for the continuum) or the line data (for
#>>> projects where the sensitivity is defined for the line data). It
#>>> should not be attempted on 'partial' datasets (ACA and _TC datasets),
#>>> and, in the interest of reducing the time required to image, is not
#>>> recommended for datasets which do meet the rms requirement

#>>> The example here obtains solutions from the scan time to
#>>> down to times as short as per integration. Depending on the source,
#>>> you may not be able to find solution on timescales that short and
#>>> may need to adjust the solint parameter.

contvis = 'calibrated_final_cont.ms'         
contimagename = '' # Grab from continuum imaging step above if needed.

refant = 'DV09' # reference antenna.

#>>> Choose a reference antenna that's in the array. The tasks plotants
#>>> and listobs/vishead can tell you what antennas are in the array. For
#>>> data sets with multiple executions, you will want to choose an antenna
#>>> that's present in all the executions. The task au.commonAntennas()
#>>> can help with this.

#>>> Indicate the spectral window mapping below. The spwmap map
#>>> variable is a list that consists of n entries where n is the
#>>> number of spectral windows. The value of each entry tells you
#>>> which solution to apply to the indicated spectral window. For a
#>>> single execution, you would apply the calibration from spw 0 to
#>>> all spws since you are combining all spws in the solution below,
#>>> so spwmap=[0,0,0,0]. For multiple executions, you want to apply
#>>> the solution for each execution to itself. For example, if you had 4
#>>> spectral windows in a data set and two executions resulting in 8
#>>> total spectral windows, you would want to apply the solutions for
#>>> the first execution to itself and the solution to the second
#>>> execution to itself, so spwmap=[0,0,0,0,4,4,4,4]

spwmap = [0,0,0,0] # mapping self-calibration solutions to individual spectral windows.

# save initial flags in case you don't like the final
# self-calibration. The task applycal will flag data that doesn't have
# solutions.
flagmanager(vis=contvis,mode='save',versionname='before_selfcal',merge='replace')

# Get rid of any models that might be hanging around in the image header
delmod(vis=contvis,otf=True,scr=True)

# If you are re-doing your self-cal, uncomment the next line to reset
# your corrected data column back to its original state and get rid of
# the old model. You can check the contents of the model and corrected
# data columns by plotting them using plotms. For example, 
# plotms(vis=contvis, xaxis='uvwave', yaxis='amplitude', ydatacolumn='model',field=field)

# clearcal(vis=contvis)
# delmod(vis=contvis,otf=True,scr=True)

# shallow clean on the continuum

for ext in ['.image','.mask','.model','.image.pbcor','.psf','.residual','.pb','.sumwt','.weight']:
    rmtables(contimagename + '_p0'+ ext)

tclean(vis=contvis,
       imagename=contimagename + '_p0',
       field=field,
       #phasecenter=phasecenter, # uncomment if mosaic or imaging an ephemeris object
       # mosweight = True, # uncomment if mosaic
       specmode='mfs',
       deconvolver='hogbom',
       # Uncomment the below to image with nterms>1.
       #deconvolver='mtmfs',
       #nterms=2,
       imsize = imsize, 
       cell= cell, 
       weighting = weighting, 
       robust=robust,
       niter=niter, 
       threshold=threshold, 
       interactive=True,
       gridder=gridder,
       savemodel='modelcolumn',
       usepointing=False)

#>>> Note number of iterations performed.

#>>> Note: before proceeding, you should verify that TCLEAN saved the MODEL 
#>>> column. If it did, you will see a message like
#>>> INFO .... ------ Predict Model ------
#>>> INFO ... Saving model column
#>>> if you don't see this in the CASA logger, you should call TCLEAN again 
#>>> with niter=0, calcpsf=False, calcres=False to populate the modelcolumn

# per scan solution
rmtables('pcal1')
gaincal(vis=contvis,
        caltable='pcal1',
        field=field,
        gaintype='T',
        refant=refant, 
        calmode='p',
        combine='spw', 
        solint='inf',
        minsnr=3.0,
        minblperant=6)

#>>> If many of the above many solutions are flagged, consider setting
#>>> minsnr=1.5 and comparing the solutions. For low (<~500) dynamic
#>>> range cases, including a bit more random noise in the solution
#>>> has only a small effect on the image.

# Check the solution
plotms(vis='pcal1',
       xaxis='time',
       yaxis='phase',
       iteraxis='antenna',
       plotrange=[0,0,-180,180])

# apply the calibration to the data for next round of imaging
applycal(vis=contvis,
         field=field,
         spwmap=spwmap, 
         gaintable=['pcal1'],
         gainfield='',
         calwt=False, 
         flagbackup=False,
         interp='linearperobs')

# Save the flags in case you need to go back to this step. 
flagmanager(vis=contvis,mode='save',versionname='after_pcal1')

# clean deeper
for ext in ['.image','.mask','.model','.image.pbcor','.psf','.residual','.pb','.sumwt','.weight']:
    rmtables(contimagename + '_p1'+ ext)

tclean(vis=contvis,
       imagename=contimagename + '_p1',
       field=field,
       # phasecenter=phasecenter, # uncomment if mosaic or imaging an ephemeris object
       # mosweight = True, # uncomment if mosaic
       specmode='mfs',
       deconvolver='hogbom',
       # Uncomment the below to image with nterms>1.
       #deconvolver='mtmfs',
       #nterms=2,
       imsize = imsize, 
       cell= cell, 
       weighting = weighting, 
       robust=robust,
       niter=niter, 
       threshold=threshold, 
       interactive=True,
       gridder=gridder,
       #pbcor = True, #if final image
       savemodel='modelcolumn',
       usepointing=False)

#>>> Note number of iterations performed.


#>>> Note: before proceeding, you should verify that TCLEAN saved the MODEL 
#>>> column. If it did, you will see a message like
#>>> INFO .... ------ Predict Model ------
#>>> INFO ... Saving model column
#>>> if you don't see this in the CASA logger, you should call TCLEAN again 
#>>> with niter=0, calcpsf=False, calcres=False to populate the modelcolumn

# shorter solution
rmtables('pcal2')
gaincal(vis=contvis,
        field=field,
        caltable='pcal2',
        gaintype='T',
        refant=refant, 
        calmode='p',
        combine='spw', 
        solint='30.25s', # solint=30.25s gets you five 12m integrations, while solint=50.5s gets you five 7m integration
        minsnr=3.0,
        minblperant=6)

# Check the solution
plotms(vis='pcal2',
       xaxis='time',
       yaxis='phase',
       iteraxis='antenna',
       plotrange=[0,0,-180,180])

# apply the calibration to the data for next round of imaging
applycal(vis=contvis,
         spwmap=spwmap, 
         field=field,
         gaintable=['pcal2'],
         gainfield='',
         calwt=False, 
         flagbackup=False,
         interp='linearperobs')

# Save the flags in case you need to go back to this step. 
flagmanager(vis=contvis,mode='save',versionname='after_pcal2')

# clean deeper
for ext in ['.image','.mask','.model','.image.pbcor','.psf','.residual','.pb','.sumwt','.weight']:
    rmtables(contimagename + '_p2'+ ext)

tclean(vis=contvis,
       imagename=contimagename + '_p2',
       field=field,
       # phasecenter=phasecenter, # uncomment if mosaic or imaging an ephemeris object
       # mosweight = True, # uncomment if mosaic
       specmode='mfs',
       deconvolver='hogbom',
       # Uncomment the below to image with nterms>1.
       #deconvolver='mtmfs',
       #nterms=2,
       imsize = imsize, 
       cell= cell, 
       weighting = weighting, 
       robust=robust,
       niter=niter, 
       threshold=threshold, 
       interactive=True,
       gridder=gridder,
       #pbcor = True, #if final image
       savemodel='modelcolumn',
       usepointing=False)

#>>> Note number of iterations performed.

#>>> Note: before proceeding, you should verify that TCLEAN saved the MODEL 
#>>> column. If it did, you will see a message like
#>>> INFO .... ------ Predict Model ------
#>>> INFO ... Saving model column
#>>> if you don't see this in the CASA logger, you should call TCLEAN again 
#>>> with niter=0, calcpsf=False, calcres=False to populate the modelcolumn

# shorter solution
rmtables('pcal3')
gaincal(vis=contvis,
        field=field,
        caltable='pcal3',
        gaintype='T',
        refant=refant, 
        calmode='p',
        combine='spw', 
        solint='int',
        minsnr=3.0,
        minblperant=6)

# Check the solution
plotms(vis='pcal3',
       xaxis='time',
       yaxis='phase',
       iteraxis='antenna',
       plotrange=[0,0,-180,180])

# apply the calibration to the data for next round of imaging
applycal(vis=contvis,
         spwmap=spwmap,
         field=field,
         gaintable=['pcal3'],
         gainfield='',
         calwt=False, 
         flagbackup=False,
         interp='linearperobs')

# Save the flags in case you need to go back to this step. 
flagmanager(vis=contvis,mode='save',versionname='after_pcal3')

# do the amplitude self-calibration.
for ext in ['.image','.mask','.model','.image.pbcor','.psf','.residual','.pb','.sumwt','.weight']:
    rmtables(contimagename + '_p3'+ ext)

tclean(vis=contvis,
       imagename=contimagename + '_p3',
       field=field,
       # phasecenter=phasecenter, # uncomment if mosaic or imaging an ephemeris object
       # mosweight = True, # uncomment if mosaic
       specmode='mfs',
       deconvolver='hogbom',
       # Uncomment the below to image with nterms>1.
       #deconvolver='mtmfs',
       #nterms=2,
       imsize = imsize, 
       cell= cell, 
       weighting = weighting, 
       robust=robust,
       niter=niter, 
       threshold=threshold, 
       interactive=True,
       gridder=gridder,
       #pbcor = True, #if final image
       savemodel='modelcolumn',
       usepointing=False)


#>>> Note number of iterations performed.

#>>> Note: before proceeding, you should verify that TCLEAN saved the MODEL 
#>>> column. If it did, you will see a message like
#>>> INFO .... ------ Predict Model ------
#>>> INFO ... Saving model column
#>>> if you don't see this in the CASA logger, you should call TCLEAN again 
#>>> with niter=0, calcpsf=False, calcres=False to populate the modelcolumn


rmtables('apcal')
gaincal(vis=contvis,
        field=field,
        caltable='apcal',
        gaintype='T',
        refant=refant,
        calmode='ap',
        combine='spw',
        solint='inf',
        minsnr=3.0,
        minblperant=6,
#        uvrange='>50m', # may need to use to exclude extended emission
        gaintable='pcal3',
        spwmap=spwmap,
        solnorm=True)

plotms(vis='apcal',
       xaxis='time',
       yaxis='amp',
       iteraxis='antenna',
       plotrange=[0,0,0.2,1.8])

applycal(vis=contvis,
         spwmap=[spwmap,spwmap], # select which spws to apply the solutions for each table
         field=field,
         gaintable=['pcal3','apcal'],
         gainfield='',
         calwt=False,
         flagbackup=False,
         interp=['linearperobs','linearperobs'])

# Save the flags in case you need to go back to this step. 
flagmanager(vis=contvis,mode='save',versionname='after_apcal')

# Make amplitude and phase self-calibrated image.
for ext in ['.image','.mask','.model','.image.pbcor','.psf','.residual','.pb','.sumwt','.weight']:
    rmtables(contimagename + '_ap'+ ext)


tclean(vis=contvis,
       imagename=contimagename + '_ap',
       field=field,
       # phasecenter=phasecenter, # uncomment if mosaic or imaging an ephemeris object
       # mosweight = True, # uncomment if mosaic
       specmode='mfs',
       deconvolver='hogbom',
       # Uncomment the below to image with nterms>1.
       #deconvolver='mtmfs',
       #nterms=2,
       imsize = imsize, 
       cell= cell, 
       weighting = weighting, 
       robust=robust,
       niter=niter, 
       threshold=threshold, 
       interactive=True,
       gridder=gridder,
       savemodel='modelcolumn',
       pbcor=True,
       usepointing=False) # apply the primary beam correction since this is the last image.

#>>> Note final RMS and number of clean iterations. Compare the RMS to
#>>> the RMS from the earlier, pre-selfcal image.

# Save results of self-cal in a new ms
split(vis=contvis,
      outputvis=contvis+'.selfcal',
      datacolumn='corrected')

# reset the corrected data column in the  ms to the original calibration.

#>>> This can also be used to return your ms to it's original
#>>> pre-self-cal state if you are unhappy with your self-calibration.
clearcal(vis=contvis)

#>>> The applycal task will automatically flag data without good
#>>> gaincal solutions. If you are unhappy with your self-cal and wish to
#>>> return the flags to their original state, run the following command
#>>> flagmanager(vis=contvis, mode='restore',versionname='before_selfcal')


########################################
# Continuum Subtraction for Line Imaging

#>>> If you have observations that include both line and strong (>3 sigma
#>>> per final line image channel) continuum emission, you need to
#>>> subtract the continuum from the line data. You should not continuum
#>>> subtract if the line of interest is in absorption.

#>>> You can use au.invertChannelRanges(flagchannels,vis=finalvis) to
#>>> get the fitspw below. You will need to insert any continuum spws
#>>> that weren't included in flagchannels. For example, if your continuum
#>>> spws are '0,1,2' and flagchannels='1:260~500', au.invertChannelRanges will return
#>>> '1:0~259,1:501~3839'. The fitspw parameter should be '0,1:0~259,1:501~3839,2'
#>>> Make sure to cut and paste the output in fitspw below since PIs don't have
#>>> analysisUtilities by default.

fitspw = '2:0~1200;1500~3839,3:0~1200;1500~3839' # *line-free* channels for fitting continuum
linespw = '2,3' # line spectral windows. You can subtract the continuum from multiple spectral line windows at once. The linespw must be in the fitspw unless combine='spw'

finalvis='calibrated_final.ms'

uvcontsub(vis=finalvis,
          spw=linespw, # spw to do continuum subtraction on
          fitspw=fitspw, # regions without lines.
          excludechans=False, # fit the regions in fitspw
          #combine='spw', uncomment if there are no line-free channels in the line spectral window.  
          solint='int',
          fitorder=1,
          want_cont=False) # This value should not be changed.

#>>> Note that the continuum subtraction is done for each field in 
#>>> turn. However, if the fields have different line-free channels, you
#>>> will need to do the continuum subtraction separately for each field.

# NOTE: Imaging the continuum produced by uvcontsub with
# want_cont=True will lead to extremely poor continuum images because
# of bandwidth smearing effects. For imaging the continuum, you should
# always create a line-free continuum data set using the process
# outlined above.

#########################################################
# Apply continuum self-calibration to line data [OPTIONAL]

# uncomment one  of the following
#linevis = finalvis+'.contsub' # if continuum subtracted
# linevis = finalvis  #  if not continuum subtracted

# save original flags in case you don't like the self-cal
flagmanager(vis=linevis,mode='save',versionname='before_selfcal',merge='replace')

spwmap_line = [0] # Mapping self-calibration solution to the individual line spectral windows.
applycal(vis=linevis,
         spwmap=[spwmap_line, spwmap_line], # entering the appropriate spwmap_line value for each spw in the input dataset
         field=field,
         gaintable=['pcal3','apcal'],
         gainfield='',
         calwt=False,
         flagbackup=False,
         interp=['linearperobs','linearperobs'])

# Save results of self-cal in a new ms and reset the image name.
split(vis=linevis,
      outputvis=linevis+'.selfcal',
      datacolumn='corrected')

# reset the corrected data column in the  ms to the original calibration
#>>> This can also be used to return your ms to it's original
#>>> pre-self-cal state if you are unhappy with your self-calibration.
clearcal(linevis)

#>>> The applycal task will automatically flag data without good
#>>> gaincal solutions. If you are unhappy with your self-cal and wish to
#>>> return the flags to their original state, run the following command
#>>> flagmanager(vis=linevis, mode='restore',versionname='before_selfcal')

linevis=linevis+'.selfcal'

##############################################
# Image line emission [REPEAT AS NECESSARY]

#>>> If you did an mstransform/cvel, use the same velocity parameters in
#>>> the clean that you did for the regridding. If you did not do an
#>>> mstransform and have multiple executions of a scheduling block,
#>>> select the spws with the same rest frequency using the spw parameter
#>>> (currently commented out below). DO NOT INCLUDE SPWS WITH DIFFERENT
#>>> REST FREQUENCIES IN THE SAME RUN OF CLEAN: THEY WILL SLOW DOWN
#>>> IMAGING CONSIDERABLY.

finalvis = 'calibrated_final.ms'
# linevis = finalvis # uncomment if you neither continuum subtracted nor self-calibrated your data.
#linevis = finalvis + '.contsub' # uncomment if continuum subtracted
# linevis = finalvis + '.contsub.selfcal' # uncommment if both continuum subtracted and self-calibrated
# linevis = finalvis + '.selfcal' # uncomment if just self-calibrated (no continuum subtraction)



restfreq='115.27120GHz' # Typically the rest frequency of the line of
                        # interest. If the source has a significant
                        # redshift (z>0.2), use the observed sky
                        # frequency (nu_rest/(1+z)) instead of the
                        # rest frequency of the
                        # line.
#spw='1' # uncomment and replace with appropriate spw 

#>>> To specify a spws from multiple executions that had not been regridded using cvel, use
#>>>       import numpy as np
#>>>       spw = str.join(',',map(str,np.arange(0,n,nspw)))
#>>>
#>>> where n is the total number of windows x executions and nspw is the
#>>> number of spectral windows per execution. Note that the spectral
#>>> windows need to have the same order in all data sets for this code
#>>> to work. Add a constant offset (i.e., +1,+2,+3) to the array
#>>> generated by np.arange to get the other sets of windows.

#>>> Automatically generate image name for PI and paste info below
#>>>  aU.genImageName(vis=finalvis,spw=map(int,spw.split(',')),field=int(field.split('~')[0]),imtype='cube',targettype='sci',stokes='I',mous='',modtext='manual',spwmap=map(int,sciencespws.split(',')))

#>>>  If you don't have the variable sciencespws already set, you can
#>>>  determine it via a listobs on the original ms:
#>>>     listobs(vis='uid___A002_Xc3412f_X53ff.ms',intent='OBSERVE_TARGET*',spw='*FULL_RES*')
#>>>  Note that sciencespws needs to be a list of integers.

lineimagename =  ''

start='-100km/s' # start velocity. See science goals for appropriate value.
width='2km/s' # velocity width. See science goals.
nchan = 100  # number of channels. See science goals for appropriate value.

# If necessary, run the following commands to get rid of older clean
# data.

#clearcal(vis=linevis)
#delmod(vis=linevis)

for ext in ['.image','.mask','.model','.image.pbcor','.psf','.residual','.pb','.sumwt','.weight']:
    rmtables(lineimagename + ext)

tclean(vis=linevis,
       imagename=lineimagename, 
       field=field,
       spw=spw,
       # phasecenter=phasecenter, # uncomment if mosaic or imaging an ephemeris object   
       # mosweight = True, # uncomment if mosaic      
       specmode='cube', # comment this if observing an ephemeris source
       # specmode='cubesource', #uncomment this line if observing an ephemeris source
       perchanweightdensity=True, 
       start=start,
       width=width,
       nchan=nchan, 
       outframe=outframe,
       veltype=veltype, 
       restfreq=restfreq, 
       niter=niter,  
       threshold=threshold, 
       interactive=True,
       cell=cell,
       imsize=imsize, 
       weighting='briggsbwtaper',# a modified version of Briggs weighting for cubes, see CASA documentation for more info
       robust=robust,
       gridder=gridder,
       pbcor=True,
       restoringbeam='common',
       usepointing=False) 

#>>> If interactively cleaning (interactive=True), then note number of
#>>> iterations at which you stop for the PI. This number will help the
#>>> PI replicate the delivered images. Do not clean empty
#>>> images. Just click the red X to stop the interactive and note the
#>>> RMS.

# If you'd like to redo your clean, but don't want to make a new mask
# use the following commands to save your original mask. This is an
# optional step.
# linemaskname = 'line.mask'
## rmtables(linemaskname) # uncomment if you want to overwrite the mask.
# os.system('cp -ir ' + lineimagename + '.mask ' + linemaskname)

##############################################
# Export the images

import glob

myimages = glob.glob("*.pbcor")
for image in myimages:
    exportfits(imagename=image, fitsimage=image+'.fits',overwrite=True)

myimages = glob.glob("*.pb")
for image in myimages:
    exportfits(imagename=image, fitsimage=image+'.fits',overwrite=True) 

##############################################
# Create Diagnostic PNGs

os.system("rm -rf *.png")
mycontimages = glob.glob("*mfs*manual.image")
for cimage in mycontimages:
    mymax=imstat(cimage)['max'][0]
    mymin=-0.1*mymax
    outimage = cimage+'.png'
    os.system('rm -rf '+outimage)
    imview(raster={'file':cimage,'range':[mymin,mymax]},out=outimage)

mylineimages = glob.glob("*cube*manual.image")
for limage in mylineimages:
    mom8=limage+'.mom8'
    os.system("rm -rf "+mom8)
    immoments(limage,moments=[8],outfile=mom8)
    mymax=imstat(mom8)['max'][0]
    mymin=-0.1*mymax
    os.system("rm -rf "+mom8+".png")
    imview(raster={'file':mom8,'range':[mymin,mymax]},out=mom8+'.png')


##############################################
# Analysis

# For examples of how to get started analyzing your data, see
#     https://casaguides.nrao.edu/index.php/TWHydraBand7_Imaging_4.3
#