ACA Simulation (CASA 4.0): Difference between revisions

From CASA Guides
Jump to navigationJump to search
Jbraatz (talk | contribs)
Jbraatz (talk | contribs)
 
(11 intermediate revisions by the same user not shown)
Line 9: Line 9:


=====Set simobserve as current task=====
=====Set simobserve as current task=====
Reset all parameters to default, and then set the project name to m51c
Reset all parameters to default, and then set the project name to ''m51c''
<source lang="python">
<source lang="python">
# Set simobserve to default parameters
# Set simobserve to default parameters
Line 15: Line 15:
# Our project name will be m51c, and all simulation products will be placed in a subdirectory m51c/
# Our project name will be m51c, and all simulation products will be placed in a subdirectory m51c/
project="m51c"
project="m51c"
</source>  
</source>
 
=====Specify sky model image=====
=====Specify sky model image=====
We'll use an Halpha image of M51 as an example model sky. Download [[File:M51ha.fits.txt]] and place in your working directory, or use the curl command in the script.
We'll use an Halpha image of M51 as a model of the sky, for this example.   The ''curl'' command will copy the file from a URL and rename it.


simobserve will make a copy m51c/m51c.skymodel, and not modify your input image.
<source lang="python">
<source lang="python">
# Model sky = Halpha image of M51  
# Model sky = Halpha image of M51  
os.system('curl http://casaguides.nrao.edu/images/3/3f/M51ha.fits.txt -f -o M51ha.fits.txt')
os.system('curl http://casaguides.nrao.edu/images/3/3f/M51ha.fits.txt -f -o M51ha.fits')
skymodel        =  "M51ha.fits.txt"
skymodel        =  "M51ha.fits"
</source>   
</source>   
Although the image has a world coordinate system, we want to override most of the parameters.
 
* We'll place the source in the southern hemisphere with the indirection parameter,  
Note that '''simobserve''' will not modify your original input image.  Rather, it will make a copy ''m51c/m51c.skymodel''.
* set the pixel size to 0.1arcsec, effectively moving the galaxy further away (M51 itself would require a quite large mosaic, and in any case we need for the input model pixels to be significantly smaller than the synthesized beam that we'll be simulating, or else we won't be learning anything)
 
