Open in Colab: https://colab.research.google.com/github/casangi/examples/blob/master/community/simulation_script_demo.ipynb
Simulation in CASA
Original Author: rurvashi@aoc.nrao.edu
Description
Get creative with data sets to be used for test scripts and characterization of numerical features/changes. This notebook goes beneath the simobserve task and illustrates simple ways in which developers and test writers can make full use of the flexibility offered by our tools and the imager framework. It also exercises some usage modes that our users regularly encounter and exposes some quirks of our scripting interface(s). Rudimentary image and data display routines are included below.
Topics Covered below
- Install CASA 6 and Import required libraries
- Make an empty MS with the desired sub-structure
- Make a true sky model
- Predict visibilities onto the DATA column of the MS
- Add noise and other errors
- A few example use cases
- Image one channel
- Cube imaging with a spectral line
- Continuum wideband imaging with model subtraction
- Self-calibration and imaging
- Ideas for CASA developers and test writers to do beyond these examples.
Installation
[ ]:
import os
print("installing casa...")
os.system('pip install casaconfig')
os.system('pip install casatools==6.6.4.34')
os.system('pip install casatasks==6.6.4.34')
print("complete")
[ ]:
print("configuring casa...")
## For google colab
mydatapath = '/content/.casa/data'
## For a local install (choose any writeable path for upto 1GB of metadata)
#mydatapath = '<PATHNAME>/data'
import pathlib
from casaconfig import config
if not pathlib.Path(mydatapath).exists():
pathlib.Path(mydatapath).mkdir(parents=True)
config.measurespath=mydatapath
print("complete")
Import Libraries
[ ]:
print("install external dependencies")
os.system('pip install astropy')
## For google colab (py 3.10)
os.system('pip install numpy==1.24.4')
print("complete")
[80]:
# Import required tools/tasks
from casatools import simulator, image, table, coordsys, measures, componentlist, quanta, ctsys, vpmanager
from casatasks import tclean, ft, imhead, listobs, exportfits, flagdata, bandpass, applycal
from casatasks.private import simutil
import os
import pylab as pl
import numpy as np
from astropy.io import fits
from astropy.wcs import WCS
# Instantiate all the required tools
sm = simulator()
ia = image()
tb = table()
cs = coordsys()
me = measures()
qa = quanta()
cl = componentlist()
vp = vpmanager()
qa = quanta()
mysu = simutil.simutil()
Make an empty MS with the desired uvw/scan/field/ddid setup
Construct an empty Measurement Set that has the desired observation setup. This includes antenna configuration, phase center direction, spectral windows, date and timerange of the observation, structure of scans/spws/obsidd/fieldids (and all other MS metadata). Evaluate UVW coordinates for the entire observation and initialize the DATA column to zero.
Methods
Make an empty MS
[2]:
def makeMSFrame(msname = 'sim_data.ms'):
"""
Construct an empty Measurement Set that has the desired observation setup.
"""
os.system('rm -rf '+msname)
## Open the simulator
sm.open(ms=msname);
## Read/create an antenna configuration.
## Canned antenna config text files are located here : /home/casa/data/trunk/alma/simmos/*cfg
antennalist = os.path.join( ctsys.resolve("alma/simmos") ,"vla.d.cfg")
## Read x,y,z locations, antenna diameter and name, and telescope name/position.
## Fictitious telescopes can be simulated by specifying x, y, z, d, an, telname, antpos.
## x,y,z are locations in meters in ITRF (Earth centered) coordinates.
## d, an are lists of antenna diameter and name.
## telname and obspos are the name and coordinates of the observatory.
(x,y,z,d,an,an2,telname, obspos) = mysu.readantenna(antennalist)
obsposname=telname
## Set the antenna configuration
sm.setconfig(telescopename=telname,
x=x,
y=y,
z=z,
dishdiameter=d,
mount=['alt-az'],
antname=an,
coordsystem='global',
referencelocation=me.observatory(obsposname));
## Set the polarization mode (this goes to the FEED subtable)
sm.setfeed(mode='perfect R L', pol=['']);
## Set the spectral window and polarization (one data-description-id).
## Call multiple times with different names for multiple SPWs or pol setups.
sm.setspwindow(spwname="LBand",
freq='1.0GHz',
deltafreq='0.1GHz',
freqresolution='0.2GHz',
nchannels=10,
stokes='RR LL');
## Setup source/field information (i.e. where the observation phase center is)
## Call multiple times for different pointings or source locations.
sm.setfield( sourcename="fake",
sourcedirection=me.direction(rf='J2000', v0='19h59m28.5s',v1='+40d44m01.5s'));
## Set shadow/elevation limits (if you care). These set flags.
sm.setlimits(shadowlimit=0.01, elevationlimit='1deg');
## Leave autocorrelations out of the MS.
sm.setauto(autocorrwt=0.0);
## Set the integration time, and the convention to use for timerange specification
## Note : It is convenient to pick the hourangle mode as all times specified in sm.observe()
## will be relative to when the source transits.
sm.settimes(integrationtime='2000s',
usehourangle=True,
referencetime=me.epoch('UTC','2019/10/4/00:00:00'));
## Construct MS metadata and UVW values for one scan and ddid
## Call multiple times for multiple scans.
## Call this with different sourcenames (fields) and spw/pol settings as defined above.
## Timesteps will be defined in intervals of 'integrationtime', between starttime and stoptime.
sm.observe(sourcename="fake",
spwname='LBand',
starttime='-5.0h',
stoptime='+5.0h');
## Close the simulator
sm.close()
## Unflag everything (unless you care about elevation/shadow flags)
flagdata(vis=msname,mode='unflag')
Plot columns of the MS
[3]:
def plotData(msname='sim_data.ms', myplot='uv'):
"""
Options : myplot='uv'
myplot='data_spectrum'
"""
from matplotlib.collections import LineCollection
tb.open(msname)
# UV coverage plot
if myplot=='uv':
pl.figure(figsize=(4,4))
pl.clf()
uvw = tb.getcol('UVW')
pl.plot( uvw[0], uvw[1], '.')
pl.plot( -uvw[0], -uvw[1], '.')
pl.title('UV Coverage')
# Spectrum of chosen column. Make a linecollection out of each row in the MS.
if myplot=='data_spectrum' or myplot=='corr_spectrum' or myplot=='resdata_spectrum' or myplot=='rescorr_spectrum' or myplot=='model_spectrum':
dats=None
if myplot=='data_spectrum':
dats = tb.getcol('DATA')
if myplot=='corr_spectrum':
dats = tb.getcol('CORRECTED_DATA')
if myplot=='resdata_spectrum':
dats = tb.getcol('DATA') - tb.getcol('MODEL_DATA')
if myplot=='rescorr_spectrum':
dats = tb.getcol('CORRECTED_DATA') - tb.getcol('MODEL_DATA')
if myplot=='model_spectrum':
dats = tb.getcol('MODEL_DATA')
xs = np.zeros((dats.shape[2],dats.shape[1]),'int')
for chan in range(0,dats.shape[1]):
xs[:,chan] = chan
npl = dats.shape[0]
fig, ax = pl.subplots(1,npl,figsize=(10,4))
for pol in range(0,dats.shape[0]):
x = xs
y = np.abs(dats[pol,:,:]).T
data = np.stack(( x,y ), axis=2)
ax[pol].add_collection(LineCollection(data))
ax[pol].set_title(myplot + ' \n pol '+str(pol))
ax[pol].set_xlim(x.min(), x.max())
ax[pol].set_ylim(y.min(), y.max())
pl.show()
tb.close()
Examples
Make a Measurement Set and inspect it
[4]:
makeMSFrame()
[5]:
plotData(myplot='uv')
[6]:
listobs(vis='sim_data.ms', listfile='obslist.txt', verbose=False, overwrite=True)
## print(os.popen('obslist.txt').read()) # ?permission denied?
fp = open('obslist.txt')
for aline in fp.readlines():
print(aline.replace('\n',''))
fp.close()
================================================================================
MeasurementSet Name: /home/vega/rurvashi/TestCASA/Test2024/DemoBooks/examples/community/sim_data.ms MS Version 2
================================================================================
Observer: CASA simulator Project: CASA simulation
Observation: VLA(27 antennas)
Data records: 6318 Total elapsed time = 36000 seconds
Observed from 03-Oct-2019/21:21:40.2 to 04-Oct-2019/07:21:40.2 (UTC)
Fields: 1
ID Code Name RA Decl Epoch SrcId nRows
0 fake 19:59:28.500000 +40.44.01.50000 J2000 0 6318
Spectral Windows: (1 unique spectral windows and 1 unique polarization setups)
SpwID Name #Chans Frame Ch0(MHz) ChanWid(kHz) TotBW(kHz) CtrFreq(MHz) Corrs
0 LBand 10 TOPO 1000.000 100000.000 1000000.0 1450.0000 RR LL
Antennas: 27 'name'='station'
ID= 0-5: 'W01'='P', 'W02'='P', 'W03'='P', 'W04'='P', 'W05'='P', 'W06'='P',
ID= 6-11: 'W07'='P', 'W08'='P', 'W09'='P', 'E01'='P', 'E02'='P', 'E03'='P',
ID= 12-17: 'E04'='P', 'E05'='P', 'E06'='P', 'E07'='P', 'E08'='P', 'E09'='P',
ID= 18-23: 'N01'='P', 'N02'='P', 'N03'='P', 'N04'='P', 'N05'='P', 'N06'='P',
ID= 24-26: 'N07'='P', 'N08'='P', 'N09'='P'
Make a True Sky Model (component list and/or image)
Construct a true sky model for which visibilities will be simulated and stored in the DATA column. This could be a component list (with real-world positions and point or gaussian component types), or a CASA image with a real-world coordinate system and pixels containing model sky values. It is possible to also evaluate component lists onto CASA images.
Methods
Make a source list
Once made,it can be used either for direction evaluation of simulated visibilities, or first evaluated onto a CASA image before visibility prediction.
[7]:
def makeCompList(clname_true='sim_onepoint.cl',
dirn='J2000 19h59m28.5s +40d44m01.5s',
flux=5.0,
index=-1.0):
# Make sure the cl doesn't already exist. The tool will complain otherwise.
os.system('rm -rf '+clname_true)
cl.done()
# Add sources, one at a time.
# Call multiple times to add multiple sources. ( Change the 'dir', obviously )
cl.addcomponent(dir=dirn,
flux=flux, # For a gaussian, this is the integrated area.
fluxunit='Jy',
freq='1.5GHz',
shape='point', ## Point source
# shape='gaussian', ## Gaussian
# majoraxis="5.0arcmin",
# minoraxis='2.0arcmin',
spectrumtype="spectral index",
index=index)
# Print out the contents of the componentlist
#print('Contents of the component list')
#print(cl.torecord())
# Save the file
cl.rename(filename=clname_true)
cl.done()
Make an empty CASA image
[8]:
def makeEmptyImage(imname_true='sim_onepoint_true.im'):
## Define the center of the image
radir = '19h59m28.5s'
decdir = '+40d44m01.5s'
## Make the image from a shape
ia.close()
ia.fromshape(imname_true,[256,256,1,10],overwrite=True)
## Make a coordinate system
cs=ia.coordsys()
cs.setunits(['rad','rad','','Hz'])
cell_rad=qa.convert(qa.quantity('8.0arcsec'),"rad")['value']
cs.setincrement([-cell_rad,cell_rad],'direction')
cs.setreferencevalue([qa.convert(radir,'rad')['value'],qa.convert(decdir,'rad')['value']],type="direction")
cs.setreferencevalue('1.0GHz','spectral')
cs.setreferencepixel([0],'spectral')
cs.setincrement('0.1GHz','spectral')
## Set the coordinate system in the image
ia.setcoordsys(cs.torecord())
ia.setbrightnessunit("Jy/pixel")
ia.set(0.0)
ia.close()
### Note : If there is an error in this step, subsequent steps will give errors of " Invalid Table Operation : SetupNewTable.... imagename is already opened (is in the table cache)"
## The only way out of this is to restart the kernel (equivalent to exit and restart CASA).
## Any other way ?
Evaluate the component list onto the image cube
[9]:
def evalCompList(clname='sim_onepoint.cl', imname='sim_onepoint_true.im'):
## Evaluate a component list
cl.open(clname)
ia.open(imname)
ia.modify(cl.torecord(),subtract=False)
ia.close()
cl.done()
Edit pixel values directly
[10]:
def editPixels(imname='sim_onepoint_true.im'):
## Edit pixel values directly
ia.open(imname)
pix = ia.getchunk()
shp = ia.shape()
#pix.fill(0.0)
#pix[ int(shp[0]/2), int(shp[1]/2), 0, :] = 4.0 # A flat spectrum unpolarized source of amplitude 1 Jy and located at the center of the image.
pix[ int(shp[0]/2), int(shp[1]/2), 0, 6] = pix[ int(shp[0]/2), int(shp[1]/2), 0, 6] + 2.0 # Add a spectral line in channel 1
ia.putchunk( pix )
ia.close()
View an Image Cube
Use some image viewer, or just pull the pixels out and use matplotlib
[11]:
# Display an image using AstroPy, with coordinate system rendering.
def dispAstropy(imname='sim_onepoint_true.im'):
exportfits(imagename=imname, fitsimage=imname+'.fits', overwrite=True)
hdu = fits.open(imname+'.fits')[0]
wcs = WCS(hdu.header,naxis=2)
fig = pl.figure()
fig.add_subplot(121, projection=wcs)
pl.imshow(hdu.data[0,0,:,:], origin='lower', cmap=pl.cm.viridis)
pl.xlabel('RA')
pl.ylabel('Dec')
# Display an image cube or a single plane image.
# For a Cube, show the image at chan 0 and a spectrum at the location of the peak in chan0.
# For a Single plane image, show the image.
def dispImage(imname='sim_onepoint_true.im', useAstropy=False):
ia.open(imname)
pix = ia.getchunk()
shp = ia.shape()
ia.close()
pl.figure(figsize=(10,4))
pl.clf()
if shp[3]>1:
pl.subplot(121)
if useAstropy==False:
pl.imshow(pix[:,:,0,0])
pl.title('Image from channel 0')
else:
dispAstropy(imname)
if shp[3]>1:
pl.subplot(122)
ploc = np.where( pix == pix.max() )
pl.plot(pix[ploc[0][0], ploc[1][0],0,:])
pl.title('Spectrum at source peak')
pl.xlabel('Channel')
Examples
Make a component list and evaluate it onto a CASA image
[12]:
## Make the component list
makeCompList()
## Make an empty CASA image
makeEmptyImage()
## Evaluate the component list onto the CASA image
evalCompList()
## Display
dispImage()
[13]:
## Edit the pixels of the CASA image directly (e.g. add a spectral line)
editPixels()
## Display
dispImage()
Simulate visibilities from the sky model into the DATA column of the MS
Simulate visibilities for the true sky model, applying a variety of instrumental effects. This step either evaluates the DFT of a component model, or uses an imaging (de)gridder. Instrumental effects can be applied either by pre-processing the sky model before ‘standard’ degridding, or by invoking one of the wide-field imaging gridders to apply W-term, A-term and mosaicing effects. Noise, extra spectral lines or RFI may be added at this point, as well as gain errors via the application of carefully constructed calibration tables.
Methods
Use the simulator tool
Visibilities are predicted and saved in the DATA column of the MS. It is preferable to use the simulator only when the standard gridder is desired. Prediction can be done from an input model image or a component list
[14]:
def predictSim(msname='sim_data.ms',
imname='sim_onepoint_true.im',
clname='sim_onepoint.cl',
usemod='im',
usepb=False,
vptable=''):
"""
usemod = 'im' : use the imname image
usemod = 'cl' : use the clname component list
usepb = True : to include static primary beams in the simulation.
"""
## Open an existing MS Frame
sm.openfromms(msname)
# Include primary Beams
if usepb==True:
if vptable=='':
sm.setvp( dovp = True, usedefaultvp = True )
else:
sm.setvp( dovp = True, usedefaultvp = False, vptable=vptable )
if usemod=='im':
# Predict from a model image
sm.predict( imagename = imname, incremental=False)
else:
# Predict from a component list
sm.predict( complist = clname ,incremental=False)
# Close the tool
sm.close()
Use imager (or ft)
Visibilities are predicted and saved in the MODEL_DATA column of the MS. The values must then be copied to the DATA column. Use this approach when non-standard gridders are required, typically when instrument-dependent effects are included, or when Taylor-coefficient wideband image models are to be used for visibility prediction.
Step 1 : Simulate visibilities into the MODEL column using tclean
tclean can be used for model prediction with all gridders (‘standard’, ‘wproject’, ‘mosaic’, ‘awproject’). Wide-field and full-beam effects along with parallactic angle rotation may be included with appropriate settings. tclean can predict model visibilities only from input images and not component lists.
[15]:
## Use an input model sky image - widefield gridders
def predictImager(msname='sim_data.ms',
imname_true='sim_onepoint_true.im',
gridder='standard'):
os.system('rm -rf sim_predict.*')
# Run tclean in predictModel mode.
tclean(vis=msname,
startmodel=imname_true,
imagename='sim_predict',
savemodel='modelcolumn',
imsize=256,
cell='8.0arcsec',
specmode='cube',
interpolation='nearest',
start='1.0GHz',
width='0.1GHz',
nchan=10,
reffreq='1.5Hz',
gridder=gridder,
normtype='flatsky', # sky model is flat-sky
cfcache='sim_predict.cfcache',
wbawp=True, # ensure that gridders='mosaic' and 'awproject' do freq-dep PBs
pblimit=0.05,
conjbeams=False,
calcres=False,
calcpsf=True,
niter=0,
wprojplanes=1)
Step 1 : Simulate visibilities into the MODEL column using ft
The ‘ft’ task implements the equivalent of gridder=’standard’ in tclean. Wide-field effects cannot be simulated.
In addition, it offers the ability to predict visibilities from component lists (which tclean does not).
[16]:
def predictFt(msname='sim_data.ms',
imname='sim_onepoint_true.im',
clname='sim_onepoint.cl',
usemod='im'):
if usemod=='im':
## Use an image name and the ft task
ft(vis = msname, model = imname, incremental = False, usescratch=True)
else:
## Use a component list and the ft task
ft(vis = msname, complist = clname, incremental = False, usescratch=True)
Step 2 : Copy contents of the MODEL column to the DATA column
[17]:
### Copy visibilities from the MODEL column to the data columns
### This is required when predicting using tclean or ft as they will only write to the MODEL column
def copyModelToData(msname='sim_data.ms'):
tb.open(msname,nomodify=False);
moddata = tb.getcol(columnname='MODEL_DATA');
tb.putcol(columnname='DATA',value=moddata);
#tb.putcol(columnname='CORRECTED_DATA',value=moddata);
moddata.fill(0.0);
tb.putcol(columnname='MODEL_DATA',value=moddata);
tb.close();
Examples
If the above commands were run in order, the component list contains only a steep-spectrum continuum source, but the model image cube contains an additional spectral line in it.
Option 1 : Predict using the simulator and a componentlist
[18]:
# Predict Visibilities
predictSim(usemod='cl')
# Plot
plotData(myplot='data_spectrum')
Option 2 : Predict using the simulator and an input image
[19]:
# Predict visibilities
predictSim(usemod='im')
# Plot
plotData(myplot='data_spectrum')
Option 3 : Predict using tclean and a model image with gridder=’standard’
[20]:
predictImager()
copyModelToData()
plotData(myplot='data_spectrum')
Option 4 : Predict using ft and a component list
[21]:
# Predict using ft
predictFt(usemod='cl')
copyModelToData()
# Plot
plotData(myplot='data_spectrum')
Option 5 : Predict using ft and an input image
[22]:
# Predict using ft
predictFt(usemod='im')
copyModelToData()
# Plot
plotData(myplot='data_spectrum')
Add Noise and other errors to the simulated visibilities
Methods
Add Visibility noise
[23]:
## Add Gaussian random noise
def addNoiseSim(msname='sim_data.ms'):
sm.openfromms(msname);
sm.setseed(50)
sm.setnoise(mode='simplenoise',simplenoise='0.05Jy');
sm.corrupt();
sm.close();
Add random numbers
[24]:
def addNoiseRand(msname = 'sim_data.ms'):
## Add noise and other variations
tb.open( msname, nomodify=False )
dat = tb.getcol('DATA')
## Add noise to the first few channels only. ( Ideally, add separately to real and imag parts... )
from numpy import random
dat[:,0:4,:] = dat[:,0:4,:] + 0.5 * random.random( dat[:,0:4,:].shape )
## Add some RFI in a few rows and channels....
#dat[ :, :, 1 ] = dat[ :, :, 1] + 2.0
tb.putcol( 'DATA', dat )
tb.close()
Add antenna gain errors
[48]:
## Add antenna gain errors.
def addGainErrors(msname='sim_data.ms'):
sm.openfromms(msname);
sm.setseed(50)
sm.setgain(mode='fbm',amplitude=0.1)
sm.corrupt()
sm.close();
## Note : This step sometimes produces NaN/Inf in the visibilities and plotData() will complain ! If so, just run it again. I thought that setting the seed will control this, but apparently not.
Examples
Use the simulator to add Gaussian random noise (1 Jy rms noise)
[26]:
addNoiseSim()
plotData(myplot='data_spectrum')
A few Imaging and Calibration examples
Image one channel
[27]:
# Call tclean
os.system('rm -rf try0.*')
tclean(vis='sim_data.ms',
imagename='try0',
datacolumn='data',
spw='0:5', # pick channel 5 and image it
imsize=300,
cell='8.0arcsec',
specmode='mfs',
gridder='standard',
niter=200,
gain=0.3,
interactive=False,
usemask='auto-multithresh');
[28]:
# Display the output restored image
dispImage('try0.image')
[29]:
# Display the point spread function
dispImage('try0.psf')
Cube Imaging of a spectral-line dataset
This is a spectral line dataset with noise
[30]:
# Call tclean
os.system('rm -rf try1.*')
tclean(vis='sim_data.ms',
imagename='try1',
datacolumn='data',
imsize=300,
cell='8.0arcsec',
specmode='cube',
interpolation='nearest',
gridder='standard',
niter=200,
gain=0.3,
savemodel='modelcolumn');
[31]:
# Display the output restored image
dispImage('try1.image')
[32]:
# Display the residual image
dispImage('try1.residual')
Continuum imaging with model subtraction
Pick the line-free channels (all but chan 6) and fit a 2nd order polynomial to the spectrum.
[33]:
# Call tclean
os.system('rm -rf try2.*')
tclean(vis='sim_data.ms',
imagename='try2',
datacolumn='data',
spw='0:0~5,0:7~9', # Select line-free channels
imsize=300,
cell='8.0arcsec',
specmode='mfs',
deconvolver='mtmfs',
nterms=3,
gridder='standard',
niter=150,
gain=0.3);
[34]:
# Display the output restored image
dispImage('try2.image.tt0')
[35]:
# Predict the tclean mtmfs model onto all channels.
tclean(vis='sim_data.ms',
imagename='try2',
datacolumn='data',
spw='', # Select all channels to predict onto.
imsize=300,
cell='8.0arcsec',
specmode='mfs',
deconvolver='mtmfs',
nterms=3,
gridder='standard',
niter=0,
calcres=False,
calcpsf=False,
savemodel='modelcolumn');
[36]:
# Plot residual data.
plotData(myplot='resdata_spectrum')
This shows the continuum power law emission subtracted out and only the spectral line remaining in the data (if the model is subtracted).
If the ‘uvsub’ task is run, this is what would get saved in corrected_data. It is also a form of continuum modeling and subtraction.
Imaging with Gain Errors and Self Calibration
First, re-simulate by starting from ideal visibilities, and adding gain errors and noise.
[61]:
# Predict visibilities
predictSim(usemod='im')
[62]:
# Simulate antenna gain errors
addGainErrors()
# Add noise on top of the gain-corrupted data
addNoiseSim()
# Display
plotData(myplot='data_spectrum')
Image the corrupted data
[63]:
# Call tclean
os.system('rm -rf try3.*')
tclean(vis='sim_data.ms',
imagename='try3',
datacolumn='data',
imsize=300,
cell='8.0arcsec',
specmode='cube',
interpolation='nearest',
gridder='standard',
niter=150, # Don't go too deep since the data are corrupted
gain=0.3,
mask='circle[[150pix,150pix],3pix]', # Give it a mask to help. Without this, the self-cal isn't as good.
savemodel='modelcolumn');
[64]:
# Display the output restored image
dispImage('try3.image')
[65]:
# Display the new residual image
dispImage('try3.residual')
This image shows artifacts from gain errors (different from the pure noise-like errors in the previous simulation)
Calculate gain solutions (since we have already saved the model)
[66]:
bandpass(vis='sim_data.ms',
caltable='sc.tab',solint='int');
Apply gain solutions
[67]:
applycal(vis='sim_data.ms',
gaintable='sc.tab')
[68]:
## Plot Calibrated data
plotData(myplot='corr_spectrum')
Compare with the above uncalibrated data. Also, compare with visibilities simulated just with noise. Subsequent imaging should use the corrected_data column.
[69]:
# Call tclean to image the corrected data
os.system('rm -rf try4.*')
tclean(vis='sim_data.ms',
imagename='try4',
datacolumn='corrected',
imsize=300,
cell='8.0arcsec',
specmode='cube',
interpolation='nearest',
gridder='standard',
niter=200, # Go deeper now. Also, no mask needed
gain=0.3,
savemodel='modelcolumn');
[70]:
# Display the output restored image
dispImage('try4.image')
[71]:
# Display the residual image
dispImage('try4.residual')
Use custom primary beams for simulation
This section illustrates how to use the vpmanager tool to create primary beam models that may be used during simulation (of multi-channel datasets).
Start with a flat spectrum sky model : To clearly assess the effect of a primary beam, let us start with a flat spectrum 1 Jy source in the sky model. The source is located away from the pointing center, to see the effect of the primary beam.
[72]:
## Make the component list
## Note : The 1Jy flat spectrum source is located slightly offset
## from the phase/pointing center
makeCompList(clname_true='sim_onepoint_offset.cl',
dirn='J2000 19h59m28.5s +40d50m01.5s', #+40d44m01.5s',
flux=1.0,index=0.0)
## Make an empty CASA image
makeEmptyImage(imname_true='sim_onepoint_offset_true.im')
## Evaluate the component list onto the CASA image
evalCompList(clname='sim_onepoint_offset.cl', imname='sim_onepoint_offset_true.im')
## Display
dispImage(imname='sim_onepoint_offset_true.im')
The above is a flat spectrum source. When visibilities are simulated with a primary beam, the PB amplitude (and frequency dependence) should be evident in a similar plot below.
Make an MS with the desired telescope name
[73]:
makeMSFrame(msname='sim_custompb.ms')
[74]:
listobs(vis='sim_custompb.ms', listfile='obslist_custom.txt', verbose=False, overwrite=True)
fp = open('obslist_custom.txt')
for aline in fp.readlines():
print(aline.replace('\n',''))
fp.close()
## Check the telescope name below.
## This is the name that must be used to make a VPTable
================================================================================
MeasurementSet Name: /home/vega/rurvashi/TestCASA/Test2024/DemoBooks/examples/community/sim_custompb.ms MS Version 2
================================================================================
Observer: CASA simulator Project: CASA simulation
Observation: VLA(27 antennas)
Data records: 6318 Total elapsed time = 36000 seconds
Observed from 03-Oct-2019/21:21:40.2 to 04-Oct-2019/07:21:40.2 (UTC)
Fields: 1
ID Code Name RA Decl Epoch SrcId nRows
0 fake 19:59:28.500000 +40.44.01.50000 J2000 0 6318
Spectral Windows: (1 unique spectral windows and 1 unique polarization setups)
SpwID Name #Chans Frame Ch0(MHz) ChanWid(kHz) TotBW(kHz) CtrFreq(MHz) Corrs
0 LBand 10 TOPO 1000.000 100000.000 1000000.0 1450.0000 RR LL
Antennas: 27 'name'='station'
ID= 0-5: 'W01'='P', 'W02'='P', 'W03'='P', 'W04'='P', 'W05'='P', 'W06'='P',
ID= 6-11: 'W07'='P', 'W08'='P', 'W09'='P', 'E01'='P', 'E02'='P', 'E03'='P',
ID= 12-17: 'E04'='P', 'E05'='P', 'E06'='P', 'E07'='P', 'E08'='P', 'E09'='P',
ID= 18-23: 'N01'='P', 'N02'='P', 'N03'='P', 'N04'='P', 'N05'='P', 'N06'='P',
ID= 24-26: 'N07'='P', 'N08'='P', 'N09'='P'
Make a VP table First, we create a vp table using polynomial models. We will call this ‘VLA’ to match with the telescope name (telname) used in makeMSFrame().
[88]:
def make_vp_table_poly(vptab='my_beams.vp'):
os.system('rm -rf '+vptab)
vp.reset()
# frequency polynomial to be given here
vp.setpbpoly(telescope ='VLA', ## Use whatever telescope name is in makeMSFrame()
coeff=np.array([1, -2.608e-3, 27.357e-7, -13.091e-10, 2.368e-13]))
# Construct an airy disk from a dish diameter
#vp.setpbairy(telescope='VLA',
# dishdiam=[qa.quantity('45m')])
# Construct a Gaussian PB
#vp.setpbgauss(telescope='VLA',
# halfwidth=[qa.quantity('0.2deg')])
vp.summarizevps()
vp.saveastable(vptab)
[89]:
make_vp_table_poly(vptab='my_beams.vp')
2024-07-19 21:24:07 WARN PBMath::pbMathInterfaceForCommonPB ATCA_L3 not yet implemented defaulting to L2 version
Predict visibilities with the vptable model
[90]:
predictSim(msname='sim_custompb.ms',
imname='sim_onepoint_offset_true.im',
clname='sim_onepoint_offset.cl',
usemod='im', ### or 'cl'
usepb=True,
vptable='my_beams.vp')
plotData(msname='sim_custompb.ms', myplot='data_spectrum')
The spectral slope seen above is from the frequency dependence of the primary beam used in the simulation. The source is at approximately the 0.8 gain level of the primary beam at the middle frequency of 1.5 GHz.
The same custom PBs may be used during imaging by supplying the vptable as an input parameter to tclean (with its ‘standard’, ‘mosaic’ and ‘wproject’ gridder options).
[ ]:
Predict visibilities with the default PB model
For ‘known’ observatories, canned PB models may also be used. The following shows how to use the canned VLA PB model. Note the difference from the custom model used above in the simulation.
[87]:
predictSim(msname='sim_custompb.ms',
imname='sim_onepoint_offset_true.im',
clname='sim_onepoint_offset.cl',
usemod='im', ### or 'cl'
usepb=True) ### Let it pick the default PB model for the telname
plotData(msname='sim_custompb.ms', myplot='data_spectrum')
[ ]:
[ ]: