### For running a selection of a few steps use: #my_steps=[0,1,2,3,4,5,6,7,8,9,10] ### For running a row of many steps use: #my_steps=range(0,10) # from 0 to 9 my_vis = 'xx' # ms-file name targetname = 'xx' # choose name for image praefix fld=targetname if 0 in my_steps : ######################################## ### Inspection for spectral line ### ######################################## ## Visually inspect amplitudes vs. frequency of our scientific target in order to identify ## emission line feature(s) os.system('mkdir '+my_vis+'.split.inspect.plots') # make new folder for plots plotms(vis=my_vis, # ms name spw='', xaxis='channel', yaxis='amp', field=fld, # target avgtime='1e8', avgscan=True, showatm=True, # Show the atmospheric transmission? iteraxis='spw', # plot for each SpW (iteration) - redo for iteration over scan (set also avgscan=False!) coloraxis='corr', # colour code for data points showgui=True, # open user interface plotfile=my_vis+'.split.inspect.plots/field5-amp-channel.png', # export the plot - iteration name will be inserted automatically - change field name manually! exprange='all', # make plot for every iteration overwrite=True) ## user_check=raw_input('When you are done with the window, close it and press enter to continue:') ## CO line: spw xx ## unknown: xx ## Spectral window xx is devoid of any line emission -> pure continuum; let's check if continuum ## source is extended plotms(vis=my_vis, # ms name spw='', xaxis='uvwave', # uvdist[m]/wavelength yaxis='amp', field=fld, # target avgtime='1e8', avgchannel='1e8', coloraxis='baseline', iteraxis='spw', # plot for each SpW (iteration) - redo for iteration over scan (set also avgscan=False!) #coloraxis='corr', # colour code for data points showgui=True, # open user interface plotfile=my_vis+'.split.inspect.plots/field5-amp-uvdist.png', # export the plot - iteration name will be inserted automatically - change field name manually! exprange='all', # make plot for every iteration overwrite=True) ## In addition, note the largest baseline in the plot in klambda ## -> corresponds to an angular resolution of theta[arcsec] ~ 206/baseline[klambda] ## ## b = xx klambda ---> theta = xx'' (spw xx) ## ## We will need this information to set up the pixel scale of the image. ## user_check=raw_input('When you are done with the window, close it and press enter to continue:') if 1 in my_steps : #################################### ### IMAGING: DIRTY continuum ### #################################### ###################### dirty continuum map ################################ ## ## Let's obtain a dirty map of the continuum emission, excluding all line emission and ## avoiding higher noise channels os.system('rm -rf '+targetname+'_cont_dirty.*') # delete old file to start from scratch, else clean will continue to work on this image ! tclean(vis=my_vis, # ms name imagename=targetname+'_cont_dirty', # image name spw='xx', # channel selection - avoid spectral features and noisy channels specmode = '', # continuum or spectral #nterms = 2, # parameter for mtmfs #deconvolver = 'mtmfs', # mtmfs needed, when total bandwidth (max-min frequency covered by all spw)/observing frequency > 10% field=targetname, # target niter=0, # number of iterations; niter=0 --> dirty image (direct transform of visibilities) imsize=xx, # number of pixels per image axis - depends primary beam size (field-of-view) - try it out or estimate cell='xxarcsec', # size of a pixel - typically a beam/PSF axis should be sampled by ~5 pixels weighting='briggs', # weighting scheme robust=0.0, # briggs weighting scheme parameter interactive=False) # work in a GUI or not? ## logger: Fitted beam used in restoration: 0.808463 by 0.486767 (arcsec) at pa -86.7345 (deg) ## ## Use the viewer to inspect the dirty continuum image and find an emission-free region viewer() if 2 in my_steps : #################################### ### IMAGING: CLEAN continuum ### #################################### ############## Obtain the statistics for the dirty continuum map ############# ## ## Get an idea of the rms-noise in an emission-free corner of the image ## ## Viewer: Use the region marker to select a part of your image (box, ellipse, polygon), ## and look at 'regions panel' (might have to activate it in 'view -> regions') ## ## rms from viewer: xx mJy ## ## ## - or - ## Script/Terminal: A more reproducible approach: stats_cont_dirty = imstat(imagename=targetname+'_cont_dirty.image', # image name box='xx') # define a box-shaped region ('x_min,y_min,x_max,y_max') where to obtain the statistics from ## imstat, like many analysis tasks in CASA, creates a dictionary of parameters - skim through it to get an impression of it (beware units): print('\n Dictionary of the statistics of the dirty continuum image: \n\n', stats_cont_dirty) # \n means 'new line' ## Get the rms-noise of dirty continuum image in mJy/beam: ## Select the first entry of the 'rms'-list rms_cont_dirty=stats_cont_dirty["rms"][0]*1e3 ## define a cleaning threshold = maximum remaining residual when cleaning clthr_cont=3*rms_cont_dirty # a source is called 'detected' when its flux > 3*sigma - still noise peaks till 5*sigma possible - use conservative limit for extended emission! clthr_cont="{:.1f}".format(clthr_cont)+'mJy' # reformat the variable entry into a CASA-digestible shape ###################### clean continuum map ################################ ## ## Let's obtain a clean map of the continuum emission, excluding all line emission and ## avoiding higher noise channels. We will use a rough clean mask os.system('rm -rf '+targetname+'_cont.*') # delete old file to start from scratch, else clean will continue to work on this image ! tclean(vis=my_vis, # ms name imagename=targetname+'_cont', # image name spw='xx', # channel selection - avoid spectral features and noisy channels field=targetname, # target specmode = 'xx', # continuum or spectral #nterms = 2, # parameter for mtmfs #deconvolver = 'mtmfs', # mtmfs needed, when total bandwidth (max-min frequency covered by all spw)/observing frequency > 10% niter=10000, # number of iterations; niter=0 --> dirty image (direct transform of visibilities) threshold=clthr_cont, # clean threshold - maximum remaining residual when cleaning mask='box[[xxpix,xxpix],[xxpix,xxpix]]', # clean mask - help CLEAN to focus on what you think is real emission; mask = [x_min, y_min, x_max, y_max] imsize=xx, # number of pixels per image axis - depends on primary beam size (field-of-view) - try it out or estimate cell='xxarcsec', # size of a pixel - typically a beam/PSF axis should be sampled by ~5 pixels weighting='briggs', # weighting scheme robust=0.0, # briggs weighting scheme parameter interactive=False) # work in a GUI or not? ############## Obtain the statistics for the clean continuum map ############# ## ## Get the rms: stats_cont = imstat(imagename=targetname+'_cont.image', box='xx') # emission-free region rms_cont=stats_cont["rms"][0] # rms-noise of continuum image in mJy/beam ## Get the peak-flux in the image: stats_cont = imstat(imagename=targetname+'_cont.image', box='xx') # reuse the clean mask is sufficient peak_cont=stats_cont["max"][0] # peak-flux of continuum image in mJy/beam ## Derive the dynamic range: dyn_cont=peak_cont/rms_cont # dynamic range in continuum image print('\nStatistics of the clean continuum image: \nRMS:', rms_cont*1e3, 'mJy - peak:', peak_cont*1e3, 'mJy - dynamic range:', dyn_cont) # see your requested statistics in the terminal ## Output a png-image of the clean continuum maps imview(raster={'file': targetname+'_cont.image', # image name 'colorwedge':True, # color bar 'range':[-rms_cont, peak_cont], # flux values from ... to ... 'scaling':0, # brightness scaling 'colormap':'Rainbow 2'}, # map color scheme out=targetname+'_cont.png', # output image zoom=3) # zoom factor if 3 in my_steps : ###################### dirty continuum map at different uv weighting (natural) ################################ ## ## Let's obtain a lower resolution dirty map of the continuum emission (using close to natural ## weighting, excluding all line emission and avoiding higher noise channels. os.system('rm -rf '+targetname+'_cont_dirty_nw*') # delete old file to start from scratch, else clean will continue to work on this image ! tclean(vis=my_vis, # ms name imagename=targetname+'_cont_dirty_nw', # image name spw='xx', # channel selection - avoid spectral features and noisy channels specmode = 'xx', # continuum or spectral #nterms = 2, # parameter for mtmfs #deconvolver = 'mtmfs', # mtmfs needed, when total bandwidth (max-min frequency covered by all spw)/observing frequency > 10% field=targetname, # target niter=0, # number of iterations; niter=0 --> dirty image (direct transform of visibilities) imsize=xx, # number of pixels per image axis - depends on primary beam size (field-of-view) - try it out or estimate cell='xxarcsec', # size of a pixel - typically a beam/PSF axis should be sampled by ~5 pixels weighting='briggs', # weighting scheme robust=xx, # briggs weighting scheme parameter to 'natural' interactive=False) # work in a GUI or not? ############## Obtain the statistics for the newly weighted dirty continuum map ############# ## ## Get an idea of the rms-noise in an emission-free corner of the image stats_cont_dirty_nw = imstat(imagename=targetname+'_cont_dirty_nw.image', box='xx') rms_cont_dirty_nw=stats_cont_dirty_nw["rms"][0]*1e3 # rms-noise of dirty continuum image in mJy/beam clthr_cont_nw=3*rms_cont_dirty_nw # a source is called 'detected' when its flux > 3*sigma - still noise peaks till 5*sigma possible - use conservative limit for extended emission! clthr_cont_nw="{:.1f}".format(clthr_cont_nw)+'mJy' ###################### clean continuum map at different uv weighting (natural)################################ ## ## Let's obtain a lower resolution clean map of the continuum emission, excluding all line emission and ## avoiding higher noise channels. We will use a rough clean mask. We will use a robust parameter close ## to natural weighting os.system('rm -rf '+targetname+'_cont_nw.*') # delete old file to start from scratch, else clean will continue to work on this image ! tclean(vis=my_vis, # ms name imagename=targetname+'_cont_nw', # image name spw='xx', # channel selection - avoid spectral features and noisy channels specmode = 'xx', # continuum or spectral #nterms = 2, # parameter for mtmfs #deconvolver = 'mtmfs', # mtmfs needed, when total bandwidth (max-min frequency covered by all spw)/observing frequency > 10% field=targetname, # target niter=500, # number of iterations; niter=0 --> dirty image (direct transform of visibilities) threshold=clthr_cont_nw, # clean threshold - maximum remaining residual when cleaning mask='box[[xxpix,xxpix],[xxpix,xxpix]]', # clean mask - help CLEAN to focus on what you think is real emission; mask = [x_min, y_min, x_max, y_max] imsize=xx, # number of pixels per image axis - depends on primary beam size (field-of-view) - try it out or estimate cell='xxarcsec', # size of a pixel - typically a beam/PSF axis should be sampled by ~5 pixels weighting='briggs', # weighting scheme robust=xx, # weighting scheme parameter to 'natural' interactive=False) # work in a GUI or not? ############## Obtain the statistics for the newly weighted clean continuum map ############# ## ## Get the rms: stats_cont_nw = imstat(imagename=targetname+'_cont_nw.image', box='xx') rms_cont_nw=stats_cont_nw["rms"][0] # rms-noise of continuum image in mJy/beam ## Get the peak-flux in the image: stats_cont_nw = imstat(imagename=targetname+'_cont_nw.image', box='xx') peak_cont_nw=stats_cont_nw["max"][0] # peak-flux of continuum image in mJy/beam ## Derive the dynamic range: dyn_cont_nw=peak_cont_nw/rms_cont_nw # dynamic range in continuum image print('\nStatistics of the naturally weighted clean continuum image: \nRMS:', rms_cont_nw*1e3, 'mJy - peak:', peak_cont_nw*1e3, 'mJy - dynamic range:', dyn_cont_nw) # see your requested statistics in the termianl ## Output a png-image of the clean continuum maps imview(raster={'file': targetname+'_cont_nw.image', # image name 'colorwedge':True, # color bar 'range':[-rms_cont_nw, peak_cont_nw], # flux values from ... to ... 'scaling':0, # brightness scaling 'colormap':'Rainbow 2'}, # map color scheme out=targetname+'_cont_nw.png', # output image zoom=3) # zoom factor ###################### dirty continuum map at different uv weighting (uniform) ################################ ## ## Let's obtain a lower resolution dirty map of the continuum emission (using close to uniform ## weighting, excluding all line emission and avoiding higher noise channels. os.system('rm -rf '+targetname+'_cont_dirty_uw*') # delete old file to start from scratch, else clean will continue to work on this image ! tclean(vis=my_vis, # ms name imagename=targetname+'_cont_dirty_uw', # image name spw='xx', # channel selection - avoid spectral features and noisy channels specmode = 'xx', # continuum or spectral #nterms = 2, # parameter for mtmfs #deconvolver = 'mtmfs', # mtmfs needed, when total bandwidth (max-min frequency covered by all spw)/observing frequency > 10% field=targetname, # target niter=0, # number of iterations; niter=0 --> dirty image (direct transform of visibilities) imsize=xx, # number of pixels per image axis - depends on primary beam size (field-of-view) - try it out or estimate cell='xxarcsec', # size of a pixel - typically a beam/PSF axis should be sampled by ~5 pixels weighting='briggs', # weighting scheme robust=xx, # briggs weighting scheme parameter to 'natural' interactive=False) # work in a GUI or not? ############## Obtain the statistics for the newly weighted dirty continuum map ############# ## ## Get an idea of the rms-noise in an emission-free corner of the image stats_cont_dirty_uw = imstat(imagename=targetname+'_cont_dirty_uw.image', box='xx') rms_cont_dirty_uw=stats_cont_dirty_uw["rms"][0]*1e3 # rms-noise of dirty continuum image in mJy/beam clthr_cont_uw=3*rms_cont_dirty_uw # a source is called 'detected' when its flux > 3*sigma - still noise peaks till 5*sigma possible - use conservative limit for extended emission! clthr_cont_uw="{:.1f}".format(clthr_cont_uw)+'mJy' ###################### clean continuum map at different uv weighting (uniform) ################################ ## ## Let's obtain a lower resolution clean map of the continuum emission, excluding all line emission and ## avoiding higher noise channels. We will use a rough clean mask. We will use a robust parameter close ## to uniform weighting os.system('rm -rf '+targetname+'_cont_uw.*') # delete old file to start from scratch, else clean will continue to work on this image ! tclean(vis=my_vis, # ms name imagename=targetname+'_cont_uw', # image name spw='xx', # channel selection - avoid spectral features and noisy channels specmode = 'xx', # continuum or spectral #nterms = 2, # parameter for mtmfs #deconvolver = 'mtmfs', # mtmfs needed, when total bandwidth (max-min frequency covered by all spw)/observing frequency > 10% field=targetname, # target niter=500, # number of iterations; niter=0 --> dirty image (direct transform of visibilities) threshold=clthr_cont_uw, # clean threshold - maximum remaining residual when cleaning mask='box[[xxpix,xxpix],[xxpix,xxpix]]', # clean mask - help CLEAN to focus on what you think is real emission; mask = [x_min, y_min, x_max, y_max] imsize=xx, # number of pixels per image axis - depends on primary beam size (field-of-view) - try it out or estimate cell='xxarcsec', # size of a pixel - typically a beam/PSF axis should be sampled by ~5 pixels weighting='briggs', # weighting scheme robust=xx, # weighting scheme parameter to 'natural' interactive=False) # work in a GUI or not? ############## Obtain the statistics for the newly weighted clean continuum map ############# ## ## Get the rms: stats_cont_uw = imstat(imagename=targetname+'_cont_uw.image', box='xx') rms_cont_uw=stats_cont_uw["rms"][0] # rms-noise of continuum image in mJy/beam ## Get the peak-flux in the image: stats_cont_uw = imstat(imagename=targetname+'_cont_uw.image', box='xx') peak_cont_uw=stats_cont_uw["max"][0] # peak-flux of continuum image in mJy/beam ## Derive the dynamic range: dyn_cont_uw=peak_cont_uw/rms_cont_uw # dynamic range in continuum image print('\nStatistics of the uniformly weighted clean continuum image: \nRMS:', rms_cont_uw*1e3, 'mJy - peak:', peak_cont_uw*1e3, 'mJy - dynamic range:', dyn_cont_uw) # see your requested statistics in the termianl ## Output a png-image of the clean continuum maps imview(raster={'file': targetname+'_cont_uw.image', # image name 'colorwedge':True, # color bar 'range':[-rms_cont_uw, peak_cont_uw], # flux values from ... to ... 'scaling':0, # brightness scaling 'colormap':'Rainbow 2'}, # map color scheme out=targetname+'_cont_uw.png', # output image zoom=3) # zoom factor