* consistent with simulating a more distant source, we'll set the peak brightness to 0.004 Jy/pixel
We will override most of the parameters in the Halpha FITS image to make the image more suitable to a sub-millimeter ALMA observation. We will:
* set the frequency to 330GHz, and since its a 2D image we'll set the single "channel" width to be 50MHz, and peak brightness of 0.004 Jy/pixel - parameters plausible for observing an emission line in a galaxy.
* place the source in the southern hemisphere with the ''indirection'' parameter,  
* set the pixel size to 0.1arcsec, to simulate an observation of a galaxy that is smaller in angular size. (M51 itself would require a quite large mosaic, and in any case we'd like the input model pixels to be significantly smaller than the synthesized beam.)
* set the peak brightness to 0.004 Jy/pixel
* set the frequency to 330GHz, and since it's a 2D image we'll set the single "channel" width to be 50MHz, and peak brightness of 0.004 Jy/pixel.  These parameters are plausible for observing a sub-mm emission line in a galaxy.
<source lang="python">
<source lang="python">
# Set model image parameters:
# Set model image parameters:
Line 44: Line 47:
[[Image:M51c.ALMA_0.5arcsec.skymodel.png|thumb|hexagonal mosaic overplotted on sky model]]
[[Image:M51c.ALMA_0.5arcsec.skymodel.png|thumb|hexagonal mosaic overplotted on sky model]]


We'll begin with the 12m ALMA array observation, and have simobserve calculate a hexagonal mosaic of pointings.
We'll begin by simulating the observation as seen by the main 12 m ALMA array.  We'll have '''simobserve''' calculate a hexagonal mosaic of pointings.


The default interface for simobserve provides an integration parameter, which is the dwell time at each mosaic pointing -- we'll set that to 10s.  A real observation would integrate a scan of ~5 min at each mosaic pointing;  we could set integration="5min" but then for data volume and speed, simobserve would only generate one measurement per 5min scanWhile thermal noise levels would be scaled correctly, corruption of the data with a phase screen, or the details of uv coverage, would not be realistic.  Rotating through the mosaic more rapidly than a real simulation will result in more representative uv coverageIf you wish to simulate the more realistic case such as 5 min scans with 5s integrations, please see [[Complex pointingtable (CASA 3.4)]] for a guide to doing that with simobserve (Its not hard, it just takes two calls to the task instead of one).
The interface for '''simobserve''' provides an ''integration'' parameter, which is the dwell time at each pointing in the mosaic.  We'll set that to 10 s.  A real observation might dwell on each pointing for ~ 20 s and record 3 sets of visibilities (i.e. 6 s integration time) for each pointing'''simobserve''' does not currently offer a capability to reproduce this behavior exactly, but the difference in UV coverage should be very minor between the real observation and the simulated oneYou could try setting ''integration'' to ~ 20 s, which would more closely match the dwell time for some real observations.


<source lang="python">
<source lang="python">
Line 58: Line 61:


=====Calculate Visibilities=====
=====Calculate Visibilities=====
simobserve can determine what array configuration to use, if you provide a desired resolution or synthesized beam size.   
When you provide a desired resolution in the ''antennalist'' parameter, '''simobserve''' will determine what array configuration to use.   
<source lang="python">
<source lang="python">
obsmode            =  "int"
obsmode            =  "int"
Line 65: Line 68:
</source>
</source>


run simobserve, displaying graphics to screen and file (files can be found in the project directory, e.g. m51c)
Run '''simobserve''', displaying graphics to the screen and also writing graphics output to files that can be found in the project directory, ''m51c''.


<source lang="python">
<source lang="python">
Line 77: Line 80:
[[Image:M51c.aca.tp.skymodel.png|thumb|rectangular map overplotted on sky model]]
[[Image:M51c.aca.tp.skymodel.png|thumb|rectangular map overplotted on sky model]]


Next we'll simulate a total power raster map of the same area, but on a more realistic square grid.  CASA simulation tools can not simulate  
Next we'll simulate a total power raster map of the same area, on a square grid.  CASA simulation tools can not simulate  
true on-the-fly mapping (with smearing on timescales smaller than an integration time), but a square grid with a short integration time  
true on-the-fly mapping (with smearing on timescales smaller than an integration time), but a square grid with a short integration time will provide a very accurate approximation.
will provide a very accurate approximation.


By virtue of CASA's global parameters, we already have project and image world coordinate system parameters set correctly.   
By virtue of CASA's global parameters, we already have project and image world coordinate system parameters set correctly.   


We need to change the pointing calculation to make it square and a bit larger than the interferometric map.
We need to change the mapping parameters to specify a square region a bit larger than the interferometric map.
<source lang="python">
<source lang="python">
integration        =  "10s"
integration        =  "10s"
Line 90: Line 92:
</source>  
</source>  


We'll observe on a different day (this doesn't really matter, but if you choose to simulate two different 12m ALMA configurations and combine them, if they're simulated on the same day with the same antenna names you will have issues concatenating the datasets, so its a good habit to change the day)
We'll observe on a different day.  This doesn't really matter, but if you choose to simulate two different 12m ALMA configurations and combine them, you will have issues concatenating the datasets if they are simulated on the same day with the same antenna names. So, it's a good habit to change the day.


<source lang="python">
<source lang="python">
Line 100: Line 102:
</source>
</source>


run simobserve, displaying graphics to screen and file
Run '''simobserve''', displaying graphics to screen and to files.


<source lang="python">
<source lang="python">
Line 120: Line 122:
</source>  
</source>  


We can specify an integral number of times to repeat the mosaic by setting totaltime to an integer string without units.
We can specify an integral number of times to repeat the mosaic by setting ''totaltime'' to an integer string without units.


<source lang="python">
<source lang="python">
Line 129: Line 131:
</source>
</source>


run simobserve, displaying graphics to screen and file
Run '''simobserve''', displaying graphics to screen and to files.


<source lang="python">
<source lang="python">
Line 139: Line 141:
====Deconvolve the visibilities back into an image====
====Deconvolve the visibilities back into an image====


We use simanalyze to take the three measurement sets and create a single image.  
Next we use '''simanalyze''' to combine the three measurement sets and create a single image.  


There are many ways to do this, and the project is NOT making any recommendation yet at this time about what is optimalPlease discuss with your ARC contact scientist if you have ALMA data now, or wait for additional recommendations to be posted here over time.
There are many ways to do this, and you may wish to discuss options with scientists at your ARC. In this example we will use the total power image as a model when deconvolving the ACA image, and then use the result as a model when deconvolving the 12m interferometric imageThis method tends to give low weight to the large spatial scales, but is simple to illustrate.


* We will use the total power image as a model when deconvolving the ACA image, and then use the result as a model when deconvolving the 12m interferometric image.  This tends to give low weight to the large spatial scales, but is simple to illustrate.
It's possible to get better results if one used multiscale clean in the clean task (again using the lower resolution image as a model when deconvolving the higher resolution one)An alternative would be to create an image independently from each dataset, and then use the CASA feather task to combine them entirely in the image plane.
* Improved results would result if one used multiscale clean, in the clean task (again using the lower resolution image as a model when deconvolving the higher resolution one)
* An alternative would be to create each image independently, and then use the CASA feather task to combine them entirely in the image plane.


simanalyze, if given a total power and interferometric measurement set, will automatically create the total power image,  
If given a total power and interferometric measurement set, '''simanalyze''' will automatically create the total power image,  
then use it as a model and deconvolve the interferometric image.  It is not recommended to do both interferometric images simultaneously.
then use it as a model and deconvolve the interferometric image.  It is not recommended to do both interferometric images simultaneously.


Line 168: Line 168:
====Next add the 12m interferometric data====
====Next add the 12m interferometric data====
[[Image:M51c.ALMA_0.5arcsec.analysis.3.4.png|thumb|TP+ACA combined (with relatively low weight) with 12m ALMA]]
[[Image:M51c.ALMA_0.5arcsec.analysis.3.4.png|thumb|TP+ACA combined (with relatively low weight) with 12m ALMA]]
Here we explicitly have to set the modelimage to the previous output.
Here we explicitly have to set the ''modelimage'' to the result from the previous run.


<source lang="python">
<source lang="python">

Latest revision as of 22:45, 28 November 2012

Simulating Observations in CASA

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

ALMA 12m + ACA + Total Power

We can simulate ALMA observations that use the main array of 12 m antennas augmented by the ACA and Total Power antennas by generating simulated observations for each component separately, and then combining the Measurement Sets in simanalyze. This technique is general and can be used to simulate observations using multiple 12 m array configurations, as well. Total power observations can be simulated either in an independent run of simobserve, or along with an interferometric simulation. Note that if you simulate total power and an interferometric observation simultaneously with simobserve, they must have the same set of pointing centers and the same integration and total time, which is probably not realistic. (For example it is generally recommended to observe a larger area by 1/2 primary beam in total power mode to combine with a 12 m ALMA mosaic).

Set simobserve as current task

Reset all parameters to default, and then set the project name to m51c

# Set simobserve to default parameters
default("simobserve")
# Our project name will be m51c, and all simulation products will be placed in a subdirectory m51c/
project="m51c"
Specify sky model image

We'll use an Halpha image of M51 as a model of the sky, for this example. The curl command will copy the file from a URL and rename it.

# Model sky = Halpha image of M51 
os.system('curl http://casaguides.nrao.edu/images/3/3f/M51ha.fits.txt -f -o M51ha.fits')
skymodel         =  "M51ha.fits"

Note that simobserve will not modify your original input image. Rather, it will make a copy m51c/m51c.skymodel.

We will override most of the parameters in the Halpha FITS image to make the image more suitable to a sub-millimeter ALMA observation. We will:

  • place the source in the southern hemisphere with the indirection parameter,
  • set the pixel size to 0.1arcsec, to simulate an observation of a galaxy that is smaller in angular size. (M51 itself would require a quite large mosaic, and in any case we'd like the input model pixels to be significantly smaller than the synthesized beam.)
  • set the peak brightness to 0.004 Jy/pixel
  • set the frequency to 330GHz, and since it's a 2D image we'll set the single "channel" width to be 50MHz, and peak brightness of 0.004 Jy/pixel. These parameters are plausible for observing a sub-mm emission line in a galaxy.
# Set model image parameters:
indirection="B1950 23h59m59.96s -34d59m59.50s"
incell="0.1arcsec"
inbright="0.004"
incenter="330.076GHz"
inwidth="50MHz"

Simulate 12m interferometric observation

hexagonal mosaic overplotted on sky model

We'll begin by simulating the observation as seen by the main 12 m ALMA array. We'll have simobserve calculate a hexagonal mosaic of pointings.

The interface for simobserve provides an integration parameter, which is the dwell time at each pointing in the mosaic. We'll set that to 10 s. A real observation might dwell on each pointing for ~ 20 s and record 3 sets of visibilities (i.e. 6 s integration time) for each pointing. simobserve does not currently offer a capability to reproduce this behavior exactly, but the difference in UV coverage should be very minor between the real observation and the simulated one. You could try setting integration to ~ 20 s, which would more closely match the dwell time for some real observations.

# have simobserve calculate mosaic pointing locations:
setpointings       =  True
integration        =  "10s"
mapsize            =  "1arcmin"
maptype            =  "hex"
pointingspacing    =  "9arcsec"      # this could also be specified in units of the primary beam e.g. "0.5PB"
Calculate Visibilities

When you provide a desired resolution in the antennalist parameter, simobserve will determine what array configuration to use.

obsmode            =  "int"
antennalist        =  "ALMA;0.5arcsec"
totaltime          =  "3600s"

Run simobserve, displaying graphics to the screen and also writing graphics output to files that can be found in the project directory, m51c.

graphics           =  "both"
simobserve()

Simulate 12m total power observation

rectangular map overplotted on sky model

Next we'll simulate a total power raster map of the same area, on a square grid. CASA simulation tools can not simulate true on-the-fly mapping (with smearing on timescales smaller than an integration time), but a square grid with a short integration time will provide a very accurate approximation.

By virtue of CASA's global parameters, we already have project and image world coordinate system parameters set correctly.

We need to change the mapping parameters to specify a square region a bit larger than the interferometric map.

integration        =  "10s"
mapsize            =  "1.3arcmin"
maptype            =  "square"

We'll observe on a different day. This doesn't really matter, but if you choose to simulate two different 12m ALMA configurations and combine them, you will have issues concatenating the datasets if they are simulated on the same day with the same antenna names. So, it's a good habit to change the day.

obsmode            = "sd"
sdantlist          = "aca.tp.cfg"
sdant              = 0
refdate            = "2012/12/01"
totaltime          =  "2h"

Run simobserve, displaying graphics to screen and to files.

simobserve()

Simulate 7m ACA observation

hexagonal map overplotted on sky model

Next we'll add an ACA mosaic, with its larger primary beam.

integration        =  "10s"
mapsize            =  "1arcmin"
maptype            =  "hex"
pointingspacing    =  "15arcsec"

We can specify an integral number of times to repeat the mosaic by setting totaltime to an integer string without units.

obsmode            = "int"
refdate            = "2012/12/02"
antennalist        =  "aca.i.cfg"
totaltime          =  "3"

Run simobserve, displaying graphics to screen and to files.

simobserve()

Deconvolve the visibilities back into an image

Next we use simanalyze to combine the three measurement sets and create a single image.

There are many ways to do this, and you may wish to discuss options with scientists at your ARC. In this example we will use the total power image as a model when deconvolving the ACA image, and then use the result as a model when deconvolving the 12m interferometric image. This method tends to give low weight to the large spatial scales, but is simple to illustrate.

It's possible to get better results if one used multiscale clean in the clean task (again using the lower resolution image as a model when deconvolving the higher resolution one). An alternative would be to create an image independently from each dataset, and then use the CASA feather task to combine them entirely in the image plane.

If given a total power and interferometric measurement set, simanalyze will automatically create the total power image, then use it as a model and deconvolve the interferometric image. It is not recommended to do both interferometric images simultaneously.

First image total power and ACA with total power as a model

Total power combined (with relatively low weight) with ACA
default("simanalyze")
project            =  "m51c"
vis                =  '$project.aca.i.ms,$project.aca.tp.sd.ms'  
imsize             =  [512,512]
cell               =  '0.2arcsec'
analyze            =  True
showpsf            =  False
showresidual       =  False
showconvolved      =  True
simanalyze()

Next add the 12m interferometric data

TP+ACA combined (with relatively low weight) with 12m ALMA

Here we explicitly have to set the modelimage to the result from the previous run.

default("simanalyze")
project            =  "m51c"
vis                =  '$project.ALMA_0.5arcsec.ms'
imsize             =  [512,512]
cell               =  '0.1arcsec'
modelimage         =  "$project.aca.i.image"
analyze            =  True
showpsf            =  False
showresidual       =  False
showconvolved      =  True
simanalyze()

Simulating Observations in CASA

Last checked on CASA Version 4.0.0.