IRAS16293 Band9 - Imaging for CASA 4.5: Difference between revisions
Created page with "<pre style="background-color: #ffa07a;"> WARNING: On June 15, 2012 the calibration guide and the final data products (calibrated science data: IRAS16293_Band9_CalibratedMS_FI..." |
No edit summary |
||
(6 intermediate revisions by one other user not shown) | |||
Line 46: | Line 46: | ||
# listobs to check data | # listobs to check data | ||
listobs(vis='IRAS16293_Band9.fixed.rebin.ms', | listobs(vis='IRAS16293_Band9.fixed.rebin.ms', | ||
listfile='IRAS16293_Band9.fixed.rebin.ms.listobs',verbose= | listfile='IRAS16293_Band9.fixed.rebin.ms.listobs',verbose=False) | ||
</source> | </source> | ||
Line 73: | Line 73: | ||
</pre> | </pre> | ||
''' | '''Unless you have calibrated yourself the measurement set you are using for imaging, it is essential to do the following.''' | ||
<source lang="python"> | <source lang="python"> | ||
Line 118: | Line 118: | ||
</source> | </source> | ||
We need to note the channels that do not present line emission. To identify them we can plot the amplitude vs frequency for all the spectral windows. Figure 4 shows the output from the next command for spw | We need to note the channels that do not present line emission. To identify them we can plot the amplitude vs frequency for all the spectral windows. Figure 4 shows the output from the next command for spw 1 (field 3). You are dealing with a set of data which were split so the calibrated data now reside in the "data" column of the new measurement set. Therefore the following call to plotms must specify ydatacolumn = "data". If you had never performed the split, the calibrated data would have been accessed from the "corrected" column of the old, unsplit, measurement set. Those who have used previous versions of CASA should note that in those older versions, plotms would default to using the "data" column if there were no data in the "corrected" column and one specified ydatacolumn="corrected". | ||
[[File:amp_vs_channel.spw1.png|200px|thumb|right|'''Fig. 4.''' Frequency-amplitude plot for spw=1 of field=3. ]] | [[File:amp_vs_channel.spw1.png|200px|thumb|right|'''Fig. 4.''' Frequency-amplitude plot for spw=1 of field=3. ]] | ||
Line 129: | Line 128: | ||
plotms(vis='IRAS16293_Band9.fixed.rebin.ms', | plotms(vis='IRAS16293_Band9.fixed.rebin.ms', | ||
xaxis='channel', yaxis='amp',field='3', | xaxis='channel', yaxis='amp',field='3', | ||
spw='', avgtime='1e8',avgscan= | spw='', avgtime='1e8',avgscan=True,coloraxis='corr', | ||
avgchannel='',iteraxis='spw',ydatacolumn='data', | avgchannel='',iteraxis='spw',ydatacolumn='data', | ||
xselfscale=True,yselfscale=True) | xselfscale=True,yselfscale=True) | ||
Line 147: | Line 146: | ||
# flagging the line channels | # flagging the line channels | ||
flagdata(vis='IRAS16293_Band9.fixed.rebin.ms', | flagdata(vis='IRAS16293_Band9.fixed.rebin.ms', | ||
spw='0:1700~2500,1:400~500;1100~1900;3000~3450,2:1700~2200;3100~3839,3:0~800;1600~2900;3300~3839',action='apply', flagbackup= | spw='0:1700~2500,1:400~500;1100~1900;3000~3450,2:1700~2200;3100~3839,3:0~800;1600~2900;3300~3839',action='apply', flagbackup=False) | ||
</source> | </source> | ||
Line 170: | Line 169: | ||
mode='mfs',imagermode='mosaic', | mode='mfs',imagermode='mosaic', | ||
imsize=500,cell='0.04arcsec',minpb=0.25, | imsize=500,cell='0.04arcsec',minpb=0.25, | ||
interactive= | interactive=True, | ||
weighting='briggs',robust=0.5, | weighting='briggs',robust=0.5, | ||
niter=10000, threshold='0.3mJy',usescratch= | niter=10000, threshold='0.3mJy',usescratch=False) | ||
</source> | </source> | ||
Now we proceed with {{gaincal}} to perform self calibration (for phases) on the data. Even with minsnr=3.0 there will be some failed solutions. | Now we proceed with {{gaincal}} to perform self calibration (for phases) on the data. Even with minsnr=3.0 there will be some failed solutions. | ||
Line 213: | Line 202: | ||
# application of the table | # application of the table | ||
applycal(vis='IRAS16293_Band9.fixed.CONT.msAVG',spwmap=[],spw='', | applycal(vis='IRAS16293_Band9.fixed.CONT.msAVG',spwmap=[],spw='', | ||
gaintable=['cont_pcal1'],calwt= | gaintable=['cont_pcal1'],calwt=False,flagbackup=False) | ||
</source> | </source> | ||
Line 226: | Line 215: | ||
mode='mfs',imagermode='mosaic', | mode='mfs',imagermode='mosaic', | ||
imsize=500,cell='0.04arcsec',minpb=0.25, | imsize=500,cell='0.04arcsec',minpb=0.25, | ||
interactive= | interactive=True, | ||
weighting='briggs',robust=0.5, | weighting='briggs',robust=0.5, | ||
mask='IRAS16293_Band9.fixed.CONTIN.mask', | mask='IRAS16293_Band9.fixed.CONTIN.mask', | ||
niter=10000, threshold='0.3mJy', usescratch= | niter=10000, threshold='0.3mJy', usescratch=False) | ||
# visualize the image | # visualize the image | ||
Line 256: | Line 245: | ||
# Applying results to the data | # Applying results to the data | ||
applycal(vis='IRAS16293_Band9.fixed.CONT.msAVG',spwmap=[],spw='', | applycal(vis='IRAS16293_Band9.fixed.CONT.msAVG',spwmap=[],spw='', | ||
gaintable=['cont_pcal2'],calwt= | gaintable=['cont_pcal2'],calwt=False,flagbackup=False) | ||
</source> | </source> | ||
Line 269: | Line 258: | ||
mode='mfs',imagermode='mosaic', | mode='mfs',imagermode='mosaic', | ||
imsize=500,cell='0.04arcsec',minpb=0.25, | imsize=500,cell='0.04arcsec',minpb=0.25, | ||
interactive= | interactive=True, | ||
weighting='briggs',robust=0.5, | weighting='briggs',robust=0.5, | ||
mask='IRAS16293_Band9.fixed.CONTIN.mask', | mask='IRAS16293_Band9.fixed.CONTIN.mask', | ||
niter=10000, threshold='0.3mJy', usescratch= | niter=10000, threshold='0.3mJy', usescratch=False) | ||
# to view the resulting image | # to view the resulting image | ||
Line 295: | Line 284: | ||
solint='8000s',minsnr=3.0,minblperant=4,gaintable='cont_pcal2') | solint='8000s',minsnr=3.0,minblperant=4,gaintable='cont_pcal2') | ||
# | # plotting the table | ||
plotcal(caltable='cont_apcal',xaxis='time',yaxis='amp', | plotcal(caltable='cont_apcal',xaxis='time',yaxis='amp', | ||
spw='',iteration='antenna',subplot=511,plotrange=[0,0,0.5,2], | spw='',iteration='antenna',subplot=511,plotrange=[0,0,0.5,2], | ||
Line 302: | Line 291: | ||
# and applying to the data | # and applying to the data | ||
applycal(vis='IRAS16293_Band9.fixed.CONT.msAVG',spwmap=[],spw='', | applycal(vis='IRAS16293_Band9.fixed.CONT.msAVG',spwmap=[],spw='', | ||
gaintable=['cont_pcal2','cont_apcal'],calwt= | gaintable=['cont_pcal2','cont_apcal'],calwt=False,flagbackup=False) | ||
# A large correction is present for the 2nd dataset. This amplitude | # A large correction is present for the 2nd dataset. This amplitude | ||
Line 319: | Line 308: | ||
mode='mfs',imagermode='mosaic', | mode='mfs',imagermode='mosaic', | ||
imsize=500,cell='0.04arcsec',minpb=0.25, | imsize=500,cell='0.04arcsec',minpb=0.25, | ||
interactive= | interactive=True, | ||
weighting='briggs',robust=0.5, | weighting='briggs',robust=0.5, | ||
mask='IRAS16293_Band9.fixed.CONTIN_P2.mask', | mask='IRAS16293_Band9.fixed.CONTIN_P2.mask', | ||
niter=10000, threshold='0.3mJy', usescratch= | niter=10000, threshold='0.3mJy', usescratch=False) | ||
# and to view the image | # and to view the image | ||
Line 352: | Line 341: | ||
In Figure 7 we show the final image of the self-calibration process. Source A is the weaker source to the SE and Source B is the stronger source to the NW. The rms noise is about 18 mJy/beam, and the restoring beam is 0.29 x 0.17 arcsec, PA=-71 deg with robust=0.5 weighting. The peak intensities for the two components are about 3.19 (B) and 1.07 (A) Jy/beam, while the flux densities are about 11.9 (B) and 9 (A) Jy. | In Figure 7 we show the final image of the self-calibration process. Source A is the weaker source to the SE and Source B is the stronger source to the NW. The rms noise is about 18 mJy/beam, and the restoring beam is 0.29 x 0.17 arcsec, PA=-71 deg with robust=0.5 weighting. The peak intensities for the two components are about 3.19 (B) and 1.07 (A) Jy/beam, while the flux densities are about 11.9 (B) and 9 (A) Jy. | ||
[[File:IRAS16293_Band9.fixed.CONTIN_AP.image.pbcor.subim_c.png|200px|thumb|right|'''Fig. 7.''' Continuum image resulted from the complete process of self calibration for the mosaic of IRAS16293 Band 9.]] | [[File:IRAS16293_Band9.fixed.CONTIN_AP.image.pbcor.subim_c.png|200px|thumb|right|'''Fig. 7.''' Continuum image resulted from the complete process of self calibration for the mosaic of IRAS16293 Band 9.]] | ||
Line 372: | Line 358: | ||
</source> | </source> | ||
Then we | Then we subtract the continuum, using the line free channels | ||
<source lang="python"> | <source lang="python"> | ||
Line 385: | Line 371: | ||
# Apply final self-calibration to the line data | # Apply final self-calibration to the line data | ||
applycal(vis='IRAS16293_Band9.fixed.rebin.ms.contsub',spwmap=[],spw='', | applycal(vis='IRAS16293_Band9.fixed.rebin.ms.contsub',spwmap=[],spw='', | ||
gaintable=['cont_pcal2','cont_apcal'],calwt= | gaintable=['cont_pcal2','cont_apcal'],calwt=False,flagbackup=False) | ||
</source> | </source> | ||
Some notes for this sections: the threshold was chosen after doing a couple of smaller test images. This is about 3 sigma in channels with weak line emission. The mask was also chosen from a few test images. | Some notes for this sections: the threshold was chosen after doing a couple of smaller test images. This is about 3 sigma in channels with weak line emission. The mask was also chosen from a few test images. | ||
=== Create a velocity cube for H13CN (J=8-7) === | === Create a velocity cube for H13CN (J=8-7) === | ||
Line 446: | Line 389: | ||
plotms(vis='IRAS16293_Band9.fixed.rebin.ms.contsub', | plotms(vis='IRAS16293_Band9.fixed.rebin.ms.contsub', | ||
xaxis='velocity', yaxis='amp',field='', | xaxis='velocity', yaxis='amp',field='', | ||
spw='1', avgtime='1e8',avgscan= | spw='1', avgtime='1e8',avgscan=True,coloraxis='field', | ||
transform=True,freqframe='LSRK',restfreq='690.55207GHz', | transform=True,freqframe='LSRK',restfreq='690.55207GHz', | ||
avgchannel='8', | avgchannel='8', | ||
Line 454: | Line 397: | ||
Now do an interactive clean, making clean boxes for the emission. Stop | Now do an interactive clean, making clean boxes for the emission. Stop | ||
cleaning when noise like residuals remain. Number of iterations | cleaning when noise like residuals remain. Number of iterations | ||
between showing interactive window can be controlled within the | between showing each interactive window can be controlled within the | ||
interactive viewer. We used about 400 iterations. | interactive viewer. We used about 400 iterations. | ||
Line 466: | Line 409: | ||
imagermode='mosaic', | imagermode='mosaic', | ||
imsize=500,cell='0.04arcsec',minpb=0.25, | imsize=500,cell='0.04arcsec',minpb=0.25, | ||
interactive= | interactive=True, | ||
weighting='briggs',robust=0.5, | weighting='briggs',robust=0.5, | ||
niter=500, threshold='220.0mJy') | niter=500, threshold='220.0mJy') | ||
Line 475: | Line 418: | ||
<source lang="python"> | <source lang="python"> | ||
impbcor(imagename='IRAS16293_Band9.fixed.H13CN_8_7.image', | impbcor(imagename='IRAS16293_Band9.fixed.H13CN_8_7.image', | ||
Line 499: | Line 429: | ||
<source lang="python"> | <source lang="python"> | ||
exportfits(imagename='IRAS16293_Band9.fixed.H13CN_8_7.image.pbcor.subim', | exportfits(imagename='IRAS16293_Band9.fixed.H13CN_8_7.image.pbcor.subim', |
Latest revision as of 13:31, 4 April 2017
WARNING: On June 15, 2012 the calibration guide and the final data products (calibrated science data: IRAS16293_Band9_CalibratedMS_FIXED.tgz and reference images: IRAS16293_Band9_ReferenceImages_FIXED.tgz)) were changed to correct for a 1.2" position error in the phase calibrator (1625-254). Without correction, the science images will suffer from a similar offset.
Overview
This section of the casa guide continues from IRAS16293_Band9_-_Calibration_for_CASA_4.5 . If you completed the calibration section, then you already have the final calibrated dataset. If you just started to work on this IRAS16293 B9 data, you can download the fully calibrated dataset IRAS16293_Band9_CalibratedMS from
This casa guide covers imaging of both continuum and spectral lines of IRAS16293 Band 9. This guide requires CASA 4.5. In this casa guide you will also need the Analysis Utils package, so make sure you have it imported into CASA.
Confirm your version of CASA
This guide has been written for CASA release 4.5.0. Please confirm your version before proceeding.
# In CASA
version = casadef.casa_version
print "You are using " + version
if (version < '4.5.0'):
print "YOUR VERSION OF CASA IS TOO OLD FOR THIS GUIDE."
print "PLEASE UPDATE IT BEFORE PROCEEDING."
else:
print "Your version of CASA is appropriate for this guide."
Install Analysis Utilities
Analysis Utilities (or analysisUtils for short) is a small set of Python scripts that provide a number of analysis and plotting utilities for ALMA data reduction. This guide uses a few of these utilities. They are very easy to install (just download and untar). See
http://casaguides.nrao.edu/index.php?title=Analysis_Utilities
for a full description and download instructions. If you do not wish to do this, see the CASA 3.3 version of the guide linked at the top of this page for alternative (but slow) plotting options. Analysis Utilities are updated frequently so if its been a while since you installed it, its probably worth doing it again. If you are at an ALMA site or ARC, the analysis utilities are probably already installed and up to date.
List the contents of the Data
Before we start with the imaging itself, do a listobs to the calibrated dataset, so you can check that everything is ok.
# listobs to check data
listobs(vis='IRAS16293_Band9.fixed.rebin.ms',
listfile='IRAS16293_Band9.fixed.rebin.ms.listobs',verbose=False)
You should see the next information in the .listobs file.
Fields: 7 ID Code Name RA Decl Epoch SrcId nRows 0 none IRAS16293-2422-a 16:32:22.99200 -24.28.36.0000 J2000 0 7811 1 none IRAS16293-2422-a 16:32:22.47925 -24.28.36.0000 J2000 0 7829 2 none IRAS16293-2422-a 16:32:22.73563 -24.28.36.0000 J2000 0 7829 3 none IRAS16293-2422-a 16:32:22.73563 -24.28.32.5000 J2000 0 7829 4 none IRAS16293-2422-a 16:32:22.47925 -24.28.29.0000 J2000 0 7429 5 none IRAS16293-2422-a 16:32:22.73563 -24.28.29.0000 J2000 0 7151 6 none IRAS16293-2422-a 16:32:22.99200 -24.28.29.0000 J2000 0 7151 Spectral Windows: (4 unique spectral windows and 1 unique polarization setups) SpwID #Chans Frame Ch1(MHz) ChanWid(kHz) TotBW(kHz) Corrs 0 3840 TOPO 703312.744 488.28125 1875000 XX YY 1 3840 TOPO 692237.256 -488.28125 1875000 XX YY 2 3840 TOPO 690437.256 -488.28125 1875000 XX YY 3 3840 TOPO 688437.256 -488.28125 1875000 XX YY Antennas: 13 'name'='station' ID= 0-3: 'DA41'='A003', 'DA43'='A075', 'DA47'='A026', 'DV02'='A077', ID= 4-8: 'DV03'='A137', 'DV07'='A076', 'DV09'='A046', 'DV10'='A071', ID= 9-13: 'DV12'='A011', 'DV13'='A072', 'DV14'='A025', 'DV17'='A138',
Unless you have calibrated yourself the measurement set you are using for imaging, it is essential to do the following.
# In CASA
tb.open('IRAS16293_Band9.fixed.rebin.ms/POINTING', nomodify = False)
a = tb.rownumbers()
tb.removerows(a)
tb.close()
Continuum Emission Imaging
As shown in Figure 2 of IRAS16293Band9, these observations have a 7 pointing mosaic, which you can visualize with the next command (see also Figure 1).
# Make a plot of the mosaic pattern with field ids identified
au.plotmosaic(vis='IRAS16293_Band9.fixed.rebin.ms',sourceid='0',
doplot=True,figfile='mosaic_pattern.png')
You might wonder which of these pointings or fields have stronger line emission. You can determine this by running the next command. We set spw='1' which contains CO (6-5), and you can check all the plots for all the fields by just hitting the arrow for "next". Also see Figure 2 for an example of the output.
# See which mosaic fields have the brightest line emission
plotms(vis='IRAS16293_Band9.fixed.rebin.ms',
field='',xaxis='uvdist', yaxis='amp',
spw='1', avgchannel='3840',coloraxis='field',
iteraxis='field')
Another way to visualize the uv-amp distribution is with the next plotms. This will load all the fields for pointing into the plot, and you can change interactively between the spws using the appropriate button in plotms. See Figure 3 for a plot of spw 0.
# Iterate over spws all fields
plotms(vis='IRAS16293_Band9.fixed.rebin.ms',
field='',xaxis='uvdist', yaxis='amp',
spw='', avgchannel='3840',coloraxis='field',
iteraxis='spw')
We need to note the channels that do not present line emission. To identify them we can plot the amplitude vs frequency for all the spectral windows. Figure 4 shows the output from the next command for spw 1 (field 3). You are dealing with a set of data which were split so the calibrated data now reside in the "data" column of the new measurement set. Therefore the following call to plotms must specify ydatacolumn = "data". If you had never performed the split, the calibrated data would have been accessed from the "corrected" column of the old, unsplit, measurement set. Those who have used previous versions of CASA should note that in those older versions, plotms would default to using the "data" column if there were no data in the "corrected" column and one specified ydatacolumn="corrected".
# Select Channels for uvcontsub. Average every 8 channels for
# increased S/N.
plotms(vis='IRAS16293_Band9.fixed.rebin.ms',
xaxis='channel', yaxis='amp',field='3',
spw='', avgtime='1e8',avgscan=True,coloraxis='corr',
avgchannel='',iteraxis='spw',ydatacolumn='data',
xselfscale=True,yselfscale=True)
Using the plots from the previous command you can select a block of channels that are free of strong line emission. To do this in the most rigorous way, you would first need to make dirty image cubes and then examine the spectra on and near the continuum sources, because many more weaker lines will become apparent in the cubes. Here we provide you with the channel selection we used for continuum emission.
The next flagmanager command will save the flags state in case you might need to re-do the steps from this point.
flagmanager(vis='IRAS16293_Band9.fixed.rebin.ms',mode='save',versionname='Original')
Now we will flag the channels identified to have line emission, so we will be ready to proceed with the continuum imaging.
# flagging the line channels
flagdata(vis='IRAS16293_Band9.fixed.rebin.ms',
spw='0:1700~2500,1:400~500;1100~1900;3000~3450,2:1700~2200;3100~3839,3:0~800;1600~2900;3300~3839',action='apply', flagbackup=False)
To help clean to be faster, you can average 200 channels, as you do not need to have a high spectral resolution for continuum imaging.
split(vis='IRAS16293_Band9.fixed.rebin.ms',
outputvis='IRAS16293_Band9.fixed.CONT.msAVG',field='',
datacolumn='data',width=200,spw='0~3:20~3819')
Clean and self-calibrate the continuum
Now we proceed with the cleaning itself. For the clean task you will need to choose the right cell size, map size, etc. Below are the parameters that were used in this casaguide. Note that the command will be interactive and you will have to set the cleaning boxes.
# cleaning 1
os.system('rm -rf *CONTIN*')
clean(vis='IRAS16293_Band9.fixed.CONT.msAVG',
imagename='IRAS16293_Band9.fixed.CONTIN',
spw='0~3',field='',phasecenter='3',
mode='mfs',imagermode='mosaic',
imsize=500,cell='0.04arcsec',minpb=0.25,
interactive=True,
weighting='briggs',robust=0.5,
niter=10000, threshold='0.3mJy',usescratch=False)
Now we proceed with gaincal to perform self calibration (for phases) on the data. Even with minsnr=3.0 there will be some failed solutions.
# get one scan per field
os.system('rm -rf cont_pcal1')
gaincal(vis='IRAS16293_Band9.fixed.CONT.msAVG',caltable='cont_pcal1',gaintype='T',
refant='DV14',calmode='p',combine='',spw='',field='',
solint='inf',minsnr=3.0,minblperant=4)
The solutions for five antennas are shown in Figure 6, and are produced with the next plotcal. Be sure to check all the antennas.
plotcal(caltable='cont_pcal1',xaxis='time',yaxis='phase',
spw='',iteration='antenna',subplot=511,plotrange=[0,0,-80,80],
antenna='',timerange='')
# The phases are remarkably stable for Band 9!
We then apply the table to the data to proceed right after with the proper cleaning, to check the improvement of the data.
# application of the table
applycal(vis='IRAS16293_Band9.fixed.CONT.msAVG',spwmap=[],spw='',
gaintable=['cont_pcal1'],calwt=False,flagbackup=False)
Now you are ready for the second cleaning. For this one we are using the mask resulting from the previous cleaning. In the next blocks of commands we will repeat this process (clean, viewer, gaincal, plotcal, applycal, clean, etc.) to improve the overall quality of the image.
# cleaning 2
os.system('rm -rf *CONTIN_P1*')
clean(vis='IRAS16293_Band9.fixed.CONT.msAVG',
imagename='IRAS16293_Band9.fixed.CONTIN_P1',
spw='0~3',field='',phasecenter='3',
mode='mfs',imagermode='mosaic',
imsize=500,cell='0.04arcsec',minpb=0.25,
interactive=True,
weighting='briggs',robust=0.5,
mask='IRAS16293_Band9.fixed.CONTIN.mask',
niter=10000, threshold='0.3mJy', usescratch=False)
# visualize the image
imview (raster=[{'file':'IRAS16293_Band9.fixed.CONTIN.image',
'range': [-0.01,0.5]},
{'file': 'IRAS16293_Band9.fixed.CONTIN_P1.image',
'range': [-0.01,0.5]}])
The image "IRAS16293_Band9.fixed.CONTIN_P1.image" should have 2 very clearly defined regions of emission, moreso then the first clean run "IRAS16293_Band9.fixed.CONTIN.image," following self-calibration.
Second gaincal step for self calibration.
# second self-cal
os.system('rm -rf cont_pcal2')
gaincal(vis='IRAS16293_Band9.fixed.CONT.msAVG',caltable='cont_pcal2',gaintype='T',
refant='DV14',calmode='p',combine='',spw='',field='',
solint='inf',minsnr=3.0,minblperant=4)
#Checking the table
plotcal(caltable='cont_pcal2',xaxis='time',yaxis='phase',
spw='',iteration='antenna',subplot=511,plotrange=[0,0,-80,80],
antenna='',timerange='')
# Applying results to the data
applycal(vis='IRAS16293_Band9.fixed.CONT.msAVG',spwmap=[],spw='',
gaintable=['cont_pcal2'],calwt=False,flagbackup=False)
We continue with the next cleaning.
# Cleaning 3
os.system('rm -rf *CONTIN_P2*')
clean(vis='IRAS16293_Band9.fixed.CONT.msAVG',
imagename='IRAS16293_Band9.fixed.CONTIN_P2',
spw='0~3',field='',phasecenter='3',
mode='mfs',imagermode='mosaic',
imsize=500,cell='0.04arcsec',minpb=0.25,
interactive=True,
weighting='briggs',robust=0.5,
mask='IRAS16293_Band9.fixed.CONTIN.mask',
niter=10000, threshold='0.3mJy', usescratch=False)
# to view the resulting image
imview (raster=[{'file':'IRAS16293_Band9.fixed.CONTIN_P1.image',
'range': [-0.01,0.5]},
{'file': 'IRAS16293_Band9.fixed.CONTIN_P2.image',
'range': [-0.01,0.5]}])
# There is very little change, the bold could try another iteration using a shorter
# solint phase self-cal. For now, moving on to amplitude.
Next iteration with gaincal, now using amplitude and phase for "calmode".
# For amplitude self-cal we only want to get one solution per spw, per dataset
os.system('rm -rf cont_apcal')
gaincal(vis='IRAS16293_Band9.fixed.CONT.msAVG',caltable='cont_apcal',gaintype='T',
refant='DV14',calmode='ap',combine='scan,field',spw='',field='',
solint='8000s',minsnr=3.0,minblperant=4,gaintable='cont_pcal2')
# plotting the table
plotcal(caltable='cont_apcal',xaxis='time',yaxis='amp',
spw='',iteration='antenna',subplot=511,plotrange=[0,0,0.5,2],
antenna='',timerange='')
# and applying to the data
applycal(vis='IRAS16293_Band9.fixed.CONT.msAVG',spwmap=[],spw='',
gaintable=['cont_pcal2','cont_apcal'],calwt=False,flagbackup=False)
# A large correction is present for the 2nd dataset. This amplitude
# error was causing many of the artifacts noticeable in the previous
# images.
Continuing the final cleaning step.
# Cleaning 4
os.system('rm -rf *CONTIN_AP*')
clean(vis='IRAS16293_Band9.fixed.CONT.msAVG',
imagename='IRAS16293_Band9.fixed.CONTIN_AP',
spw='0~3',field='',phasecenter='3',
mode='mfs',imagermode='mosaic',
imsize=500,cell='0.04arcsec',minpb=0.25,
interactive=True,
weighting='briggs',robust=0.5,
mask='IRAS16293_Band9.fixed.CONTIN_P2.mask',
niter=10000, threshold='0.3mJy', usescratch=False)
# and to view the image
imview (raster=[{'file':'IRAS16293_Band9.fixed.CONTIN_P2.image',
'range': [-0.01,0.5]},
{'file': 'IRAS16293_Band9.fixed.CONTIN_AP.image',
'range': [-0.01,0.5]}])
We use now impbcor to produce an image corrected by the primary beam response. We are outputting a sub-region of the cube because the area covered is large compared to the area of interest, and we want to keep the file size of the result as small as possible.
# correcting for primary beam
impbcor(imagename='IRAS16293_Band9.fixed.CONTIN_AP.image',
pbimage='IRAS16293_Band9.fixed.CONTIN_AP.flux',
outfile='IRAS16293_Band9.fixed.CONTIN_AP.image.pbcor.subim',
box='150,104,345,299')
# and checking the image
imview (raster=[{'file': 'IRAS16293_Band9.fixed.CONTIN_AP.image.pbcor.subim',
'range': [-0.01,1.5],'colorwedge':True}],
out='IRAS16293_Band9.fixed.CONTIN_AP.image.pbcor.subim.png')
In Figure 7 we show the final image of the self-calibration process. Source A is the weaker source to the SE and Source B is the stronger source to the NW. The rms noise is about 18 mJy/beam, and the restoring beam is 0.29 x 0.17 arcsec, PA=-71 deg with robust=0.5 weighting. The peak intensities for the two components are about 3.19 (B) and 1.07 (A) Jy/beam, while the flux densities are about 11.9 (B) and 9 (A) Jy.
Now you can export your final image in FITS format, so you can analyze it in other software packages if desired.
exportfits(imagename='IRAS16293_Band9.fixed.CONTIN_AP.image.pbcor.subim',
fitsimage='IRAS16293_Band9.fixed.CONTIN_AP.image.pbcor.subim.fits')
Spectral Line Imaging
# Restore the flags made for the continuum imaging
flagmanager(vis='IRAS16293_Band9.fixed.rebin.ms',mode='restore',versionname='Original')
Then we subtract the continuum, using the line free channels
# Subtract the continuum from the line data
contspw='0:20~1700;2500~3820,1:20~400;500~1100;1900~3000;3450~3800,2:20~1700;2200~3100,3:800~1600;2900~3300'
uvcontsub(vis='IRAS16293_Band9.fixed.rebin.ms',fitspw=contspw,fitorder=1,combine='spw')
Before we commence with the cleaning we also need to apply the tables we produced in the self calibration steps.
# Apply final self-calibration to the line data
applycal(vis='IRAS16293_Band9.fixed.rebin.ms.contsub',spwmap=[],spw='',
gaintable=['cont_pcal2','cont_apcal'],calwt=False,flagbackup=False)
Some notes for this sections: the threshold was chosen after doing a couple of smaller test images. This is about 3 sigma in channels with weak line emission. The mask was also chosen from a few test images.
Create a velocity cube for H13CN (J=8-7)
In this section we provide an example for how to create a Doppler corrected velocity cube for one of the interesting lines in the data: H13CN (J=8-7).
A good reference for finding rest frequencies is http://www.splatalogue.net/. Another good reference for Band 9 specifically is the 450 micron CSO survey of Orion: http://iopscience.iop.org/0067-0049/132/2/281/
First, determine the velocity range for H13CN 690.55207 GHz with the help of the next plotms
plotms(vis='IRAS16293_Band9.fixed.rebin.ms.contsub',
xaxis='velocity', yaxis='amp',field='',
spw='1', avgtime='1e8',avgscan=True,coloraxis='field',
transform=True,freqframe='LSRK',restfreq='690.55207GHz',
avgchannel='8',
ydatacolumn='data',plotfile='Vel_H13CN.png')
Now do an interactive clean, making clean boxes for the emission. Stop cleaning when noise like residuals remain. Number of iterations between showing each interactive window can be controlled within the interactive viewer. We used about 400 iterations.
os.system('rm -rf IRAS16293_Band9.fixed.H13CN_8_7*')
clean(vis='IRAS16293_Band9.fixed.rebin.ms.contsub',
imagename='IRAS16293_Band9.fixed.H13CN_8_7',
spw='1',field='',phasecenter='3',
mode='velocity',outframe='LSRK',restfreq='690.55207GHz',
nchan=60,start='-10km/s',width='0.5km/s',
imagermode='mosaic',
imsize=500,cell='0.04arcsec',minpb=0.25,
interactive=True,
weighting='briggs',robust=0.5,
niter=500, threshold='220.0mJy')
You need to correct for the primary beam attenuation and subimage the cubes to a smaller more manageable size since most of the mosaic is empty.
impbcor(imagename='IRAS16293_Band9.fixed.H13CN_8_7.image',
pbimage='IRAS16293_Band9.fixed.H13CN_8_7.flux',
outfile='IRAS16293_Band9.fixed.H13CN_8_7.image.pbcor.subim',
box='150,104,345,299')
To enable analysis in other software packages, you can export the data from CASA to FITS format using exportfits
exportfits(imagename='IRAS16293_Band9.fixed.H13CN_8_7.image.pbcor.subim',
fitsimage='IRAS16293_Band9.fixed.H13CN_8_7.image.pbcor.subim.fits')
Moment Maps for H13CN (J=8-7)
To proceed with the moment maps generation, you will need to determine the channel range and pixel range for inclusion for higher order moments.
viewer('IRAS16293_Band9.fixed.H13CN_8_7.image.pbcor.subim')
Generation of moment 0
os.system('rm -rf IRAS16293_Band9.fixed.H13CN_8_7.image.pbcor.subim.mom0')
immoments(imagename='IRAS16293_Band9.fixed.H13CN_8_7.image.pbcor.subim',moments=[0],
outfile='IRAS16293_Band9.fixed.H13CN_8_7.image.pbcor.subim.mom0',
chans='10~44')
which you can see in Figure 8. To produce that figure run the following command
imview( raster=[ {'file':'IRAS16293_Band9.fixed.H13CN_8_7.image.pbcor.subim.mom0',
'range': [-1.,5.5]} ],
contour={'file':'IRAS16293_Band9.fixed.CONTIN_AP.image.pbcor.subim',
'base':0, 'unit':0.025, 'levels':[3]},
out='IRAS16293_Band9.fixed.H13CN_8_7.mom0.png')
One cannot make higher order moments that simultaneously include both positive and negative emission. So we separate them below. Also, we need to restrict the higher order moments to regions of real signal. Below we demonstrate using the clean mask for this purpose. We must first subimage it to match the primary beam corrected image. We also set the included pixels to about 5sigma.
immath(imagename='IRAS16293_Band9.fixed.H13CN_8_7.mask',
mode='evalexpr',expr='IM0',box='150,104,345,299',
outfile='IRAS16293_Band9.fixed.H13CN_8_7.mask.subim')
Producing moment 1 and 2
os.system('rm -rf IRAS16293_Band9.fixed.H13CN_8_7.image.pbcor.subim.Positive*')
immoments(imagename='IRAS16293_Band9.fixed.H13CN_8_7.image.pbcor.subim',moments=[1,2],
outfile='IRAS16293_Band9.fixed.H13CN_8_7.image.pbcor.subim.Positive',
mask='IRAS16293_Band9.fixed.H13CN_8_7.mask.subim',
chans='10~44',includepix=[0.4,100])
os.system('rm -rf IRAS16293_Band9.fixed.H13CN_8_7.image.pbcor.subim.Negative*')
immoments(imagename='IRAS16293_Band9.fixed.H13CN_8_7.image.pbcor.subim',moments=[1,2],
outfile='IRAS16293_Band9.fixed.H13CN_8_7.image.pbcor.subim.Negative',
mask='IRAS16293_Band9.fixed.H13CN_8_7.mask.subim',
chans='10~44',includepix=[-100,-0.4])
Visualizing the maps
imview( raster=[ {'file':'IRAS16293_Band9.fixed.H13CN_8_7.image.pbcor.subim.Negative.weighted_coord',
'colorwedge':True}],
contour={'file':'IRAS16293_Band9.fixed.CONTIN_AP.image.pbcor.subim',
'base':0, 'unit':0.025, 'levels':[3]},
out='IRAS16293_Band9.fixed.H13CN_8_7.Negative_mom1.png')
imview( raster=[ {'file':'IRAS16293_Band9.fixed.H13CN_8_7.image.pbcor.subim.Positive.weighted_coord',
'colorwedge':True}],
contour={'file':'IRAS16293_Band9.fixed.CONTIN_AP.image.pbcor.subim',
'base':0, 'unit':0.025, 'levels':[3]},
out='IRAS16293_Band9.fixed.H13CN_8_7.Positive_mom1.png')
And finally export the result to FITS format
os.system('rm -rf IRAS16293_Band9.fixed.H13CN_8_7.image.pbcor.subim.fits')
exportfits(imagename='IRAS16293_Band9.fixed.H13CN_8_7.image.pbcor.subim',
fitsimage='IRAS16293_Band9.fixed.H13CN_8_7.image.pbcor.subim.fits')
os.system('rm -rf IRAS16293_Band9.fixed.H13CN_8_7.image.pbcor.subim.mom0.fits')
exportfits(imagename='IRAS16293_Band9.fixed.H13CN_8_7.image.pbcor.subim.mom0',
fitsimage='IRAS16293_Band9.fixed.H13CN_8_7.image.pbcor.subim.mom0.fits')
os.system('rm -rf IRAS16293_Band9.fixed.H13CN_8_7.image.pbcor.subim.Negative.weighted_coord.fits')
exportfits(imagename='IRAS16293_Band9.fixed.H13CN_8_7.image.pbcor.subim.Negative.weighted_coord',
fitsimage='IRAS16293_Band9.fixed.H13CN_8_7.image.pbcor.subim.Negative.weighted_coord.fits')
os.system('rm -rf IRAS16293_Band9.fixed.H13CN_8_7.image.pbcor.subim.Positive.weighted_coord.fits')
exportfits(imagename='IRAS16293_Band9.fixed.H13CN_8_7.image.pbcor.subim.Positive.weighted_coord',
fitsimage='IRAS16293_Band9.fixed.H13CN_8_7.image.pbcor.subim.Positive.weighted_coord.fits')
os.system('rm -rf IRAS16293_Band9.fixed.H13CN_8_7.image.pbcor.subim.Negative.weighted_dispersion_coord.fits')
exportfits(imagename='IRAS16293_Band9.fixed.H13CN_8_7.image.pbcor.subim.Negative.weighted_dispersion_coord',
fitsimage='IRAS16293_Band9.fixed.H13CN_8_7.image.pbcor.subim.Negative.weighted_dispersion_coord.fits')
os.system('rm -rf IRAS16293_Band9.fixed.H13CN_8_7.image.pbcor.subim.Positive.weighted_dispersion_coord.fits')
exportfits(imagename='IRAS16293_Band9.fixed.H13CN_8_7.image.pbcor.subim.Positive.weighted_dispersion_coord',
fitsimage='IRAS16293_Band9.fixed.H13CN_8_7.image.pbcor.subim.Positive.weighted_dispersion_coord.fits')
Last checked on CASA Version 4.5.0.