View source for Fit an arbitrary sky (image) model to an existing MS


'''Simulating and Fitting a Sky Model to an Existing Measurement Set in the UV Domain'''


''version 20may2015, written with CASA v4.3.1'' 

Versioning note (Jan 2022): In CASA 6, sm is imported as smtool in the standard NRAO CASA installation and simutil is imported as 'from casatasks.private.simutil import simutil'

This casaguide illustrates how to fit an arbitrary sky image model to an existing MS. Along the way it illustrates how to "get your hands" on the raw visibility data in the python interface and how to compute the Chi-squared of the visibility data to your given model. This could be useful if you have a set of source models you would like to compare in the UV domain, but the models are too complicated to be handled by the set of models provided in CASA's uvmodelfit() task (note also that the Nordic ARC has developed a more sophisticated UV model fitting tool set called UVMULTIFIT, documented at [http://www.nordic-alma.se/support/software-tools] and [https://arxiv.org/abs/1401.4984]). The methods illustrated here can be used to estimate the <b>scaling</b> and <b>goodness-of-fit</b> of your model to the MS (although for absolute goodness of fit please note the caveat at the end of this page!). They could also be the starting point for a more comprehensive and sophisticated exploration of a set of possible models using the Chi-squared statistic, for instance. The focus of this guide is an ALMA dataset but the techniques are generally applicable.  

The example assumes we have an ALMA 7m MS; we are interested in the continuum data only; and we have a sky model as a FITS image which we wish to compare with the UV data (i.e., the MS). As written, this tutorial requires you to process one spectral window (SPW) at a time. '''It also assumes each SPW has been averaged down to only a single channel''' through appropriate use of the CASA routine split().

First suppose we have a calibrated ALMA 7m MS called <b>calibrated_final_cont.ms</b> and a sky model image in a FITS file called <b>myskymap.fits</b>. We need to import it into CASA as a CASA image. Furthermore, this image must have <b>four</b> coordinate axes: two spatial coordinates, a FREQUENCY axis, and a STOKES axis.  This is true even if the data and simulation are single-polarization continuum data, in which case the FREQUENCY and STOKES axes can be "singleton" axes. We assume our FITS input model image does not have these extra axes.  There are two ways you can add them: 1) the CASA importfits() task will fill these in for you if you set <b>defaultaxes=False</b> with the appropriate sub-parameters; 2) the simutil() module will allow you to do it. We illustrate the second case:

<pre>
default importfits
fitsimage='myskymap.fits'
imagename='myskymap.image'
overwrite=True
defaultaxes=False
inp
go

from simutil import simutil
util=simutil()

# Set the INCENTER and INWIDTH for the SPW of interest from the listobs()
#   specifically INCENTER +/- 0.5xINWIDTH should encompass all channels in the SPW of interest.
util.modifymodel(inimage='myskymap.image',
                 outimage='myskymapFixed.image',
                 inbright='',indirection='',incell='',incenter='85.015425GHz',
                 inwidth='1937.5MHz',innchan=1)
</pre>
In the call to util.modifymodel() we have taken advantage of the fact that modifymodel() does not change axes for which the corresponding function argument is a null string. This is handy, we don't wish to have to redundantly specify the pixel size and image center ("incell" and "indirection") if they are already correct.


Next we split out the SPW and field(s) of interest --- in this example there happens to be only a single field, but the method would work with a mosaic--- and do some time averaging to reduce the data volume (the averaging by default will be limited by scan and field ID). Two copies of the MS are made; one with the split data, the other will be overwritten to contain the calculated visibilities for the model. Note that in principle you can also do this data selection step using the simulator toolkit's setdata() method, but for the time being the recommended method is to create a fully split MS that contains only the data you want to simulate in a given call to the simulator (i.e. to sm.predict()).  The simulated visibilities will go in the DATA column of the MS.
<pre>
# split out one SPW and the field(s) of interest-
#  make two copies, one of which will get overwritten by the MOCK
spwid='0'
realmsname='calibrated_final_cont_spw'+spwid+'_AVG.ms'
split(vis='calibrated_final_cont.ms',outputvis=realmsname,field='4',
      spw=spwid,combine='scan',timebin='1e8',keepflags=False)
