ACA Simulation (CASA 3.3): Difference between revisions

From CASA Guides
Jump to navigationJump to search
Rindebet (talk | contribs)
No edit summary
 
(30 intermediate revisions by 3 users not shown)
Line 3: Line 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]].
 
__NOTOC__
== ALMA 12m + ACA + Total Power ==
== ALMA 12m + ACA + Total Power ==


Line 10: Line 10:
Measurement Sets to run sim_analyze.  Total power observations can be simulated either in an independent run of sim_observe, or along with an interferometric simulation.  Note that if you simulate total power and an interferometric observation simultaneously with sim_observe, 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 12m ALMA mosaic).
Measurement Sets to run sim_analyze.  Total power observations can be simulated either in an independent run of sim_observe, or along with an interferometric simulation.  Note that if you simulate total power and an interferometric observation simultaneously with sim_observe, 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 12m ALMA mosaic).


*Simdata version for CASA 3.3
=====Set sim_observe as current task=====
__NOTOC__
Reset all parameters to default, and then set the project name to m51c
<br>
====Explanation of the script====
 
=====Set sim_observe as current task and reset all parameters=====
<source lang="python">
<source lang="python">
# Setting everything in sim_observe to original defaults
# Set sim_observe to default parameters
default("sim_observe")
default("sim_observe")
# Our project name will be m51c, and all simulation products will be placed in a subdirectory 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.
sim_observe will make a copy m51c/m51c.skymodel, and not modify your input image.
<source lang="python">
<source lang="python">
# We'll use an Halpha image of M51 as an illustrative model
# Model sky = Halpha image of M51  
# Make sure you are running CASA in the same directory as this file
os.system('curl http://casaguides.nrao.edu/images/3/3f/M51ha.fits.txt -f -o M51ha.fits.txt')
modelimage         =  "input50pc_672GHz.fits"
skymodel         =  "M51ha.fits.txt"
</source> 
=====Image coordinate system can be verified=====
<source lang="python">
# This reports image header parameters in the Log Messages window
imhead("input50pc_672GHz.fits")
</source>   
</source>   
=====Image center can be identified=====
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,
* 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
* 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.
<source lang="python">
<source lang="python">
# These are tools found in the CASA toolkit
# Set model image parameters:
# They are very useful, but the interface is not as straightforward as the tasks
indirection="B1950 23h59m59.96s -34d59m59.50s"
# You can find the tool reference manual here: http://casa.nrao.edu/docs/CasaRef/CasaRef.html
incell="0.1arcsec"
# The following command is used to open an image (this is part of the image analysis toolkit)
inbright="0.004"
# Whenever data are being viewed/manipulated by tools, what is being operated on needs to be explicitly opened
incenter="330.076GHz"
# and closed (i.e. an image, a table, etc.)
inwidth="50MHz"
ia.open("input50pc_672GHz.fits")
#  Out[9]: True
#
# Reports the length of each axis in the opened image
ia.shape()
#  Out[11]: [257L, 257L, 1L, 1L]
#
# This command converts from pixel (our source file) to world coordinates (something usable by simdata)
ia.toworld([128.5,128.5])
#  Out[12]:
#{'numeric': array([  4.71239120e+00,  -4.01423802e-01,  1.00000000e+00,
#        6.72000001e+11])}
#
# Formats the coordinate just converted into hms
qa.formxxx("4.71239120rad",format='hms',prec=5)
#  Out[13]: '18:00:00.03052'
#
# Formats one of the other coordinates into dms
qa.formxxx("-0.401423802rad",format='dms',prec=5)
#  Out[14]: '-022.59.59.602743'
#
# Final housekeeping by closing the image tool
# The image tool will now be detached from the image
ia.close()
</source>
</source>


=====Brightness scale can be viewed with 'imstat' task=====
<source lang="python">
# Default parameters are adequate for this
imstat("input50pc_672GHz.fits")
# ...
#  'max': array([  6.52469971e-05]),
# ...
# that's 0.0652 mJy/pixel.
</source> 
=====Let's call our project psim2=====
<source lang="python">
# This defines the root prefix for any output files from simdata
project            =  "psim2"
</source>


=====We'll leave the sky model the way it is: simdata will create psim2.skymodel CASA image since this model is a fits file, and most but not all of CASA routines can operate directly on fits=====
----
<source lang="python">
 
skymodel          =  "input50pc_672GHz.fits"
====Simulate 12m interferometric observation====
</source> 
[[Image:M51c.ALMA_0.5arcsec.skymodel.png|thumb|hexagonal mosaic overplotted on sky model]]
=====We need to decide where to point the telescopeThe image is 2/3 arcsec in size, so we only need one pointing. We could put that in a text file ourself, or let simdata create the ascii pointing file for us.=====
 
We'll begin with the 12m ALMA array observation, and have sim_observe calculate a hexagonal mosaic of pointings.
 
The default interface for sim_observe provides an integration parameter, which is the dwell time at each mosaic pointing -- we'll set that to 10sA 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, sim_observe would only generate one measurement per 5min scan. While 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 coverage.  If you wish to simulate the more realistic case such as 5 min scans with 5s integrations, please see [[Complex pointingtable]] for a guide to doing that with sim_observe (Its not hard, it just takes two calls to the task instead of one).
 
<source lang="python">
<source lang="python">
# have sim_observe calculate mosaic pointing locations:
setpointings      =  True
setpointings      =  True
direction          =  "J2000 18h00m00.031s -22d59m59.6s"
integration        =  "10s"
mapsize            =  "0.76arcsec"
mapsize            =  "1arcmin"
maptype            =  "hex"
pointingspacing    =  "9arcsec"      # this could also be specified in units of the primary beam e.g. "0.5PB"
</source>  
</source>  
=====The default pointingspacing is fine: we'll only fit one pointing in the small mapsize the default calculation maptype hexagonal is ok too since only one will fit anyway.=====


=====We do want to calculate visibilities in a measurement set: let's do a 20 min snapshot observation using the "out20" ALMA antenna configuration:=====
=====Calculate Visibilities=====
sim_observe can determine what array configuration to use, if you provide a desired resolution or synthesized beam size. 
<source lang="python">
<source lang="python">
observe            =  True
observe            =  True
totaltime          =  "1200s"
antennalist        =  "ALMA;0.5arcsec"
totaltime          =  "3600s"
</source>
</source>


=====Use appropriate antenna configurations based on desired angular resolution (configuration 20 - alma.out20.cfg in this case - is the largest "compact" configuration)=====
run sim_observe, displaying graphics to screen and file (files can be found in the project directory, e.g. m51c)
 
<source lang="python">
<source lang="python">
# It might be helpful to confirm the alma.out20.cfg file exists in the path defined below
graphics          = "both"
# If you have a problem, this might be the first thing to check, if you haven't already
go
repodir=os.getenv("CASAPATH").split(' ')[0]
antennalist        =  repodir+"/data/alma/simmos/alma.out20.cfg"
</source>
</source>


=====Deconvolve the visibilities back into an image=====
 
----
 
====Simulate 12m total power observation====
[[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
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 pointing calculation to make it square and a bit larger than the interferometric map.
<source lang="python">
<source lang="python">
image              True
integration        "10s"
vis                =  "$project.ms"
mapsize            =  "1.3arcmin"
imsize            [192, 192]  
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, 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)
 
<source lang="python">
sdantlist          = "aca.tp.cfg"
sdant              = 0
refdate            = "2012/12/01"
# SET interferometric antennalist to default, or else we'll redo the interferometric simulation too
antennalist        = ""
totaltime          =  "2h"
</source>
</source>


=====Specify number of iteration of cleaning task with proper threshold and weighting=====
run sim_observe, displaying graphics to screen and file
 
<source lang="python">
<source lang="python">
niter              =  10000
go
threshold          =  "1e-7Jy"
weighting          =  "natural"   
</source>
</source>


=====We'd like to calculate a difference and fidelity image, and see some diagnostics:=====
 
----
 
====Simulate 7m ACA observation====
[[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.
 
<source lang="python">
<source lang="python">
analyze           =  True  
integration        =  "10s"
mapsize            =  "1arcmin"
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.
 
<source lang="python">
# remember to set this back to empty
sdantlist          = ""
refdate           = "2012/12/02"
antennalist        "aca.i.cfg"
totaltime          = "3"
</source>
</source>


=====And see the array but not the UV coverage:=====
run sim_observe, displaying graphics to screen and file
 
<source lang="python">
<source lang="python">
showarray          =  True
go
showuv            =  False 
</source>
</source>


=====Plot both to the screen and the png files with lots of messages:=====
 
----
 
====Deconvolve the visibilities back into an image====
 
We use sim_analyze to take 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 optimal.  Please discuss with your ARC contact scientist if you have ALMA data now, or wait for additional recommendations to be posted here over time.
 
* 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.
* 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.
 
sim_analyze, if given a total power and interferometric measurement set, 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:M51c.aca.i.analysis.png|thumb|Total power combined (with relatively low weight) with ACA]]
 
<source lang="python">
<source lang="python">
graphics          =  "both"
default("sim_analyze")
verbose           =  True
project            =  "m51c"
overwrite = True
vis                =  '$project.aca.i.ms,$project.aca.tp.sd.ms' 
imsize            =  [512,512]
cell              =  '0.2arcsec'
analyze           =  True
showpsf            =  False
showresidual      = False
showconvolved      =  True
go()
</source>
</source>


===Run simdata===
 
 
 
 
====Next add the 12m interferometric data====
[[Image:M51c.ALMA_0.5arcsec.analysis.png|thumb|TP+ACA combined (with relatively low weight) with 12m ALMA]]
Here we explicitly have to set the modelimage to the previous output.
 
<source lang="python">
<source lang="python">
# This commands CASA to execute simdata
default("sim_analyze")
simdata()
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
go()
</source>
</source>
<br>
*Output results:


{| style="border:1px solid #3366FF; " cellspacing=2
|Input:<br> [[File:Psim2.skymodel.png|300px]]
|Predict:<br> [[File:Psim2.predict.png|300px]]
|-
|Image:<br> [[File:Psim2.image.png|300px]]
|Analyze:<br> [[File:Psim2.analysis.png|300px]]
|}


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

Latest revision as of 22:51, 26 March 2024

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

Multiple configurations, or an observation in one 12m configuration plus one ACA observation, or other combinations, can be simulated by running sim_observe for each, and then combining the Measurement Sets to run sim_analyze. Total power observations can be simulated either in an independent run of sim_observe, or along with an interferometric simulation. Note that if you simulate total power and an interferometric observation simultaneously with sim_observe, 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 12m ALMA mosaic).

Set sim_observe as current task

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

# Set sim_observe to default parameters
default("sim_observe")
# 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 an example model sky. Download File:M51ha.fits.txt and place in your working directory, or use the curl command in the script.

sim_observe will make a copy m51c/m51c.skymodel, and not modify your input image.

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

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,
  • 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
  • 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.
# 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 with the 12m ALMA array observation, and have sim_observe calculate a hexagonal mosaic of pointings.

The default interface for sim_observe 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, sim_observe would only generate one measurement per 5min scan. While 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 coverage. If you wish to simulate the more realistic case such as 5 min scans with 5s integrations, please see Complex pointingtable for a guide to doing that with sim_observe (Its not hard, it just takes two calls to the task instead of one).

# have sim_observe 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

sim_observe can determine what array configuration to use, if you provide a desired resolution or synthesized beam size.

observe            =  True
antennalist        =  "ALMA;0.5arcsec"
totaltime          =  "3600s"

run sim_observe, displaying graphics to screen and file (files can be found in the project directory, e.g. m51c)

graphics           =  "both"
go



Simulate 12m total power observation

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 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 pointing calculation to make it square and 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, 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)

sdantlist          = "aca.tp.cfg"
sdant              = 0
refdate            = "2012/12/01"
# SET interferometric antennalist to default, or else we'll redo the interferometric simulation too
antennalist        =  ""
totaltime          =  "2h"

run sim_observe, displaying graphics to screen and file

go



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.

# remember to set this back to empty
sdantlist          = ""
refdate            = "2012/12/02"
antennalist        =  "aca.i.cfg"
totaltime          =  "3"

run sim_observe, displaying graphics to screen and file

go



Deconvolve the visibilities back into an image

We use sim_analyze to take 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 optimal. Please discuss with your ARC contact scientist if you have ALMA data now, or wait for additional recommendations to be posted here over time.

  • 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.
  • 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.

sim_analyze, if given a total power and interferometric measurement set, 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("sim_analyze")
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
go()



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 previous output.

default("sim_analyze")
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
go()


Simulating Observations in CASA