Simalma (CASA 4) original: Difference between revisions

From CASA Guides
Jump to navigationJump to search
Rindebet (talk | contribs)
Created page with "Category: Simulations To create a script of the Python code on this page see Extracting scripts from these tutorials. __NOTOC__ == The simalma task == ALMA consists ..."
 
m Estarr moved page Simalma to Simalma (CASA 4) original: making this a redirect page for newest version
 
(10 intermediate revisions by one other user not shown)
Line 13: Line 13:
New in 4.1 is the '''simalma''' task, which takes one set of parameters describing the region of the sky to observe, and makes the appropriate calls to '''simobserve''' and '''simanalyze''' for the user.  How this works is described here.  
New in 4.1 is the '''simalma''' task, which takes one set of parameters describing the region of the sky to observe, and makes the appropriate calls to '''simobserve''' and '''simanalyze''' for the user.  How this works is described here.  


As of CASA 4.1, ALMA is still optimizing the algorithms for combining total power and interferometric data, so the parameters used here are likely to change.
<font color="purple">As of CASA 4.1, ALMA is still optimizing the algorithms for combining total power and interferometric data, so the parameters used here are likely to change.</font>




UNDER CONSTRUCTION
=====Set simalma as current task=====
 
Reset all parameters to default, and then set the project name to ''m51''
=====Set simobserve as current task=====
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 simalma to default parameters
default("simobserve")
default("simalma")
# Our project name will be m51c, and all simulation products will be placed in a subdirectory m51c/
# Our project name will be m51, and all simulation products will be placed in a subdirectory m51/
project="m51c"
project="m51"
</source>
</source>


Line 36: Line 34:
</source>   
</source>   


Note that '''simobserve''' will not modify your original input image.  Rather, it will make a copy ''m51c/m51c.skymodel''.
Note that '''simalma''' will not modify your original input image.  Rather, it will make a copy ''m51/m51.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:
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:
Line 45: Line 43:
<source lang="python">
<source lang="python">
# Set model image parameters:
# Set model image parameters:
indirection="B1950 23h59m59.96s -34d59m59.50s"
indirection="J2000 23h59m59.96s -34d59m59.50s"
incell="0.1arcsec"
incell="0.1arcsec"
inbright="0.004"
inbright="0.004"
Line 54: Line 52:
----
----


