Simulation Guide Component Lists (CASA 3.3): Difference between revisions

From CASA Guides
Jump to navigationJump to search
Sschnee (talk | contribs)
Created page with "{{Simulations Intro}} ''This guide is applicable to CASA version 3.3. For simulation instructions for older versions of CASA see Simdata New Users Guide (CASA 3.2).'' To c..."
 
Sbooth (talk | contribs)
 
(59 intermediate revisions by 3 users not shown)
Line 1: Line 1:
{{Simulations Intro}}
{{Simulations Intro}}


''This guide is applicable to CASA version 3.3.  For simulation instructions for older versions of CASA see [[Simdata New Users Guide (CASA 3.2)]].''
''This guide is applicable to CASA version 3.3.''


To create a script of the Python code on this page see [[Extracting scripts from these tutorials]].
To create a script of the Python code on this page see [[Extracting scripts from these tutorials]].
Line 7: Line 7:
==Explanation of the guide==
==Explanation of the guide==


When planning an interferometric observation it is useful to simulate the output of the instrument under a variety of conditions.  In CASA, interferometric observations can be simulated using task <tt>sim_observe</tt> and quickly analyzed using task <tt>sim_analyze</tt>.  This guide will demonstrate how to simulate an ALMA observation using these tasks. (Task <tt>sim_observe</tt> can also be used to simulate observations with other interferometers, but this currently requires advanced techniques not covered here.)
When writing an interferometric proposal it is often useful to simulate observations of very simple objects, like point sources, Gaussians, and disks.  In CASA, observations can be simulated using task <tt>sim_observe</tt> and analyzed using task <tt>sim_analyze</tt>.  This guide will demonstrate how to simulate ALMA observations of a Gaussian and some point sources using these tasks as well as the [http://casa.nrao.edu/docs/CasaRef/CasaRef.html CASA Toolkit].


We begin with an image similar to something that might be observed with ALMAWe show how to rescale the image and specify the ALMA antenna configuration. We produce a simulated set of visibility measurements and then produce an image from the visibility data.  In the process, we also generate useful figures to help us analyze the simulated output. 
We begin by using component lists in the Toolkit to create an image of a Gaussian flux distribution, and this will be saved as a FITS fileThe fits file will them be "observed" using <tt>sim_observe</tt> and <tt>sim_analyze</tt> along with four point sources, added via the <tt>componentlist</tt> parameter.   Finally, we show how the same observations could have been done without any skymodel in <tt>sim_observe</tt>, instead using only component lists.
 
We will assume throughout only a general knowledge of interferometry and no specific knowledge of CASA.


==Getting Started==
==Getting Started==


The two things you need to get started are
To get started you need CASA version 3.3  
# The image we will work with, and
# CASA version 3.3.
 
We will use the Spitzer IRAC 8 micron image of 30 Doradus, which you can download from the [http://casaguides.nrao.edu/index.php?title=Sim_Inputs Simulation Inputs CASA Guide page].


To install CASA, follow the instructions given on the [http://casa.nrao.edu/casa_obtaining.shtml Obtaining CASA page].
To install CASA, follow the instructions given on the [http://casa.nrao.edu/casa_obtaining.shtml Obtaining CASA page].
Line 25: Line 19:
==CASA Basics==
==CASA Basics==


CASA is the post-processing package for ALMA and EVLA and can handle both interferometric and single dish data.  Because <tt>sim_observe</tt> and <tt>sim_analyze</tt> are tasks within CASA, we start here with a brief introduction to some CASA basics. To learn much more about CASA, go to the [http://casa.nrao.edu CASA homepage].  Walk-throughs of CASA data reduction for a variety of data sets can be found on the [http://casaguides.nrao.edu CASA Guides website].
CASA is the post-processing package for ALMA and EVLA and can handle both interferometric and single dish data.  To get a brief introduction to <tt>sim_observe</tt> and <tt>sim_analyze</tt>, the tasks within CASA that we will use here, go to the [http://casaguides.nrao.edu/index.php?title=Simulation_Guide_for_New_Users_(CASA_3.3) Simulation Guide for New Users]. To learn more about CASA in general, go to the [http://casa.nrao.edu CASA homepage].  Walk-throughs of CASA data reduction for a variety of data sets can be found on the [http://casaguides.nrao.edu CASA Guides website].


Once you have installed CASA, you can launch it by typing "casapy" at the prompt or by double-clicking on the icon, depending on your system and preferences.
Once you have installed CASA, you can launch it by typing "casapy" at the prompt or by double-clicking on the icon, depending on your system and preferences.


[[Image:Simdata_inp_example.png|thumb| Screen shot of <tt>inp sim_observe</tt>.  Note that default values are shown in black, valid and non-default values are shown in blue, and invalid values are shown in red. Expandable parameters have a grey background.]]
==Making a Simple FITS Image==


To see a list of all available CASA tasks, at the CASA prompt type
Here we show how to create a simple FITS image using the CASA tasks and the toolkit.  The example here will be that of a Gaussian flux distribution.  Enter the following lines at the CASA prompt:
tasklist


To look at the inputs for an available task, use "inp". For instance
<source lang="python">
inp sim_observe
# In CASA
Any parameter with a grey background is expandable. Black text shows a default value. Red text shows an invalid value. Blue text shows a valid value different from the default.
direction = "J2000 10h00m00.0s -30d00m00.0s"
cl.done()
cl.addcomponent(dir=direction, flux=1.0, fluxunit='Jy', freq='230.0GHz', shape="Gaussian",
                majoraxis="0.1arcmin", minoraxis='0.05arcmin', positionangle='45.0deg')
ia.fromshape("Gaussian.im",[256,256,1,1],overwrite=True)
cs=ia.coordsys()
cs.setunits(['rad','rad','','Hz'])
cell_rad=qa.convert(qa.quantity("0.1arcsec"),"rad")['value']
cs.setincrement([-cell_rad,cell_rad],'direction')
cs.setreferencevalue([qa.convert("10h",'rad')['value'],qa.convert("-30deg",'rad')['value']],type="direction")
cs.setreferencevalue("230GHz",'spectral')
cs.setincrement('1GHz','spectral')
ia.setcoordsys(cs.torecord())
ia.setbrightnessunit("Jy/pixel")
ia.modify(cl.torecord(),subtract=False)
exportfits(imagename='Gaussian.im',fitsimage='Gaussian.fits',overwrite=True)
</source>


To get help on a given task, use "help" (to exit from help, hit the "q" key). For instance
[[Image:Gauss_fits.jpg|thumb| Gaussian flux distribution created by the CASA Toolkit via component list.]]
help sim_analyze


To reset a task to its default values, use "default". For instance
The first line defines a string "direction" which will be the center of the Gaussian flux distribution.
default sim_observe


To run a task using the current global values of its input parameters, just type its name at the CASA prompt. For instance
<tt>cl.done</tt> closes any open component lists, if any.
sim_observe


It's a good idea to double-check these values (i.e., re-run <tt>inp sim_observe</tt>) immediately before running a task.
<tt>cl.addcomponent</tt> creates a new component centered at "direction", with a flux of 1 Jy at a frequency of 230 GHz, a Gaussian shape of 0.1 by 0.05 arcminutes with a position angle of 45 degrees.


==Using <tt>sim_observe</tt>==
<tt>ia.fromshape</tt> creates a new, empty CASA image with the name and dimensions given.


The task <tt>sim_observe</tt> uses a model image along with a number of input parameters that define the observational conditions and alter the input model.  The output of <tt>sim_observe</tt> is a visibility measurement set and various diagnostic plots and images.
<tt>cs.coordsys</tt> gets the coordinate system of the image.


===Getting Your Input Image Into <tt>sim_observe</tt>===
<tt>cs.setunits</tt> defines the units of the four axes of the new CASA image.


First, we'll tell <tt>sim_observe</tt> where to find the model (input) image and how to scale it appropriately for our purposes.  Just to be safe, we'll begin by restoring the default values of <tt>sim_observe</tt> and then set the 30 Doradus image as the skymodel (Note: you might need to include the path to your model image, if you are not currently in the working directory where the image is).
<tt>cell_rad</tt> will be the cell size and units in this CASA image, 0.1"
 
<source lang="python">
#In CASA
default sim_observe
skymodel = '30dor.fits'
</source>


====Angular Scale====
<tt>cs.setincrement</tt> tells CASA that RA increases to the right, Dec increases going up, and in a few lines that the one channel is 1 GHz wide.  
If you open the FITS image of 30 Doradus in your favorite viewer (e.g. <tt>viewer</tt> in CASA), you will see that it covers quite a large footprint on the sky, about 10' on a side. We are going to tell <tt>sim_observe</tt> to rescale the pixels to shrink the image by roughly a factor of 15 (from  2.3" to 0.15" pixels through the <tt>incell</tt> parameter) so that the model is approximately 40" on a side.  This rescaled model will fit within a small mosaic of 6 pointings.  Although we do this primarily for convenience in this example, a scientific motivation for this type of rescaling would be to approximate what a super-giant HII region like 30 Doradus would look like if moved from the Large Magellanic Cloud to the distance of M33 or M31.  For the sake of demonstration, we will also change the coordinates of the center of the map (using the <tt>indirection</tt> parameter). Units are case sensitive, so please take care to enter them in verbatim.  CASA will throw an error otherwise.


<source lang="python">
<tt>cs.setreferencevalue</tt> sets the center of the image in RA, Dec, and frequency.
#In CASA
incell = '0.15arcsec'
indirection = 'J2000 10h00m00 -40d00m00'
</source>


====Observed Wavelength====
<tt>ia.setcoordsys</tt> puts the coordinates and frequencies into the image header.


The model image of 30 Doradus shows 8 micron continuum emission.  ALMA does not observe at wavelengths this short, so we will tell <tt>sim_observe</tt> that this is actually a 230 GHz (1.3 mm) continuum map.  We will also tell <tt>sim_observe</tt> that the observations were taken with a 2-GHz bandwidth. Although for this particular example the channel width is not a critical number, it would be very important if we were modifying a spectral cube instead of a continuum image.
<tt>ia.setbrightnessunit</tt> defines the brightness unit (Jy per pixel) of the CASA image.


<source lang="python">
<tt>ia.modify</tt> puts the Gaussian component into the image
#In CASA
incenter = '230GHz'
inwidth = '2GHz'
</source>


====Brightness Scale====
<tt>exportfits</tt> writes the resultant CASA image as a FITS file (not strictly necessary for this guide, but useful to know in general).


The 8 micron emission is probably not a great approximation for the millimeter emission from 30 Doradus.  For a true science case, one would want to calculate what the expected 230 GHz emission would be from an object like this at the distance of about 750 kiloparsec.  For the sake of simplicity, we will rescale the image so that the brightest pixel in the map has a flux density of 0.06 mJy.  This number is chosen such that the extended emission is a factor of a few brighter than the expected noise in a 2 hour observation.  The [https://almascience.nrao.edu/call-for-proposals/sensitivity-calculator ALMA sensitivity calculator] can be used to determine the expected noise for an observation.
As usual, more information can be found via the help in CASA


<source lang="python">
<source lang="python">
#In CASA
# In CASA
inbright = '0.06mJy/pixel'
help(ia.modify)  # syntax for help with toolkit or CASA tasks
help("exportfits") # syntax for help with CASA tasks, but not the toolkit
</source>
</source>


Remember that you can check your task inputs at any time by typing <tt>inp</tt> at the casa prompt.
==Simulating Observations with a FITS Image and a Component List==
 
===Defining the Mock Observations===
 
Now that <tt>sim_observe</tt> knows how to interpret the input image, the next step is to define the simulated observations. 


====Pointings and Scan Time====
One use for component lists would be to simulate the effect of having one or more point sources added to an input image, with the goal of finding out the effect on the simulated observations.  For instance, one might want to know if a faint point source would be detectable if there is extended emission around it.  Conversely, one might want to know if the artifacts from imaging a field with a bright point source would make a project tricky to carry out.  In this example, we will use component lists in the <tt>sim_observe</tt> task add four point sources to an input FITS image.  The input image will be the Gaussian flux distribution created above.


We will first change the parameters within <tt>setpointings</tt>. 
First we create the four point sources using the CASA toolkit.
 
The default value for <tt>integration</tt>, 10 seconds, might be appropriate to simulate real observations.  However, <tt>sim_observe</tt> will run much faster if the value for "integration" is increased, reducing the number of data points to be generated.  In this demonstration we will set it to 600 seconds.  When <tt>sim_observe</tt> is used for scientific purposes it may be best to set <tt>integration</tt> to a large value at first to make sure that <tt>sim_observe</tt> runs as expected, and then decrease <tt>integration</tt> to a more realistic time for the final run.  You will get a more accurate simulation with 10 second integrations than with 600 second integrations, especially in Early Science observations with a limited number of baselines.


<source lang="python">
<source lang="python">
#In CASA
# In CASA
integration = '600s'
os.system('rm -rf point.cl')
cl.done()
cl.addcomponent(dir="J2000 10h00m00.08s -30d00m02.0s", flux=0.1, fluxunit='Jy', freq='230.0GHz', shape="point")
cl.addcomponent(dir="J2000 09h59m59.92s -29d59m58.0s", flux=0.1, fluxunit='Jy', freq='230.0GHz', shape="point")
cl.addcomponent(dir="J2000 10h00m00.40s -29d59m55.0s", flux=0.1, fluxunit='Jy', freq='230.0GHz', shape="point")
cl.addcomponent(dir="J2000 09h59m59.60s -30d00m05.0s", flux=0.1, fluxunit='Jy', freq='230.0GHz', shape="point")
cl.rename('point.cl')
cl.close()
</source>
</source>


Note that the integration time is different than the total (on-source) time of the observations. The integration time, set here, is the averaging time for each data pointThe total time spent on-source is set below.  Each pair of antennas will generate a number of data points equal to the total observing time divided by the integration time.
First we delete any previous version of the file 'point.cl', which will be the output file created in a few linesThen we use <tt>cl.done</tt> to begin and end this sequence to close any open component list.  The <tt>cl.addcomponent</tt> commands create point sources that are 0.1 Jy at 230 GHz at the coordinates given in each lineThe <tt>cl.rename</tt> command tells CASA the name of output component list file.
 
We will keep <tt>direction</tt> at the default (blank) value to center the observations on the model coordinates, as given in <tt>indirection</tt> above.
 
We will keep <tt>mapsize</tt> at the default value so that the mosaic will automatically cover the entire imageIn our case, this will require a mosaic of 6 pointings (as we will see later on).  One could also set an exact output image size via <tt>mapsize = ['10arcmin','10arcmin']</tt> for example.
 
The default mosaic pattern, <tt>maptype = 'ALMA'</tt>, tells <tt>sim_observe</tt> to use the same hexagonal algorithm as the ALMA OT. We will leave this unchanged.
 
Finally, we will also leave <tt>pointingspacing</tt> to its default value (blank), which automatically sets the pointings to be half a primary beam apart, corresponding to Nyquist sampling.
 
====Antenna Positions and Total Observation Time====


We will now consider the subparameters available when <tt>observe = True</tt>.  We will keep the default values for every parameter (including <tt>totaltime = 7200s</tt>) except <tt>antennalist</tt>, which tells <tt>sim_observe</tt> the locations and sizes of each antenna in the array.  We will simulate an Early Science observation, so we will first find the location where CASA has stored the ALMA configuration files.  Then we will tell <tt>sim_observe</tt> to use the configuration file designed for Early Science.
We use <tt>sim_observe</tt> and <tt>sim_analyze</tt> to make the simulated observations of these point sources and the Gaussian flux distribution given in the FITS file we made previously in this guide.


<source lang="python">
<source lang="python">
#In CASA
# In CASA
default("sim_observe")
project = "FITS_list"
skymodel = "Gaussian.fits"
inwidth = "1GHz"
complist = 'point.cl'
compwidth = '1GHz'
direction = "J2000 10h00m00.0s -30d00m00.0s"
repodir = os.getenv("CASAPATH").split(' ')[0]
repodir = os.getenv("CASAPATH").split(' ')[0]
antennalist = repodir+"/data/alma/simmos/alma.cycle0.compact.cfg"
antennalist = repodir+'/data/alma/simmos/alma.cycle0.compact.cfg'
totaltime = "28800s"
mapsize = ["20arcsec","20arcsec"]
sim_observe()
 
default("sim_analyze")
project = "FITS_list"
vis="FITS_list.alma.cycle0.compact.ms"
imsize = [256,256]
imdirection = "J2000 10h00m00.0s -30d00m00.0s"
cell = '0.1arcsec'
niter = 5000
threshold = '10.0mJy/beam'
analyze = True
sim_analyze()
</source>
</source>


The first line, run within CASA, uses a Linux command to determine the path for CASAThe second line sets the <tt>antennalist</tt> parameter to the Early Science configuration.  Other <tt>.cfg</tt> files in the same directory exist for various ALMA Full Science array configurations as well as the configuration files for other radio interferometers.  
To learn more about <tt>sim_observe</tt> and <tt>sim_analyze</tt>, look at the [http://casaguides.nrao.edu/index.php?title=Simulation_Guide_for_New_Users_(CASA_3.3) Simulation Guide for New Users]Here we simulate observations of the Gaussian flux distribution (given with the <tt>skymodel</tt> parameter) with the point sources (given with the <tt>complist</tt> parameter) using <tt>sim_observe</tt>.  We center the observations at the center of the Gaussian flux distribution and the point sources, and the field to be imaged is 20" by 20".  We simulate 8 hours (28800 seconds) of observations with ALMA in the Cycle 0 compact array configuration.  The <tt>sim_observe</tt> task puts the u-v data in a measurement set called "FITS_list.alma.cycle0.compact.ms", and this will be the input to make a map of the emission.


Most of the other parameters are not relevant for this simulation as we are not using a component list to describe the sky emission, we do not want to simulate observations of a calibrator, and we do not want to simulate single dish observationsNote that at this stage of simulation development, the time portion of the <tt>refdate</tt> parameter is ignored and all observations are instead centered around transit on the date specified.
Inverting the u-v data and making a clean (deconvolved) map of the flux distribution is done with the <tt>sim_analyze</tt> task in CASA.  The purpose of this guide is to illustrate the use of component lists, not to make the best possible image of this simple flux distribution, so we don't take much care to do the best possible job with the cleaning.  For instance, we don't use a mask or clean boxes, nor do we clean interactively to make sure that there is no point in cleaning furtherTo see a thorough explanation of imaging and deconvolution, see the CASA guides for ALMA science verification data, for instance [http://casaguides.nrao.edu/index.php?title=Antennae_Band7_-_Imaging Antennae Band 7] or [http://casaguides.nrao.edu/index.php?title=TWHydraBand7_Imaging TW Hydra Band 7].


====Thermal Noise====
<gallery widths="280px" heights="180">
File:input_Gauss_point.jpg|''skymodel Gaussian flux distribution with four point sources added as a component list, created using sim_observe.
File:analyze_fits_list.jpg|''simulated observations of a Gaussian (skymodel) and four point sources (component list).
</gallery>


For this simple simulation, we will not include any thermal noise in the observations, so we can leave <tt>thermalnoise</tt> at its default (blank) value.
==Simulating Observations with Just a Component List==


===<tt>sim_observe</tt> Execution and Output===
The CASA task <tt>sim_observe</tt> can be run without an image given entirely by the <tt>complist</tt> parameter with nothing given in the <tt>skymodel</tt> parameter.  A simulation of observations of a Gaussian plus four point sources can be accomplished with only component lists, as shown in the example below.
 
With all the input parameters set, we are ready to execute the task:


<source lang="python">
<source lang="python">
#In CASA
# In CASA
sim_observe()
os.system('rm -rf Gauss_point.cl')
cl.done()
cl.addcomponent(dir="J2000 10h00m00.00s -30d00m00.0s", flux=1.0, fluxunit='Jy', freq='230.0GHz', shape="Gaussian",
                majoraxis="0.1arcmin", minoraxis='0.05arcmin', positionangle='45.0deg')
cl.addcomponent(dir="J2000 10h00m00.08s -30d00m02.0s", flux=0.1, fluxunit='Jy', freq='230.0GHz', shape="point")
cl.addcomponent(dir="J2000 09h59m59.92s -29d59m58.0s", flux=0.1, fluxunit='Jy', freq='230.0GHz', shape="point")
cl.addcomponent(dir="J2000 10h00m00.40s -29d59m55.0s", flux=0.1, fluxunit='Jy', freq='230.0GHz', shape="point")
cl.addcomponent(dir="J2000 09h59m59.60s -30d00m05.0s", flux=0.1, fluxunit='Jy', freq='230.0GHz', shape="point")
cl.rename('Gauss_point.cl')
cl.close()
</source>
</source>


All <tt>sim_observe</tt> output will be written to a directory whose name was given by the <tt>project</tt> parameter, in our case, <tt>project = 'sim'</tt>Inside this directory, you will find
Here we have created a Gaussian flux distribution with the same properties as in the FITS image created aboveSurrounding the Gaussian are four point sources with the same positions and brightness as before.
[[File:Sim.observe.png|thumb|right|The 2x2 summary plot output by <tt>sim_observe</tt> showing source elevation vs. time, antenna position, uv coverage, and the point spread function.]]
# The simulated measurement set (<tt>sim.ms</tt>),
# a CASA image of the point spread function (<tt>sim.quick.psf</tt>),
# a CASA image of the input sky model rescaled according to the <tt>skymodel</tt> sub-parameters (<tt>sim.skymodel</tt>),
# a flattened CASA image of the input sky model rescaled (<tt>sim.skymodel.flat</tt>),
# a flattened PNG image of the input sky model rescaled and overlaid with the mosaic pattern specified by the <tt>setpointings</tt> sub-parameters (<tt>sim.skymodel.png</tt>),
# a 2x2 PNG summary plot (<tt>sim.observe.png</tt>) showing,
## source elevation vs. time,
## antenna position,
## uv coverage,
## and the point spread function,
# and an ASCII text listing of mosaic pointings.
When <tt>sim_observe</tt> is executed, the 2x2 PNG plot will be displayed in the CASA plotter.
 
==Using <tt>sim_analyze</tt>==
 
Now that <tt>sim_observe</tt> has created a simulated set of visibility measurements, we are ready to image the visibilities and analyze the result.  The task <tt>sim_analyze</tt> has been created to easily perform the imaging and subsequent generation of useful analysis plots and figures.  We begin by resetting the <tt>sim_analyze</tt> inputs to their default values.


<source lang="python">
<source lang="python">
#In CASA
# In CASA
default sim_analyze
default("sim_observe")
</source>
project = "complist_only"
 
complist = 'Gauss_point.cl'
We next setup the imaging and analysis input parameters.
compwidth = '1GHz'
 
direction = "J2000 10h00m00.0s -30d00m00.0s"
===Image===
repodir = os.getenv("CASAPATH").split(' ')[0]
 
antennalist = repodir+'/data/alma/simmos/alma.cycle0.compact.cfg'
We want to make a deconvolved output image, but we don't want to spend too much time optimizing the cleaning. So, all we need to do is make sure the "image" parameter is set to True (which it is by default) and leave all of its sub-parameters at their default values.  Other [http://casaguides.nrao.edu/index.php?title=Main_Page data reduction guides] describe the process of imaging in greater detail.
totaltime = "28800s"
 
mapsize = ["20arcsec","20arcsec"]
For simulations intended for a proposal or scientific analysis, one would almost certainly want to choose a more appropriate cleaning threshold and define the region to be cleaned.  Instructions for how to define the region to be cleaned with the "mask" parameter can be found by typing
sim_observe()
> help clean
or by looking at the [http://casaguides.nrao.edu/index.php?title=Clean Clean CASA Guide page].


===Analyze===
default("sim_analyze")
 
project = "complist_only"
To specify how <tt>sim_analyze</tt> displays the imaged data, set the parameter "analyze" to True and then pick your favorite output formats. If the <tt>graphics</tt> parameter is set to 'screen' or 'both', up to six output formats will be displayed in the CASA plotter. More than six outputs can be written to disk (<tt>graphics = 'file'</tt> or <tt>'both'</tt>), but only six can be displayed in the plotter. In this example we will look at,
vis="complist_only.alma.cycle0.compact.ms"
 
imsize = [256,256]
# the uv coverage in the 2 hour observation,
imdirection = "J2000 10h00m00.0s -30d00m00.0s"  
# the synthesized beam (point spread function),
cell = '0.1arcsec'
# the original sky model (as defined in "modifymodel"),
niter = 5000
# the convolved model (sky model convolved with the synthesized beam),
threshold = '10.0mJy/beam'
# the clean image (the sky as observed with the interferometer after deconvolution),
# and the difference between the clean image and the convolved model.
 
To make these choices, use the following lines in CASA
 
<source lang="python">
#In CASA
analyze = True
analyze = True
showconvolved = True
showfidelity = False
</source>
===<tt>sim_analyze</tt> Execution and Output===
[[File:Sim.analysis 2.png|thumb|right| Plot of the six outputs generated by <tt>sim_analyze</tt>.]]
Execute <tt>sim_analyze</tt> by typing,
<source lang="python">
#In CASA
sim_analyze()
sim_analyze()
</source>
</source>


The six outputs we selected will be displayed in the CASA plotter and written to disk in the <tt>sim</tt> directory. 
[[Image:analyze_list_only.jpg|thumb| Simulated observations of a Gaussian and four point sources, all input via all component lists.]]
 
[[File:30Doradus.sim.image.png|thumb|right| Screen shot of the 30 Doradus image shown in the CASA viewer.]]
 
All CASA images written to disk can be opened later using the CASA viewer.  Just type
> viewer
at the CASA prompt to start the viewer tool, and navigate to the appropriate directory, in our case <tt>sim</tt>.  To display the simulated observations, first click on sim.image, then click the "raster image" button and then click "done".  The image will be shown in the viewer, and clicking on the picture of the wrench in the upper left corner will allow you to alter the image in many ways, such as changing the color scale, changing the coordinate scale, and axis labels.  The image to the right was created by:


# changing <tt>basic settings -> color map</tt> from "Rainbow 2" to "Hot Metal 1",
The only difference between this call to <tt>sim_observe</tt> and the one above is that here the Gaussian flux distribution is included in the component list instead of being defined by a FITS file.  The simulated images are essentially identical, as can be seen here.
# changing <tt>beam ellipse -> beam style</tt> from "outline" to "filled",
# and changing <tt>color wedge -> display color wedge</tt> from "No" to "Yes".  


You can learn much more about the functionality of the CASA viewer by watching this [http://casa.nrao.edu/CasaViewerDemo/casaViewerDemo.html instructional video].


{{Simulations Intro}}
{{Simulations Intro}}


[[Category: Simulations]] [[Category:ALMA]]
[[Category: Simulations]] [[Category:ALMA]]
{{Checked 3.3.0}}

Latest revision as of 17:13, 6 September 2017

Simulating Observations in CASA

This guide is applicable to CASA version 3.3.

To create a script of the Python code on this page see Extracting scripts from these tutorials.

Explanation of the guide

When writing an interferometric proposal it is often useful to simulate observations of very simple objects, like point sources, Gaussians, and disks. In CASA, observations can be simulated using task sim_observe and analyzed using task sim_analyze. This guide will demonstrate how to simulate ALMA observations of a Gaussian and some point sources using these tasks as well as the CASA Toolkit.

We begin by using component lists in the Toolkit to create an image of a Gaussian flux distribution, and this will be saved as a FITS file. The fits file will them be "observed" using sim_observe and sim_analyze along with four point sources, added via the componentlist parameter. Finally, we show how the same observations could have been done without any skymodel in sim_observe, instead using only component lists.

Getting Started

To get started you need CASA version 3.3

To install CASA, follow the instructions given on the Obtaining CASA page.

CASA Basics

CASA is the post-processing package for ALMA and EVLA and can handle both interferometric and single dish data. To get a brief introduction to sim_observe and sim_analyze, the tasks within CASA that we will use here, go to the Simulation Guide for New Users. To learn more about CASA in general, go to the CASA homepage. Walk-throughs of CASA data reduction for a variety of data sets can be found on the CASA Guides website.

Once you have installed CASA, you can launch it by typing "casapy" at the prompt or by double-clicking on the icon, depending on your system and preferences.

Making a Simple FITS Image

Here we show how to create a simple FITS image using the CASA tasks and the toolkit. The example here will be that of a Gaussian flux distribution. Enter the following lines at the CASA prompt:

# In CASA
direction = "J2000 10h00m00.0s -30d00m00.0s"
cl.done()
cl.addcomponent(dir=direction, flux=1.0, fluxunit='Jy', freq='230.0GHz', shape="Gaussian", 
                majoraxis="0.1arcmin", minoraxis='0.05arcmin', positionangle='45.0deg')
ia.fromshape("Gaussian.im",[256,256,1,1],overwrite=True)
cs=ia.coordsys()
cs.setunits(['rad','rad','','Hz'])
cell_rad=qa.convert(qa.quantity("0.1arcsec"),"rad")['value']
cs.setincrement([-cell_rad,cell_rad],'direction')
cs.setreferencevalue([qa.convert("10h",'rad')['value'],qa.convert("-30deg",'rad')['value']],type="direction")
cs.setreferencevalue("230GHz",'spectral')
cs.setincrement('1GHz','spectral')
ia.setcoordsys(cs.torecord())
ia.setbrightnessunit("Jy/pixel")
ia.modify(cl.torecord(),subtract=False)
exportfits(imagename='Gaussian.im',fitsimage='Gaussian.fits',overwrite=True)
Gaussian flux distribution created by the CASA Toolkit via component list.

The first line defines a string "direction" which will be the center of the Gaussian flux distribution.

cl.done closes any open component lists, if any.

cl.addcomponent creates a new component centered at "direction", with a flux of 1 Jy at a frequency of 230 GHz, a Gaussian shape of 0.1 by 0.05 arcminutes with a position angle of 45 degrees.

ia.fromshape creates a new, empty CASA image with the name and dimensions given.

cs.coordsys gets the coordinate system of the image.

cs.setunits defines the units of the four axes of the new CASA image.

cell_rad will be the cell size and units in this CASA image, 0.1"

cs.setincrement tells CASA that RA increases to the right, Dec increases going up, and in a few lines that the one channel is 1 GHz wide.

cs.setreferencevalue sets the center of the image in RA, Dec, and frequency.

ia.setcoordsys puts the coordinates and frequencies into the image header.

ia.setbrightnessunit defines the brightness unit (Jy per pixel) of the CASA image.

ia.modify puts the Gaussian component into the image

exportfits writes the resultant CASA image as a FITS file (not strictly necessary for this guide, but useful to know in general).

As usual, more information can be found via the help in CASA

# In CASA
help(ia.modify)  # syntax for help with toolkit or CASA tasks
help("exportfits") # syntax for help with CASA tasks, but not the toolkit

Simulating Observations with a FITS Image and a Component List

One use for component lists would be to simulate the effect of having one or more point sources added to an input image, with the goal of finding out the effect on the simulated observations. For instance, one might want to know if a faint point source would be detectable if there is extended emission around it. Conversely, one might want to know if the artifacts from imaging a field with a bright point source would make a project tricky to carry out. In this example, we will use component lists in the sim_observe task add four point sources to an input FITS image. The input image will be the Gaussian flux distribution created above.

First we create the four point sources using the CASA toolkit.

# In CASA
os.system('rm -rf point.cl')
cl.done()
cl.addcomponent(dir="J2000 10h00m00.08s -30d00m02.0s", flux=0.1, fluxunit='Jy', freq='230.0GHz', shape="point")
cl.addcomponent(dir="J2000 09h59m59.92s -29d59m58.0s", flux=0.1, fluxunit='Jy', freq='230.0GHz', shape="point")
cl.addcomponent(dir="J2000 10h00m00.40s -29d59m55.0s", flux=0.1, fluxunit='Jy', freq='230.0GHz', shape="point")
cl.addcomponent(dir="J2000 09h59m59.60s -30d00m05.0s", flux=0.1, fluxunit='Jy', freq='230.0GHz', shape="point")
cl.rename('point.cl')
cl.close()

First we delete any previous version of the file 'point.cl', which will be the output file created in a few lines. Then we use cl.done to begin and end this sequence to close any open component list. The cl.addcomponent commands create point sources that are 0.1 Jy at 230 GHz at the coordinates given in each line. The cl.rename command tells CASA the name of output component list file.

We use sim_observe and sim_analyze to make the simulated observations of these point sources and the Gaussian flux distribution given in the FITS file we made previously in this guide.

# In CASA
default("sim_observe")
project = "FITS_list"
skymodel = "Gaussian.fits"
inwidth = "1GHz"
complist = 'point.cl'
compwidth = '1GHz'
direction = "J2000 10h00m00.0s -30d00m00.0s" 
repodir = os.getenv("CASAPATH").split(' ')[0]
antennalist = repodir+'/data/alma/simmos/alma.cycle0.compact.cfg'
totaltime = "28800s"
mapsize = ["20arcsec","20arcsec"]
sim_observe()

default("sim_analyze")
project = "FITS_list"
vis="FITS_list.alma.cycle0.compact.ms"
imsize = [256,256]
imdirection = "J2000 10h00m00.0s -30d00m00.0s" 
cell = '0.1arcsec'
niter = 5000
threshold = '10.0mJy/beam'
analyze = True
sim_analyze()

To learn more about sim_observe and sim_analyze, look at the Simulation Guide for New Users. Here we simulate observations of the Gaussian flux distribution (given with the skymodel parameter) with the point sources (given with the complist parameter) using sim_observe. We center the observations at the center of the Gaussian flux distribution and the point sources, and the field to be imaged is 20" by 20". We simulate 8 hours (28800 seconds) of observations with ALMA in the Cycle 0 compact array configuration. The sim_observe task puts the u-v data in a measurement set called "FITS_list.alma.cycle0.compact.ms", and this will be the input to make a map of the emission.

Inverting the u-v data and making a clean (deconvolved) map of the flux distribution is done with the sim_analyze task in CASA. The purpose of this guide is to illustrate the use of component lists, not to make the best possible image of this simple flux distribution, so we don't take much care to do the best possible job with the cleaning. For instance, we don't use a mask or clean boxes, nor do we clean interactively to make sure that there is no point in cleaning further. To see a thorough explanation of imaging and deconvolution, see the CASA guides for ALMA science verification data, for instance Antennae Band 7 or TW Hydra Band 7.

Simulating Observations with Just a Component List

The CASA task sim_observe can be run without an image given entirely by the complist parameter with nothing given in the skymodel parameter. A simulation of observations of a Gaussian plus four point sources can be accomplished with only component lists, as shown in the example below.

# In CASA
os.system('rm -rf Gauss_point.cl')
cl.done()
cl.addcomponent(dir="J2000 10h00m00.00s -30d00m00.0s", flux=1.0, fluxunit='Jy', freq='230.0GHz', shape="Gaussian",
                majoraxis="0.1arcmin", minoraxis='0.05arcmin', positionangle='45.0deg')
cl.addcomponent(dir="J2000 10h00m00.08s -30d00m02.0s", flux=0.1, fluxunit='Jy', freq='230.0GHz', shape="point")
cl.addcomponent(dir="J2000 09h59m59.92s -29d59m58.0s", flux=0.1, fluxunit='Jy', freq='230.0GHz', shape="point")
cl.addcomponent(dir="J2000 10h00m00.40s -29d59m55.0s", flux=0.1, fluxunit='Jy', freq='230.0GHz', shape="point")
cl.addcomponent(dir="J2000 09h59m59.60s -30d00m05.0s", flux=0.1, fluxunit='Jy', freq='230.0GHz', shape="point")
cl.rename('Gauss_point.cl')
cl.close()

Here we have created a Gaussian flux distribution with the same properties as in the FITS image created above. Surrounding the Gaussian are four point sources with the same positions and brightness as before.

# In CASA
default("sim_observe")
project = "complist_only"
complist = 'Gauss_point.cl'
compwidth = '1GHz'
direction = "J2000 10h00m00.0s -30d00m00.0s" 
repodir = os.getenv("CASAPATH").split(' ')[0]
antennalist = repodir+'/data/alma/simmos/alma.cycle0.compact.cfg'
totaltime = "28800s"
mapsize = ["20arcsec","20arcsec"]
sim_observe()

default("sim_analyze")
project = "complist_only"
vis="complist_only.alma.cycle0.compact.ms"
imsize = [256,256]
imdirection = "J2000 10h00m00.0s -30d00m00.0s" 
cell = '0.1arcsec'
niter = 5000
threshold = '10.0mJy/beam'
analyze = True
sim_analyze()
Simulated observations of a Gaussian and four point sources, all input via all component lists.

The only difference between this call to sim_observe and the one above is that here the Gaussian flux distribution is included in the component list instead of being defined by a FITS file. The simulated images are essentially identical, as can be seen here.


Simulating Observations in CASA

Last checked on CASA Version 3.3.0.