SIMIO
Create a Template
Step 1: Get the ms file of the observation you want to mimic.
For this tutorial, we will create a template based on Elias 24, one of the bright DSHARP planet-forming disks (Huang et al. 2018). The original ms file can be obtained in the DSHARP database.
Download the continuum of Elias 24 (0.8GB). This file contains all the infortmation of the observation, including antenna, visibilities, times, etc.
Step 2: Create a folder for your template.
Go to you SIMIO directory, and open your templates folder. Create an empty folder with the name of your template. In this case, the template will be called Elias24.
Step 3: Create sub folders and scripts.
Create 3 empty folders, named "images", "msfiles" and "uvtables". Additionally, create 2 empty scripts named "template_casa_handling.py" and "template_info.py". In this example the template is Elias24, therefore the templates will be called: "Elias24_casa_handling.py" and "Elias24_info.py".
Move the ms file downloaded in step 1 into the folder "msfiles", and extract the ms file from the .tgz.
# Source Geometry, Elias 2-24 in Huang et al. (2018), DSHARP II
dRa = 0. # arcsec
dDec = 0. # arcsec
inc = 29. # degrees
pa = 45.7 # degrees
geom_params = [dRa, dDec, inc, pa]
# Distance
parallax = 7.180888047224439 # GAIA EDR3 mas
dist = 1000. / parallax # pc
# Sky info
temp_center_ra = '16h26m24.0782s' # J2000
temp_center_dec = '-24d16m13.884s' # J2000
temp_semimaj = 1.15 # arcsec
# Clean info
temp_robust = 0.0 # From Andrews et al. (2021)
temp_cellsize = '0.003arcsec' # DSHARP Standard
temp_imsize = 2000 # Big enough for noise estimation
Step 4: Template information.
The file "template_info.py", in this case "Elias24_info.py", contains important information to generate the images for this template. You also can copy/paste the info file of another template and replace the numbers with the ones of your template.
Source Geometry: Obtained from Huang et al. (DSHARP II, 2018). You can add a spatial offset if needed, and the inclination and position angle must be input in degrees. The ms file of Elias 27 has a phase center different from the center of the source, but we will change that in the next step. Therefore, we set dRa and dDec to 0.
Distance: We estimate the distance as the inverse of the parallax. We use GAIA EDR3 (2016, 2021) to obtain the parallax (SIMBAD name of the source is Elias 2-24).
Sky info: Here you should write the coordinates of the source center, in J2000 epoch. The coordinate we chose for Elias 24 is peak brightness pixel. The parameter temp_semimaj is the radial extent of the source in arcsec. This parameter is used to create the mask to generate its images.
Clean info: This parameters are needed to generate the template images. The robust parameters controls the weight that is given to the baselines. In this case, we know the original study was made on images generated with a robust parameter of 0.0 (Andrews et al. DSHARP I, 2021). The cellsize should be around 7x to 10x smaller than the beam size, and the image size should be big enough such that we can estimate the noise level with an annulus around the source.
# Append the au functions from other directory.
sys.path.append('../../codes/analysis_scripts/')
import analysisUtils as au
# Import uv-handling functions
execfile('../../codes/simio_ms2ascii.py')
# Import imaging functions
execfile('../../codes/simio_clean.py')
# Info
prefix = 'Elias24'
ms_cont = 'msfiles/'+prefix+'_avg.ms'
listobs('msfiles/Elias24_continuum.ms')
Step 5.1: Open CASA to prepare ms file.
Generally speaking, high angular resolution observations of disks can result in very heavy ms files. In SIMIO, we are mainly interested in checking the observability of your model under the uv-coverage of a template observation, therefore we can reduce the ms files volume by averaging in frequency and time, without considerably modifying the uv-coverage. The goal of this step is to reduce the data volume of the original observation, thus speeding up all the upcoming processes.
Open CASA in your terminal, from the folder /home/path_to_simio/simio/templates/Elias24/ . From this location, we will call the CASA analysis utils, and also the SIMIO scripts to extract uvtables and generate images.
Set the name of the template as prefix, and create the name for the template measurement set. Execute the task listobs on the downloaded ms file.
Step 5.2: Check the ms file
Here I show a small part of the information given by listobs. Here you can see that the visibilities of Elias 24 are stored in the ms file in 16 spectral windows (numbered from 0 to 15). We will average each spectral window to only 1 channel, and we will also average in time to 24s. This is an aggressive averaging, but it is ok for our comparison goal.
# Remove previous auxiliary files
os.system('rm -rf msfiles/Elias24_avg.ms')
os.system('rm -rf msfiles/Elias24_avg_aux1.ms')
os.system('rm -rf msfiles/Elias24_avg_aux2.ms')
os.system('rm -rf msfiles/Elias24_avg_aux3.ms')
# Split the data to average and combine
split(vis='msfiles/Elias24_continuum.ms', \
outputvis='msfiles/Elias24_avg_aux1.ms', \
timebin='30s', \
spw='0,1,3,4,5,6', \
width=8, \
datacolumn='data', \
intent='OBSERVE_TARGET#ON_SOURCE')
split(vis='msfiles/Elias24_continuum.ms', \
outputvis='msfiles/Elias24_avg_aux2.ms', \
timebin='30s', \
spw='11,15', \
datacolumn='data', \
intent='OBSERVE_TARGET#ON_SOURCE')
split(vis='msfiles/Elias24_continuum.ms', \
outputvis='msfiles/Elias24_avg_aux3.ms', \
timebin='30s', \
spw='2,7,8,9,10,12,13,14', \
width=16, \
datacolumn='data', \
intent='OBSERVE_TARGET#ON_SOURCE')
Step 5.3: Split to average
We start by removing previous ms files, just in case. Then, we split the original ms file we downloaded into 3 different ms files. The first ms file contains the spectral windows (spws) with 4 channels per spw, the second ms file contains the spws with 8 channels, and the third with 16 channels. In each split, we combine those channels into 1, and we also average in time.
This step can also be done with mstransform.
# Join together
concat(vis=['msfiles/Elias24_avg_aux1.ms', \
'msfiles/Elias24_avg_aux2.ms', \
'msfiles/Elias24_avg_aux3.ms'], \
concatvis='msfiles/Elias24_avg.ms', \
dirtol = '0.1arcsec', copypointing = False)
# Remove auxiliary files
os.system('rm -rf msfiles/Elias24_avg_aux1.ms')
os.system('rm -rf msfiles/Elias24_avg_aux2.ms')
os.system('rm -rf msfiles/Elias24_avg_aux3.ms')
# Check averaged Elias24
listobs('msfiles/Elias24_avg.ms')
Step 5.4: Concatenate
Now we take those 3 ms files, and we concatenate them into only 1 file. To save memory, we remove the intermediate auxiliary files.
When checking the new file with listobs, you will see that each spw has only 1 channel.
# Move phase center
fixvis(vis=ms_cont, \
outputvis=ms_cont, \
phasecenter='J2000 16h26m24.0782s -24d16m13.884s')
Step 5.5: Shift phase-center
The phase center of this observation is not in the center of the disk (you can confirm this downloading the continuum fits file form the DSHARP data release webpage).
We will move the phase center of the visibilities to match the location of the peak emission of the disk. This should be a good approximation for the disk center.
It is not mandatory to change the phase-center, and you can skip this step in your own templates.
# Mask geometry
inc, pa = 29., 45.7 # Huang et al. (2018), DSHARP II
pa_mask = str(pa)
# Mask semimajor axis size
mask_semimaj = 1.15
semimajor = str(mask_semimaj)+'arcsec'
semiminor = str(mask_semimaj*np.cos(np.deg2rad(inc)))+'arcsec'
# Mask center
center_ra = '16h26m24.0782s'
center_dec = '-24d16m13.884s'
center = center_ra + ' ' + center_dec
# Create masks
mask_LB = 'ellipse[[%s, %s], [%s, %s], %sdeg]' % \
(center_ra, center_dec, \
semimajor, semiminor, pa_mask)
res_mask_LB = 'annulus[[%s, %s], [%s, %s]]' % \
(center_ra, \
center_dec, \
'2.5arcsec', '3.5arcsec')
Step 5.6: Create a mask for your template.
Using the same information you put in the info script, create a mask for your template and to estimate the noise level. This masks are CASA regions.
# Clean
tclean_wrapper(vis='msfiles/'+prefix+'_avg.ms', \
imagename='images/'+prefix+'_im', \
imsize=2000, \
cellsize='0.003arcsec', \
mask=mask_LB, \
scales=[0, 3, 8, 13], \
robust=0.0, \
gain=0.05, \
smallscalebias=0.32, \
cyclefactor=1.75, \
threshold=str(4*2.05e-02)+'mJy', \
savemodel='modelcolumn', \
niter=10000, \
interactive=True)
Step 6.1: Image your template
Use the task tclean_wrapper to image the averaged ms file. Here you have the freedom to set the CLEAN parameters by yourself, but you should use the info script parameters (same cellsize, imsize, robust).
The parameter smallscalebias will force CLEAN to prefer the bigger scales before cleaning with point sources. The gain=0.05 and cyclefactor=1.75 are for more conservative iterations.
We will clean until we reach 4 sigma, such that afterwards we can apply the JvM correction.
Step 6.2: CLEAN
The CLEAN interactive viewer should look like this, with the mask already set. Press the green arrow iteratively, or press the blue arrow and wait until its finished.
Be patient. This step will take a long time.
Step 6.3: Estimate SNR and apply JvM correction
You have succesfully generated the images. Execute the estimate_SNR function and the JvM correction yo your newly created images. This will give you the angular resolution of your images, and the noise information.
Check your JvM corrected image with the CASA viewer task. Here we have changed the brightness scale.
exportfits(imagename='images/'+prefix+'_im.image', \
fitsimage='images/'+prefix+'_im_noJvM.fits', \
history=False, overwrite=True)
exportfits(imagename='images/'+prefix+'_im.JvMcorr.image', \
fitsimage='images/'+prefix+'_im.fits', \
history=False, overwrite=True)
exportfits(imagename='images/'+prefix+'_im.model', \
fitsimage='images/'+prefix+'_im_model.fits', \
history=False, overwrite=True)
exportfits(imagename='images/'+prefix+'_im.residual', \
fitsimage='images/'+prefix+'_im_residual.fits', \
history=False, overwrite=True)
exportfits(imagename='images/'+prefix+'_im.psf', \
fitsimage='images/'+prefix+'_im_psf.fits', \
history=False, overwrite=True)
Step 6.4: Export your images to fits files
Use the task exportfits to export your newly created images to fits files.
Here I show an example to save the JvM corrected image, the non JvM corrected image, the residuals, the CLEAN model, and the PSF.
# Initialize weight_spectrum column
initweights(vis=ms_cont, wtmode='weight', dowtsp=True)
# Extract spw info
os.system('rm -rf uvtables/'+prefix+'_spws*')
ms.open(ms_cont)
spw_info = ms.getspectralwindowinfo()
np.save('uvtables/'+prefix+'_spws', spw_info.keys())
ms.close()
os.system('rm -rf msfiles/'+prefix+'_cont_spw*')
os.system('rm -rf uvtables/'+prefix+'_cont_spw*')
for i in spw_info.keys():
# Names
ms_file = 'msfiles/' + prefix + '_cont_spw'+str(i)+'.ms'
ascii_file = 'uvtables/' + prefix + '_cont_spw'+str(i)+'.txt'
# Delete any previous spw file
os.system('rm -rf ' + ms_file)
os.system('rm -rf ' + ascii_file)
# Split
split(vis=ms_cont, outputvis=ms_file, \
keepflags=False, datacolumn='data', spw=i)
ms_to_ascii(ms_file, ascii_file, with_flags=True)
# Delete temporal msfile
os.system('rm -rf ' + ms_file)
# Combine uvtable from each spw into a single file
os.system('rm -rf uvtables/'+prefix+'_uvtable.txt')
!ls -1tr uvtables/Elias24_cont_spw*.txt | xargs cat >> uvtables/Elias24_uvtable.txt
# Check weight
os.system('du -h uvtables/'+prefix+'_uvtable.txt')
# Delete previous uvtables
for i in spw_info.keys():
ascii_file = 'uvtables/' + prefix + '_cont_spw'+str(i)+'.txt'
os.system('rm -rf ' + ascii_file)
Step 7: Write the uvtables in .txt
Execute the following piece of code to extract the uvtable of your ms file. You only need to change 1 line of this code: The name of the template, which I highlighted in bold.
The uvtable size will be printed after being concatenated.
# 93M uvtables/Elias24_uvtable.txt
Step 8: That's it!
That's all! your template is now ready to be used as a SIMIO template. The images of the template are stored in the images folder, the averaged ms file is in msfiles, and the visibilities are un the uvtables folder, in txt format.
You can now delete all the .log and .last files that were created by CASA.