fakemsname='calibrated_final_cont_spw'+spwid+'_AVG_MOCK.ms'
split(vis='calibrated_final_cont.ms',outputvis=fakemsname,field='4',
      spw=spwid,combine='scan',timebin='1e8',keepflags=False)
</pre>
NOTE (added March 2022): you need to be careful with flagged data - the keep flags argument above will eliminate any completely flagged rows (all pols/channels), such as the autocorrelations; but further steps are needed if there are partially flagged rows in the MS.

Currently there is a mismatch in the way real ALMA data are written and how the CASA simulator identifies antennas, with the result that we need to manually set the telescope name in the MS <b>for 7m data (only)</b>: for most other telescopes this will not be the case.  We do this using the table toolkit, which exists in your casa environment as "tb"-- 
<pre>
#
# set the telescope name manually- the need for this is probably a defect -
#  (sim tool infers primary beam from telescope name)
# ***DO THIS ONLY FOR ALMA 7M data!***
#
tb.open(fakemsname+'/OBSERVATION',nomodify=False)
tel=tb.getcol("TELESCOPE_NAME")
for i in range(len(tel)):
  tel[i]="ACA"

tb.putcol("TELESCOPE_NAME",tel)
tb.done()
</pre>

Also, on some occasions ALMA data may have a corrupted POINTING table within the measurement set. This can be fixed by deleting all rows in the POINTING table as follows
<pre>
tb.open(fakemsname+'/POINTING', nomodify=False)                     
tb.removerows(list(range(tb.nrows())))
tb.done()
</pre>
Note that the explicit "list()" in the removerows() call above is for python3, which uses lazy evaluation of range(), which is problematic for toolkit calls in casa6+. One symptom of the problem is zero predicted visibility values, but it is safe in general to delete the POINTING table. (the pointing table contains where each antenna was pointed as a function of time through the course of the observation -- it is important for OTF or full-stokes mosaicking, but not for simple single field or mosaic imaging cases: in those cases the antennas are all assumed to be pointed at the phase center)

Next we use the simulator toolkit--- which is instantiated in your CASA environment as "sm"--- to generate predicted model visibilities from the model image. Only three calls to sim are needed: one to open the MS; one to set the primary beam (which is done from the telescope name); and one to generate the predicted visibilities.
<pre>
sm.openfromms(fakemsname)

# you could do data selection using the sm.setdata() method call here,
#  but we have already created a split() MS by SPW and FIELD so don't need to.
#  using SPLIT() is at present the recommended method.

# this sets the primary beam ("voltage pattern")-
sm.setvp()
# check sanity-
sm.summary()

sm.predict(imagename='lyszmapjyperpixRightSignFixed.image')
# Note 1: you can also include an optional componentlist argument.
# Note 2: the result has gone in the data column, not the model column

# check sanity again  -
plotms(vis=fakemsname,xaxis='uvwave',yaxis='amp',ydatacolumn='data',avgtime='1e8')

# plot what the observed model would look like, sans noise (dirty map) - 
clean(vis=fakemsname,niter=0,cell='2.5arcsec',imsize=[128,128],imagename='model7m')
</pre>
Note that the primary beam CASA by default uses for ALMA antennas is a rough approximation to the real (known) primary beams: it uses the beam that a uniformly illuminated antenna about 15% smaller would have (an Airy disk).

Now we wish to extract the complex visibilities and their associated weights into python to see how well they match up (i.e., compute Chi-squared) and estimate a scale factor for the model from the data. We do this using the MS tool, which is instantiated in your CASA environment as "ms". For some applications, e.g. high-dynamic range imaging over large regions, this could be a limiting factor in the fit quality requiring special steps beyond the scope of this CASA guide (setting a custom primary beam).

