Create a Simulated Image

This topic is closely related to Create a Component List for Selfcal, but it is given a separate category in hopes that it will also help people find this more general topic.

There are at least two reasons that one would want to create a simulated image. 1. One wants to start with a particular sky model for self-calibration. In this case, the topic Create a Component List for Selfcal may suffice. 2. One wants to perform statistical tests on an image. In the particular case that motivated this document, I was attempting to assess whether I could believe weak sources within an image. In my case, CLEAN bias was not an issue, but there are instances in which it may be.

Assume that there is an image mysrc.image and the corresponding sky model mysrc.model, that has been imaged from the corresponding Measurement Set mysrc.ms. First, presumably one does not want to work on the actual sky model, but one wants to conduct the test on a derivative file. Outside of CASA, cp -r mysrc.model mysrc.test creates a new sky model that is initially identical to the true sky model. In this example, only a single "fake" source will be added, but the extension to multiple sources is straightforward.

For the purposes that motivated this test, I wanted to add weak sources to the visibility data to determine if I could detect them reliably. Thus, I used immath to maintain the image header, but zero out all of the values in the image.

mysrc = 'NGC4536'

imgname = mysrc + '.image'
newimgname = mysrc + '.add'

# create a new blank image
immath(imagename=imgname,outfile=newimgname,
       mode='evalexpr',expr='IM0*0.0')

After this step, the image NGC4536.add is identically zero, but has the same header (coordinate system, pointing direction, etc.) as the original file.

Now, a component list is created and associated with that file. The component list is created using the addcomponent part of the toolkit; eventually, it may be possible to use the asciitocomponentlist part of the toolkit to do this for multiple sources simultaneously.

cl.open()
cl.addcomponent(flux=1.25, fluxunit='mJy', polarization='Stokes',
                dir='J2000 19h30m00 15d00m00', shape='gaussian', majoraxis='10arcsec',
                minoraxis='6arcsec', positionangle='0deg', freq='1.25GHz',
                spectrumtype='spectral index', index=-0.8)
### you can add more components if you wish by calling addcomponent
### repeatedly with different params

See [[1]] for the full list of options to addcomponent.

If desired, one can save the component list to disk, which would allow it to be viewed or used later.

##save it to disk
cl.rename('my_1_component.cl')

Now the components are added to the image, and everything is closed up.

ia.open(newimgname)
ia.modify(cl.torecord())

cl.close()
cl.done()
ia.close()

The next step is to add these weak sources into the MODEL column of the Measurement Set. There are two ways to do this using ft, one using the image, the other the component list:

ft(vis='myms', complist='my_1_component.cl')

or

ft(vis='myms', model=newimgname)

The final step is the addition of these weak sources to the visibility data, by uvsub:

uvsub(vis='myms',reverse=True)

At this point, the visibility data has a weak source(s) added to it. One can now test whether it can be recovered reliably.

YMMV.
TJWL