C#>>> =====================================================================================# #>>> TEMPLATE IMAGING SCRIPT # #>>> =====================================================================================# #>>> #>>> Updated: Fri Sep 28 09:51:22 EDT 2018 #>>> 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. # cd ######################################## # Check CASA version import re import glob ################################################## # 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. mses = glob.glob("calibrated*.ms") #mses = ['calibrated_X164.ms','calibrated_X1bc.ms','calibrated_X1c0.ms','calibrated_X1c4.ms','calibrated_X1c8.ms']#positions 10-14 #finalvis='calibrated_minimosaic_final.ms' # This is your output ms from the data # preparation script. finalvisall= ['../minimosaic/calibrated_X164.ms', '../minimosaic/calibrated_X1bc.ms', '../minimosaic/calibrated_X1c0.ms', '../minimosaic/calibrated_X1c4.ms', '../minimosaic/calibrated_X1c8.ms', '../minimosaic/calibrated_X1cc.ms', '../semipass/calibrated_X1d0.ms', '../minimosaic/calibrated_X1d4.ms', '../minimosaic/calibrated_X1d8.ms', '../minimosaic/calibrated_X184.ms', '../minimosaic/calibrated_X1e0.ms'] # 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' ) # Set spws to be used to form continuum contspws = '17,19,21,23' # Flag the "line channels" flagchannels='19:63~72' # only applicable to Tune147 data for thisvis in finalvisall: flagmanager(vis=thisvis,mode='save', versionname='before_cont_flags') initweights(vis=thisvis,wtmode='weight',dowtsp=True) if (thisvis != '../minimosaic/calibrated_X164.ms' and thisvis != '../minimosaic/calibrated_X184.ms'): flagdata(vis=thisvis,mode='manual', spw=flagchannels,flagbackup=False) else: print 'hi' for thisvis in ['../minimosaic/calibrated_X164.ms','../minimosaic/calibrated_X184.ms']: flagmanager(vis=thisvis,mode='restore', versionname='before_cont_flags') #>>> In CASA 4.4 and higher, the behavior of the avgchannel parameter #>>> has changed. Now when you plot binned channels, plotms displays #>>> the "bin" number rather than the average channel number of each #>>> bin. A ticket has been filed to revert this back to the previous #>>> (more sensible) behavior, but this behavior hasn't been fixed as of CASA 5.1. #>>> 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. # If you have complex line emission and no dedicated continuum # windows, you will need to flag the line channels prior to averaging. #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. #contvis='calibrated_minimosaic_final_cont.ms' #rmtables(contvis) #os.system('rm -rf ' + contvis + '.flagversions') contvisall= ['calibrated_X164_cont.ms', 'calibrated_X1bc_cont.ms', 'calibrated_X1c0_cont.ms', 'calibrated_X1c4_cont.ms', 'calibrated_X1c8_cont.ms', 'calibrated_X1cc_cont.ms', 'calibrated_X1d0_cont.ms', 'calibrated_X1d4_cont.ms', 'calibrated_X1d8_cont.ms', 'calibrated_X184_cont.ms', 'calibrated_X1e0_cont.ms'] conttimevisall= ['calibrated_X164_conttime.ms', 'calibrated_X1bc_conttime.ms', 'calibrated_X1c0_conttime.ms', 'calibrated_X1c4_conttime.ms', 'calibrated_X1c8_conttime.ms', 'calibrated_X1cc_conttime.ms', 'calibrated_X1d0_conttime.ms', 'calibrated_X1d4_conttime.ms', 'calibrated_X1d8_conttime.ms', 'calibrated_X184_conttime.ms', 'calibrated_X1e0_conttime.ms'] #>>> 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. for i in range(0,11): split(vis=finalvisall[i], spw=contspws, outputvis=conttimevisall[i], width=[8,8,8,8], timebin='10s', datacolumn='data') for i in range(0,11): split(vis=finalvisall[i], spw=contspws, outputvis=contvisall[i], width=[8,8,8,8], datacolumn='data') # Check the weights. You will need to change antenna and field to # appropriate values plotms(vis=contvisall[0], yaxis='wtsp',xaxis='freq',spw='',antenna='DA43',field='0') # If you flagged any line channels, restore the previous flags for thisvis in finalvisall: flagmanager(vis=thisvis,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 # ------------------ contvis= ['calibrated_X1c8_cont.ms', 'calibrated_X1cc_cont.ms', 'calibrated_X1d0_cont.ms', 'calibrated_X1d4_cont.ms', 'calibrated_X1d8_cont.ms'] conttimevis= ['calibrated_X1c8_conttime.ms', 'calibrated_X1cc_conttime.ms', 'calibrated_X1d0_conttime.ms', 'calibrated_X1d4_conttime.ms', 'calibrated_X1d8_conttime.ms'] field=['3~13,74~84','3~13,74~84','3~13,74~84','3~13,74~84','3~13,74~84'] #field=['28~33,99~104','28~33,99~104','28~33,99~104','28~33,99~104','28~33,99~104'] # 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 = 'J2000 10h00m15.91 02d12m50.4' # 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='0.2arcsec' # cell size for imaging. imsize = [1300,1300] # 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=2.0 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=map(int,contspws.split(',')),field=int(field.split('~')[0]),imtype='mfs',targettype='sci',stokes='I',mous='',modtext='manual',spwmap=map(int,sciencespws.split(','))) contimagename = 'mosaictest_notimeav_cmc' # If necessary, run the following commands to get rid of older clean # data. for i in range(0,5): clearcal(vis=contvisall[i]) delmod(vis=contvisall[i]) 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, mosweight = True, specmode='mfs', deconvolver='hogbom', imsize = imsize, cell= cell, weighting = weighting, robust = robust, niter = niter, threshold = threshold, interactive = True, gridder = gridder, pbcor = True) #>>> 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) ############################################## ############################################## # 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 #