Tohban EOVSA Imaging Tutorial A-Z
Revision as of 22:19, 11 November 2021 by Sshaik (talk | contribs) (→Step 2: Concatenate all the 10 mins data, if there are any)
Connection details to pipeline server
Step 1: Downloading raw data (IDB) on Pipeline machine
On CASA of the Pipeline machine,
If you are not using mobaxterm, directly sshing to pipeline through linux or mac computer, to run casa start tcsh.
from astropy.time import Time
import os
trange = Time(['2017-08-21 20:15:00', '2017-08-21 20:25:00'])
#### (Optional) change output path, default current directory "./" #####
outpath = './msdata/'
if not os.path.exists(outpath):
os.makedirs(outpath)
######################################################
msfiles = importeovsa(idbfiles=trange, ncpu=1, visprefix=outpath)
OR (for a time range longer than 10 minutes)
from suncasa.tasks import task_calibeovsa as calibeovsa
from suncasa.tasks import task_importeovsa as timporteovsa
from split_cli import split_cli as split
import dump_tsys as dt
from util import Time
import numpy as np
import os
from glob import glob
from eovsapy import util
trange = Time(['2017-08-21 20:15:00', '2017-08-21 20:35:00'])
idbdir = util.get_idbdir(trange[0])
info = dt.rd_fdb(trange[0])
for k, v in info.iteritems():
info[k] = info[k][~(info[k] == '')]
sidx = np.where(
np.logical_and(info['SOURCEID'] == 'Sun', info['PROJECTID'] == 'NormalObserving') & np.logical_and(
info['ST_TS'].astype(np.float) >= trange[0].lv,
info['ST_TS'].astype(np.float) <= trange[
1].lv))
filelist = info['FILE'][sidx]
outpath = './msdata/'
if not os.path.exists(outpath):
os.makedirs(outpath)
inpath = idbdir + '{}/'.format(trange[0].datetime.strftime("%Y%m%d"))
ncpu = 1
msfiles = timporteovsa.importeovsa(idbfiles=[inpath + ll for ll in filelist], ncpu=ncpu, timebin="0s", width=1,
visprefix=outpath,
nocreatms=False, doconcat=False,
modelms="", doscaling=False, keep_nsclms=False, udb_corr=True)
If you come across errors with calibeovsa, add following lines to your ~/.casa/init.py file.
import sys
sys.path.append('/common/python')
sys.path.append('/common/python/packages/pipeline_casa')
execfile('/common/python/suncasa/tasks/mytasks.py')
Step 2: Concatenate all the 10 mins data, if there are any
Follow this step if there are more than one .ms files, if not run step 3 directly.
# This is to set the path/name for the concatenated files concatvis = os.path.basename(msfiles[0])[:11] + '_concat.ms' vis = calibeovsa(msfiles, doconcat=True, concatvis=concatvis, msoutdir=outpath)
Step 3: Calibration
This will calibrate the input visibility, write out calibration tables under /data1/eovsa/caltable/, and apply the calibration. If doimage=True, a quicklook image will be produced (by integrating over the entire time).
calibeovsa(vis='IDB20170821202020.ms', caltype=['refpha','phacal'], doimage=True)
Connection details to Inti server
For Windows, on Mobaxterm,
- Connect to any afsconnect servers (for example, afsaccess3.njit.edu) using UCID ss2273 and your UCID password.
- ssh -X UCID@inti.hpcnet.campus.njit.edu
cd /inti/data/users/YOURDIRECTORY #YOURDIRECTORY folder was created by Sijie for having enough space with EOVSA data analysis
For Linux or Mac machine,
Transferring details between servers
Software details on the servers
On Inti, when logging in for the first time, please add the following lines to your accounts .bashrc file.
>>vi .bashrc Insert the text given below and save it.
#### setting start ####
if [ $HOSTNAME == "baozi.hpcnet.campus.njit.edu" ]; then
source /srg/.setenv_baozi
fi
if [ $HOSTNAME == "inti.hpcnet.campus.njit.edu" ]; then
source /inti/.setenv_inti
fi
if [ $HOSTNAME == "guko.resource.campus.njit.edu" ]; then
source /data/data/.setenv_guko
fi
#### setting end ####
Both CASA 5 and 6 are available on Inti.
Please enter the bash environment on inti, and load the desired casa environment with the alias below. To load CASA 5: Enter bash environment by giving >>bash
>> loadcasa5 >> suncasa #This should load the software and you are ready for analysis
To load CASA 6: Enter bash environment by giving >>bash
>>loadcasa6 >>ipython #This should load the software and you are ready for analysis
Here, for example, to use clean, first start ipython as given above, then type in >>from casatasks import tclean
Step 4: Self-calibration on Inti server
https://github.com/binchensun/casa-eovsa/blob/master/slfcal_example.py
from suncasa.utils import helioimage2fits as hf
import os
import numpy as np
import pickle
from matplotlib import gridspec as gridspec
from sunpy import map as smap
from matplotlib import pyplot as plt
import time
# =========== task handlers =============
dofullsun=0 # initial full-sun imaging
domasks=0 # get masks
doslfcal=0 # main cycle of doing selfcalibration
doapply=1 # apply the results
doclean_slfcaled=1 # perform clean for self-calibrated data
# ============ declaring the working directories ============
workdir = os.getcwd()+'/' #main working directory. Using current directory in this example
slfcaldir = workdir+'slfcal/' #place to put all selfcalibration products
imagedir = slfcaldir+'images/' #place to put all selfcalibration images
maskdir = slfcaldir+'masks/' #place to put clean masks
imagedir_slfcaled = slfcaldir+'images_slfcaled/' #place to put final self-calibrated images
caltbdir = slfcaldir+'caltbs/' # place to put calibration tables
# make these directories if they do not already exist
dirs = [workdir, slfcaldir, imagedir, maskdir, imagedir_slfcaled, caltbdir]
for d in dirs:
if not os.path.exists(d):
os.makedirs(d)
# ============ Split a short time for self-calibration ===========
# input visibility
ms_in = workdir + 'IDB20170821201800-202300.4s.corrected.ms'
# output, selfcaled, visibility
ms_slfcaled = workdir + os.path.basename(ms_in).replace('corrected','slfcaled')
# intermediate small visibility for selfcalbration
# selected time range for generating self-calibration solutions
trange='2017/08/21/20:21:10~2017/08/21/20:21:30'
slfcalms = slfcaldir+'slfcalms.XX.slfcal'
slfcaledms = slfcaldir+'slfcalms.XX.slfcaled'
if not os.path.exists(slfcalms):
split(vis=ms_in,outputvis=slfcalms,datacolumn='data',timerange=trange,correlation='XX')
# ============ Prior definitions for spectral windows, antennas, pixel numbers =========
spws=[str(s+1) for s in range(30)]
antennas='0~12'
npix=512
nround=3 #number of slfcal cycles
# parameters specific to the event (found from step 1)
phasecenter='J2000 10h02m59 11d58m14'
xran=[280,480]
yran=[-50,150]
# =========== Step 1, doing a full-Sun image to find out phasecenter and appropriate field of view =========
if dofullsun:
#initial mfs clean to find out the image phase center
im_init='fullsun_init'
os.system('rm -rf '+im_init+'*')
tclean(vis=slfcalms,
antenna='0~12',
imagename=im_init,
spw='1~15',
specmode='mfs',
timerange=trange,
imsize=[npix],
cell=['5arcsec'],
niter=1000,
gain=0.05,
stokes='I',
restoringbeam=['30arcsec'],
interactive=False,
pbcor=True)
hf.imreg(vis=slfcalms,imagefile=im_init+'.image.pbcor',fitsfile=im_init+'.fits',
timerange=trange,usephacenter=False,verbose=True)
clnjunks = ['.flux', '.mask', '.model', '.psf', '.residual','.sumwt','.pb','.image']
for clnjunk in clnjunks:
if os.path.exists(im_init + clnjunk):
os.system('rm -rf '+im_init + clnjunk)
from sunpy import map as smap
from matplotlib import pyplot as plt
fig = plt.figure(figsize=(6,6))
ax = fig.add_subplot(111)
eomap=smap.Map(im_init+'.fits')
#eomap.data=eomap.data.reshape((npix,npix))
eomap.plot_settings['cmap'] = plt.get_cmap('jet')
eomap.plot(axes = ax)
eomap.draw_limb()
plt.show()
viewer(im_init+'.image.pbcor')
# =========== Step 2 (optional), generate masks =========
# if skipped, will not use any masks
if domasks:
clearcal(slfcalms)
delmod(slfcalms)
antennas='0~12'
pol='XX'
imgprefix=maskdir+'slf_t0'
# step 1: set up the clean masks
img_init=imgprefix+'_init_ar_'
os.system('rm -rf '+img_init+'*')
#spwrans_mask=['1~5','6~12','13~20','21~30']
spwrans_mask=['1~12']
#convert to a list of spws
spwrans_mask_list = [[str(i) for i in (np.arange(int(m.split('~')[0]),int(m.split('~')[1])))] for m in spwrans_mask]
masks=[]
imnames=[]
for spwran in spwrans_mask:
imname=img_init+spwran.replace('~','-')
try:
tclean(vis=slfcalms,
antenna='0~12',
imagename=imname,
spw=spwran,
specmode='mfs',
timerange=trange,
imsize=[npix],
cell=['2arcsec'],
niter=1000,
gain=0.05,
stokes='XX',
restoringbeam=['20arcsec'],
phasecenter=phasecenter,
weighting='briggs',
robust=1.0,
interactive=True,
datacolumn='data',
pbcor=True,
savemodel='modelcolumn')
imnames.append(imname+'.image')
masks.append(imname+'.mask')
clnjunks = ['.flux', '.model', '.psf', '.residual']
for clnjunk in clnjunks:
if os.path.exists(imname + clnjunk):
os.system('rm -rf '+ imname + clnjunk)
except:
print('error in cleaning spw: '+spwran)
pickle.dump(masks,open(slfcaldir+'masks.p','wb'))
# =========== Step 3, main step of selfcalibration =========
if doslfcal:
if os.path.exists(slfcaldir+'masks.p'):
masks=pickle.load(open(slfcaldir+'masks.p','rb'))
if not os.path.exists(slfcaldir+'masks.p'):
print 'masks do not exist. Use default mask'
masks=[]
os.system('rm -rf '+imagedir+'*')
os.system('rm -rf '+caltbdir+'*')
#first step: make a mock caltable for the entire database
print('Processing ' + trange)
slftbs=[]
calprefix=caltbdir+'slf'
imgprefix=imagedir+'slf'
tb.open(slfcalms+'/SPECTRAL_WINDOW')
reffreqs=tb.getcol('REF_FREQUENCY')
bdwds=tb.getcol('TOTAL_BANDWIDTH')
cfreqs=reffreqs+bdwds/2.
tb.close()
# starting beam size at 3.4 GHz in arcsec
sbeam=40.
strtmp=[m.replace(':','') for m in trange.split('~')]
timestr='t'+strtmp[0]+'-'+strtmp[1]
refantenna='0'
# number of iterations for each round
niters=[100, 300, 500]
# roburst value for weighting the baselines
robusts=[1.0, 0.5, 0.0]
# apply calibration tables? Set to true for most cases
doapplycal=[1,1,1]
# modes for calibration, 'p' for phase-only, 'a' for amplitude only, 'ap' for both
calmodes=['p','p','a']
# setting uvranges for model image (optional, not used here)
uvranges=['','','']
for n in range(nround):
slfcal_tb_g= calprefix+'.G'+str(n)
fig = plt.figure(figsize=(8.4,7.))
gs = gridspec.GridSpec(5, 6)
for s,sp in enumerate(spws):
print 'processing spw: '+sp
cfreq=cfreqs[int(sp)]
# setting restoring beam size (not very useful for selfcal anyway, but just to see the results)
bm=max(sbeam*cfreqs[1]/cfreq, 6.)
slfcal_img = imgprefix+'.spw'+sp.zfill(2)+'.slfcal'+str(n)
# only the first round uses nearby spws for getting initial model
if n == 0:
spbg=max(int(sp)-2,1)
sped=min(int(sp)+2,30)
spwran=str(spbg)+'~'+str(sped)
print('using spw {0:s} as model'.format(spwran))
if 'spwrans_mask' in vars():
for m, spwran_mask in enumerate(spwrans_mask):
if sp in spwran_mask:
mask = masks[m]
print('using mask {0:s}'.format(mask))
findmask = True
if not findmask:
print('mask not found. Do use any masks')
else:
spwran = sp
if 'spwrans_mask' in vars():
for m, spwran_mask in enumerate(spwrans_mask):
if sp in spwran_mask:
mask = masks[m]
print 'using mask {0:s}'.format(mask)
findmask = True
if not findmask:
print('mask not found. Do use any masks')
try:
tclean(vis=slfcalms,
antenna=antennas,
imagename=slfcal_img,
uvrange=uvranges[n],
spw=spwran,
specmode='mfs',
timerange=trange,
imsize=[npix],
cell=['2arcsec'],
niter=niters[n],
gain=0.05,
stokes='XX', #use pol XX image as the model
weighting='briggs',
robust=robusts[n],
phasecenter=phasecenter,
mask=mask,
restoringbeam=[str(bm)+'arcsec'],
pbcor=False,
interactive=False,
savemodel='modelcolumn')
if os.path.exists(slfcal_img+'.image'):
fitsfile=slfcal_img+'.fits'
hf.imreg(vis=slfcalms,imagefile=slfcal_img+'.image',fitsfile=fitsfile,
timerange=trange,usephacenter=False,toTb=True,verbose=False,overwrite=True)
clnjunks = ['.mask','.flux', '.model', '.psf', '.residual', '.image','.pb','.image.pbcor','.sumwt']
for clnjunk in clnjunks:
if os.path.exists(slfcal_img + clnjunk):
os.system('rm -rf '+ slfcal_img + clnjunk)
ax = fig.add_subplot(gs[s])
eomap=smap.Map(fitsfile)
eomap.plot_settings['cmap'] = plt.get_cmap('jet')
eomap.plot(axes = ax)
eomap.draw_limb()
#eomap.draw_grid()
ax.set_title(' ')
ax.get_xaxis().set_visible(False)
ax.get_yaxis().set_visible(False)
ax.set_xlim(xran)
ax.set_ylim(yran)
os.system('rm -f '+ fitsfile)
except:
print 'error in cleaning spw: '+sp
print 'using nearby spws for initial model'
sp_e=int(sp)+2
sp_i=int(sp)-2
if sp_i < 1:
sp_i = 1
if sp_e > 30:
sp_e = 30
sp_=str(sp_i)+'~'+str(sp_e)
try:
tget(tclean)
spw=sp_
print('using spw {0:s} as model'.format(sp_))
tclean()
except:
print 'still not successful. abort...'
break
gaincal(vis=slfcalms, refant=refantenna,antenna=antennas,caltable=slfcal_tb_g,spw=sp, uvrange='',\
gaintable=[],selectdata=True,timerange=trange,solint='inf',gaintype='G',calmode=calmodes[n],\
combine='',minblperant=4,minsnr=2,append=True)
if not os.path.exists(slfcal_tb_g):
print 'No solution found in spw: '+sp
figname=imagedir+'slf_t0_n{:d}.png'.format(n)
plt.subplots_adjust(left=0, bottom=0, right=1, top=1, wspace=0, hspace=0)
plt.savefig(figname)
time.sleep(10)
plt.close()
if os.path.exists(slfcal_tb_g):
slftbs.append(slfcal_tb_g)
slftb=[slfcal_tb_g]
os.chdir(slfcaldir)
if calmodes[n] == 'p':
plotcal(caltable=slfcal_tb_g,antenna='1~12',xaxis='freq',yaxis='phase',\
subplot=431,plotrange=[-1,-1,-180,180],iteration='antenna',figfile=slfcal_tb_g+'.png',showgui=False)
if calmodes[n] == 'a':
plotcal(caltable=slfcal_tb_g,antenna='1~12',xaxis='freq',yaxis='amp',\
subplot=431,plotrange=[-1,-1,0,2.],iteration='antenna',figfile=slfcal_tb_g+'.png',showgui=False)
os.chdir(workdir)
if doapplycal[n]:
clearcal(slfcalms)
delmod(slfcalms)
applycal(vis=slfcalms,gaintable=slftb,spw=','.join(spws),selectdata=True,\
antenna=antennas,interp='nearest',flagbackup=False,applymode='calonly',calwt=False)
if n < nround-1:
prompt=raw_input('Continuing to selfcal?')
#prompt='y'
if prompt.lower() == 'n':
if os.path.exists(slfcaledms):
os.system('rm -rf '+slfcaledms)
split(slfcalms,slfcaledms,datacolumn='corrected')
print 'Final calibrated ms is {0:s}'.format(slfcaledms)
break
if prompt.lower() == 'y':
slfcalms_=slfcalms+str(n)
if os.path.exists(slfcalms_):
os.system('rm -rf '+slfcalms_)
split(slfcalms,slfcalms_,datacolumn='corrected')
slfcalms=slfcalms_
else:
if os.path.exists(slfcaledms):
os.system('rm -rf '+slfcaledms)
split(slfcalms,slfcaledms,datacolumn='corrected')
print 'Final calibrated ms is {0:s}'.format(slfcaledms)
# =========== Step 4: Apply self-calibration tables =========
if doapply:
import glob
os.chdir(workdir)
clearcal(ms_in)
clearcal(slfcalms)
applycal(vis=slfcalms,gaintable=slftbs,spw=','.join(spws),selectdata=True,\
antenna=antennas,interp='linear',flagbackup=False,applymode='calonly',calwt=False)
applycal(vis=ms_in,gaintable=slftbs,spw=','.join(spws),selectdata=True,\
antenna=antennas,interp='linear',flagbackup=False,applymode='calonly',calwt=False)
if os.path.exists(ms_slfcaled):
os.system('rm -rf '+ms_slfcaled)
split(ms_in, ms_slfcaled,datacolumn='corrected')
# =========== Step 5: Generate final self-calibrated images (optional) =========
if doclean_slfcaled:
import glob
pol='XX'
print('Processing ' + trange)
img_final=imagedir_slfcaled+'/slf_final_{0}_t0'.format(pol)
vis = ms_slfcaled
spws=[str(s+1) for s in range(30)]
tb.open(vis+'/SPECTRAL_WINDOW')
reffreqs=tb.getcol('REF_FREQUENCY')
bdwds=tb.getcol('TOTAL_BANDWIDTH')
cfreqs=reffreqs+bdwds/2.
tb.close()
sbeam=30.
from matplotlib import gridspec as gridspec
from sunpy import map as smap
from matplotlib import pyplot as plt
fitsfiles=[]
for s,sp in enumerate(spws):
cfreq=cfreqs[int(sp)]
bm=max(sbeam*cfreqs[1]/cfreq,4.)
imname=img_final+'_s'+sp.zfill(2)
fitsfile=imname+'.fits'
if not os.path.exists(fitsfile):
print 'cleaning spw {0:s} with beam size {1:.1f}"'.format(sp,bm)
try:
tclean(vis=vis,
antenna=antennas,
imagename=imname,
spw=sp,
specmode='mfs',
timerange=trange,
imsize=[npix],
cell=['1arcsec'],
niter=1000,
gain=0.05,
stokes=pol,
weighting='briggs',
robust=2.0,
restoringbeam=[str(bm)+'arcsec'],
phasecenter=phasecenter,
mask='',
pbcor=True,
interactive=False)
except:
print 'cleaning spw '+sp+' unsuccessful. Proceed to next spw'
continue
if os.path.exists(imname+'.image.pbcor'):
imn = imname+'.image.pbcor'
hf.imreg(vis=vis,imagefile=imn,fitsfile=fitsfile,
timerange=trange,usephacenter=False,toTb=True,verbose=False)
fitsfiles.append(fitsfile)
junks=['.flux','.model','.psf','.residual','.mask','.image','.pb','.image.pbcor','.sumwt']
for junk in junks:
if os.path.exists(imname+junk):
os.system('rm -rf '+imname+junk)
else:
print('fits file '+fitsfile+' already exists, skip clean...')
fitsfiles.append(fitsfile)
fig = plt.figure(figsize=(8.4,7.))
gs = gridspec.GridSpec(5, 6)
for s,sp in enumerate(spws):
cfreq=cfreqs[int(sp)]
ax = fig.add_subplot(gs[s])
eomap=smap.Map(fitsfiles[s])
eomap.plot_settings['cmap'] = plt.get_cmap('jet')
eomap.plot(axes = ax)
eomap.draw_limb()
ax.set_title(' ')
ax.get_xaxis().set_visible(False)
ax.get_yaxis().set_visible(False)
ax.set_xlim(xran)
ax.set_ylim(yran)
plt.text(0.98,0.85,'{0:.1f} GHz'.format(cfreq/1e9),transform=ax.transAxes,ha='right',color='w',fontweight='bold')
plt.subplots_adjust(left=0, bottom=0, right=1, top=1, wspace=0, hspace=0)
plt.show()