<pre>
ms.open(fakemsname)
# Set up the ms tool to select the visibilities we want
ms.selectinit(reset=True)
# generally for ALMA data, datadescid is equal to the SPW number.
#  if in doubt, ms.getspectralwindowinfo() will return a python dictionary
#  showing the properties of the spectral windows and how to map them to datadescid.
spw=1
ms.selectinit(datadescid=spw)
# ms.getdata() returns a dictionary, each element of which is a numpy array
mydat=ms.getdata(['real','imaginary','u','v','weight'])
uu=mydat['u']
vv=mydat['v']
rr=mydat['real']
ii=mydat['imaginary']
# The weights for the model and the data should be identical
#   and equal to 1/error^2 for each visibility.
wt=mydat['weight']
ms.close()

# now get visibility values from the actual (not mock) data-
ms.open(realmsname)
ms.selectinit()
realdat=ms.getdata(['real','imaginary','u','v','weight'])
r2=realdat['real']
i2=realdat['imaginary']
ms.close()
</pre>


With the model and data visibility values in numpy arrays, it is then relatively straightforward to compare them, estimating scaling factors and errors in parameters of interest, etc. Some analysis steps and manipulations are illustrated here --
<pre>
# what are the dims of the data?
r2.shape
# will be [npol,nrows] if no spectral channels 
#   or [npol, nchan, nrows] if n_spectral_channels > 1

# built-in sqrt method does not work elementwise so import numpy's version -
import numpy as np
uvrad=np.sqrt(uu**2+vv**2)
ampl=np.sqrt(rr**2+ii**2)
a2=np.sqrt(r2**2+i2**2)

# first index indexes polarization (there are two --XX, YY); second is null (probably SPW or IF, nominally)-
pl.plot(uvrad,ampl[0,0,:],'.')
# the above plot should match up with the PLOTMS() plot we made just after sm.predict()

# calculate chi2-
c2= (r2[:,0,:]-rr[:,0,:])**2 * wt + (i2[:,0,:]-ii[:,0,:])**2 * wt
# typical data point contributes 1.6 to chi2-
np.median(c2)
# total chi2 per DOF-
c2.sum()/c2.size
# = 2.31 for this dataset
#  ***NOTE that the absolute normalization of the weights may or may not be correct, so take the 
#      chi-squared values with a grain of salt.

# simple single-parameter linear least squares fit of the model to the data-
numerator =  r2[:,0,:] * rr[:,0,:] * wt + i2[:,0,:]*ii[:,0,:] * wt
denominator =  (rr[:,0,:])**2 * wt + (ii[:,0,:])**2 * wt
# linear least-squares scale factor between data and model -
scalefac= numerator.sum() / denominator.sum()
# uncertainty in scalefac, if the weights are correct in an absolute sense -
escalefac=sqrt(1.0/denominator.sum())
</pre>

<b>*A caveat</b>: in current (4.2.2+) versions of CASA, ALMA data that are calibrated by the automated calibration pipeline have weights which have units of 1/Jy^2 (as they should) and values that are close to correct in an absolute sense. For reasons still be determined the absolute accuracy of these weights is about 20%, i.e. they can routinely be 20% off, depending on the correlator configuration and other variables.  Additionally, if no online channel averaging was done during the observation, adjacent spectral channels will be highly correlated due to the lag-domain smoothing which is *always* done by the correlator prior to computing the channelized spectra. This lag-domain -- e.g. Hanning -- smoothing is very common for FX correlators such as the ALMA "Baseline" (12-m array) correlator and serves to reduce spectral ringing due to the finite number of lags measured.  An easy way to deal with the channel-to-channel correlations produced by Hanning smoothing is to spectrally average your data by a factor of N >> 2, giving visibilities which are close to statistically independent.  <b> Please note also</b> that if you are working with manually calibrated data, or data from a telescope other than ALMA, the specific details of the online smoothing and absolute accuracy of the weights may be different, and you will need to understand them to trust the absolute weights if you cannot reliably estimate them from the data or fit yourself.