====Simulate 12m interferometric observation====
====Set up Observing Parameters====
[[Image:M51c.ALMA_0.5arcsec.skymodel.png|thumb|hexagonal mosaic overplotted on sky model]]
[[Image:M51.alma_cycle1_3.skymodel.png|thumb|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.
Based on the Cycle 1 capabilities, we would like to use array configuration number 3 which affords ~0.5 arcsec resolution:
<source lang="python">
antennalist="alma_cycle1_3.cfg"
</source>


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.
We'll observe for a couple of hours (this is the 12m array observing time):
<source lang="python">
totaltime="7800s"
</source>


Following the Cycle 1 convention, <tt>simalma</tt> will observe 3 times longer with the 7m array and total power dishes.
<source lang="python">
<source lang="python">
# have simobserve calculate mosaic pointing locations:
acaratio=3.0
setpointings      = True
acaconfig="aca_cycle1.cfg"
integration        =  "10s"
</source>
mapsize            =  "1arcmin"
maptype            =  "hex"
pointingspacing    = "9arcsec"     # this could also be specified in units of the primary beam e.g. "0.5PB"
</source>  


=====Calculate Visibilities=====
In nominal weather:
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"
pwv=0.6
antennalist        =  "ALMA;0.5arcsec"
totaltime          =  "3600s"
</source>
</source>


Run '''simobserve''', displaying graphics to the screen and also writing graphics output to files that can be found in the project directory, ''m51c''.
To cover the source as we've rescaled the pixel size, we'll need a 1 arcmin mosaic, and we'll let <tt>simalma</tt> calculate the pointings for us:
 
<source lang="python">
<source lang="python">
graphics          = "both"
mapsize="1arcmin"
simobserve()
</source>
</source>


----
----


====Simulate 12m total power observation====
=====What does it do?=====
[[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, on a square grid. CASA simulation tools can not simulate
[[Image:M51.alma_cycle1_3.observe.png|thumb|12m observation]]
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. 
The 12m array observation is simulated first -- <tt>simalma</tt> simply calls <tt>simobserve</tt> with your input parameters.
<tt>simobserve</tt> generates a graphic about the elevation and showing the dirty synthesized beam:


We need to change the mapping parameters to specify a square region a bit larger than the interferometric map.
The 12m-only visibilities are not currently imaged separately from the 7m visibilities, but this is an expected upgrade in a future release.  One could easily image the generated measurement set, which will be named according to the <tt>antennalist</tt> parameter above -- in this example, it is called <tt>m51.alma_cycle1_3.noisy.ms/</tt>.
<source lang="python">
integration        =  "10s"
mapsize            =  "1.3arcmin"
maptype            =  "square"
</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, 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.
<p>
----


<source lang="python">
[[Image:M51.aca_cycle1.skymodel.png|thumb|ACA hex map]]
obsmode            = "sd"
Next, the 7m ACA observation is simulated, with a second call to <tt>simobserve</tt>:
sdantlist          = "aca.tp.cfg"
<tt>simobserve</tt> follows the same conventions as the ALMA Observation Preparation Tool, and sets the mosaic pointings to cover the area requested -- it takes fewer 7m pointings to cover the region than it did 12m pointings.
sdant              = 0
refdate            = "2012/12/01"
totaltime          =  "2h"
</source>


Run '''simobserve''', displaying graphics to screen and to files.
It is useful to know that a version of the input sky model convolved to the ACA resolution is generated, in this example <tt>m51.aca_cycle1.skymodel.flat.regrid.conv/</tt>.  That image can be useful to better understand the simulation results.
----


<source lang="python">
Next, <tt>simobserve</tt> is called a third time to generate the total power image.  Again according to Cycle 1 conventions, the same size region as the 12m array mosaic are used, except that an extra pointing is added around the outside edge of the map, to guarantee that the total power map is larger than the interferometric mosaic (total power maps usually have noise and artifact at their edges).  Furthermore, a square raster pattern is used instead of the hexagonal patter of the interferometric array maps:
simobserve()
[[Image:M51.aca.tp.skymodel.png|thumb|TP map]]
</source>


----
----


====Simulate 7m ACA observation====
====Deconvolve the visibilities back into images====
[[Image:M51c.aca.i.skymodel.png|thumb|hexagonal map overplotted on sky model]]


Next we'll add an ACA mosaic, with its larger primary beam.
Next <tt>simalma</tt> uses '''simanalyze''' to combine the three measurement sets and create a single image.  


<source lang="python">
There are many ways to do this, and you may wish to discuss options with scientists at your ARC.  At time of writing, <tt>simalma</tt> concatenates the two sets of interferometric visibilities, images them, images the single dish spectra, and then uses the <tt>feather</tt> to combine the two images.
integration        =  "10s"
 
mapsize            = "1arcmin"
First, the total power image is generated using the ASAP gridding tools inside of CASA. <tt>simalma</tt> attempts to use the optimal gridding kernel to achieve maximum sensitivity and resolution of the single dish map, but the optimal parameters to use is an area of active investigation.
maptype            =  "hex"
pointingspacing    =  "15arcsec"
</source>  


We can specify an integral number of times to repeat the mosaic by setting ''totaltime'' to an integer string without units.
<font color="red">It is critical that the relative weights be set between the two different interferometric arrays.  Simulated data has weights=1, since the thermal noise is generated uniformly per baseline.  However, in reality the 7m baselines have lower sensitivity than the 12m baselines, and their weights must be decreased by that sensitivity ratio.  <tt>simalma</tt> uses the <tt>visweightscale</tt> parameter of <tt>concat</tt> to apply that lower weight of (7./12)**2 to the 7m visibilities.  If you wish to combine data manually, you must do this step yourself.</font>


<source lang="python">
The concatted visibilities are imaged, and diagnostic graphics with "concat" in their names are generated:
obsmode            = "int"
refdate            = "2012/12/02"
antennalist        =  "aca.i.cfg"
totaltime          =  "3"
</source>


Run '''simobserve''', displaying graphics to screen and to files.
[[Image:M51.concat.image.png|thumb|combined interferometric map]]


<source lang="python">
simobserve()
</source>


----
----


====Deconvolve the visibilities back into an image====
When combining the single dish and interferometric maps in the image plane using the <tt>feather</tt> task, one must use the interferometric map <it>without</it> the primary beam correction, and first multiply the total power map by the interferometric sensitivity image (".flux") -- this is ensure that noise effects are properly handled on the edges of each map, and do not grow artificially.  After running <tt>feather</tt>, the output is masked to 0.2 times the interferometric primary beam, since the total power map was created larger than the interferometric map on purpose, so the edges of the combined image do not contain any interferometric information.:
 
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====
[[Image:M51.combine.png|thumb|combined maps]]
[[Image:M51c.aca.i.analysis.png|thumb|Total power combined (with relatively low weight) with ACA]]


<source lang="python">
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()
</source>


====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]]
Here we explicitly have to set the ''modelimage'' to the result from the previous run.
<source lang="python">
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()
</source>


{{Checked 4.1.0}}
{{Checked 4.1.0}}

Latest revision as of 21:23, 21 March 2024


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

The simalma task

ALMA consists of the high resolution array of 12m antennas, the ALMA Compact Array of 7m antennas, and Total Power measurements using 12m antennas. One can simulate all of these in CASA.

One could use simobserve generate simulated observations for each component separately, and then combine 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 one of the interferometric simulations. Note that if you simulate total power and an interferometric observation ssimultaneously 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).

This "manual" method is described in the ACA_Simulation_(CASA_4.1) guide.

New in 4.1 is the simalma task, which takes one set of parameters describing the region of the sky to observe, and makes the appropriate calls to simobserve and simanalyze for the user. How this works is described here.

As of CASA 4.1, ALMA is still optimizing the algorithms for combining total power and interferometric data, so the parameters used here are likely to change.


Set simalma as current task

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

# Set simalma to default parameters
default("simalma")
# Our project name will be m51, and all simulation products will be placed in a subdirectory m51/
project="m51"
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 simalma will not modify your original input image. Rather, it will make a copy m51/m51.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="J2000 23h59m59.96s -34d59m59.50s"
incell="0.1arcsec"
inbright="0.004"
incenter="330.076GHz"
inwidth="50MHz"

Set up Observing Parameters

hexagonal mosaic overplotted on sky model

Based on the Cycle 1 capabilities, we would like to use array configuration number 3 which affords ~0.5 arcsec resolution:

antennalist="alma_cycle1_3.cfg"

We'll observe for a couple of hours (this is the 12m array observing time):

totaltime="7800s"

Following the Cycle 1 convention, simalma will observe 3 times longer with the 7m array and total power dishes.

acaratio=3.0
acaconfig="aca_cycle1.cfg"

In nominal weather:

pwv=0.6

To cover the source as we've rescaled the pixel size, we'll need a 1 arcmin mosaic, and we'll let simalma calculate the pointings for us:

mapsize="1arcmin"

What does it do?
12m observation

The 12m array observation is simulated first -- simalma simply calls simobserve with your input parameters. simobserve generates a graphic about the elevation and showing the dirty synthesized beam:

The 12m-only visibilities are not currently imaged separately from the 7m visibilities, but this is an expected upgrade in a future release. One could easily image the generated measurement set, which will be named according to the antennalist parameter above -- in this example, it is called m51.alma_cycle1_3.noisy.ms/.


ACA hex map

Next, the 7m ACA observation is simulated, with a second call to simobserve: simobserve follows the same conventions as the ALMA Observation Preparation Tool, and sets the mosaic pointings to cover the area requested -- it takes fewer 7m pointings to cover the region than it did 12m pointings.

It is useful to know that a version of the input sky model convolved to the ACA resolution is generated, in this example m51.aca_cycle1.skymodel.flat.regrid.conv/. That image can be useful to better understand the simulation results.


Next, simobserve is called a third time to generate the total power image. Again according to Cycle 1 conventions, the same size region as the 12m array mosaic are used, except that an extra pointing is added around the outside edge of the map, to guarantee that the total power map is larger than the interferometric mosaic (total power maps usually have noise and artifact at their edges). Furthermore, a square raster pattern is used instead of the hexagonal patter of the interferometric array maps:

TP map

Deconvolve the visibilities back into images

Next simalma uses 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. At time of writing, simalma concatenates the two sets of interferometric visibilities, images them, images the single dish spectra, and then uses the feather to combine the two images.

First, the total power image is generated using the ASAP gridding tools inside of CASA. simalma attempts to use the optimal gridding kernel to achieve maximum sensitivity and resolution of the single dish map, but the optimal parameters to use is an area of active investigation.

It is critical that the relative weights be set between the two different interferometric arrays. Simulated data has weights=1, since the thermal noise is generated uniformly per baseline. However, in reality the 7m baselines have lower sensitivity than the 12m baselines, and their weights must be decreased by that sensitivity ratio. simalma uses the visweightscale parameter of concat to apply that lower weight of (7./12)**2 to the 7m visibilities. If you wish to combine data manually, you must do this step yourself.

The concatted visibilities are imaged, and diagnostic graphics with "concat" in their names are generated:

combined interferometric map



When combining the single dish and interferometric maps in the image plane using the feather task, one must use the interferometric map <it>without</it> the primary beam correction, and first multiply the total power map by the interferometric sensitivity image (".flux") -- this is ensure that noise effects are properly handled on the edges of each map, and do not grow artificially. After running feather, the output is masked to 0.2 times the interferometric primary beam, since the total power map was created larger than the interferometric map on purpose, so the edges of the combined image do not contain any interferometric information.:

combined maps


Last checked on CASA Version 4.1.0.