First Look at Image Analysis: Difference between revisions

From CASA Guides
Jump to navigationJump to search
Tashton (talk | contribs)
No edit summary
Hsheets (talk | contribs)
Tag: New redirect
 
(5 intermediate revisions by 3 users not shown)
Line 1: Line 1:
The first few tutorials provided us with a spectral line data cube and a continuum image.  Next we'd like to
#REDIRECT [[First_Look_at_Image_Analysis_CASA_6.5.4]]
understand some of the properties of the images that we have produced.
 
In case you don't have the images from the first tutorials handy, let's start by copying the image files from the working/ directory to our current directory:
 
<source lang="python">
# In CASA
os.system("rm -rf sis14_twhya_cont.image")
os.system("cp -r ../working/sis14_twhya_cont.image .")
os.system("rm -rf twhya_n2hp.image")
os.system("cp -r ../working/twhya_n2hp.image .")
</source>
 
Recall that we orient ourselves with measurement sets (UV data) using the {{listobs}} task.
For images, an analogous way to get basic header information is to use the {{imhead}} task.
 
<source lang="python">
# In CASA
imhead("sis14_twhya_cont.image")
imhead("twhya_n2hp.image")
</source>
 
CASA prints the image header information to the logger:
<pre style="background-color: #fffacd;">
imhead ##########################################
imhead ##### Begin Task: imhead            #####
imhead imhead(imagename="sis14_twhya_cont.image",mode="summary",hdkey="",hdvalue="",verbose=False)
ImageMetaData  
ImageMetaData Image name      : sis14_twhya_cont.image
ImageMetaData Object name      : TW Hya
ImageMetaData Image type      : PagedImage
ImageMetaData Image quantity  : Intensity
ImageMetaData Pixel mask(s)    : mask0
ImageMetaData Region(s)        : None
ImageMetaData Image units      : Jy/beam
ImageMetaData Restoring Beam  : 0.830341 arcsec, 0.683704 arcsec, -66.1199 deg
ImageMetaData  
ImageMetaData Direction reference : J2000
ImageMetaData Spectral  reference : LSRK
ImageMetaData Velocity  type      : RADIO
ImageMetaData Rest frequency      : 3.72637e+11 Hz
ImageMetaData Pointing center    :  11:01:51.796000  -34.42.17.366000
ImageMetaData Telescope          : ALMA
ImageMetaData Observer            : cqi
ImageMetaData Date observation    : 2012/11/19/07:56:27
ImageMetaData Telescope position: [2.22514e+06m, -5.44031e+06m, -2.48103e+06m] (ITRF)
ImageMetaData  
ImageMetaData Axis Coord Type      Name            Proj Shape Tile  Coord value at pixel        Coord incr Units
ImageMetaData ----------------------------------------------------------------------------------------------------
ImageMetaData 0    0    Direction Right Ascension  SIN  250  250  11:01:51.796  125.00    -1.000000e-01 arcsec
ImageMetaData 1    0    Direction Declination      SIN  250  250 -34.42.17.366  125.00      1.000000e-01 arcsec
ImageMetaData 2    1    Stokes    Stokes                    1    1            I
ImageMetaData 3    2    Spectral  Frequency                1    1  3.72637e+11    0.00 2.34445114878e+08 Hz
ImageMetaData                     Velocity                                    0    0.00    -1.886149e+02 km/s
imhead ##### End Task: imhead              #####
imhead ##########################################
 
imhead ##########################################
imhead ##### Begin Task: imhead            #####
imhead imhead(imagename="twhya_n2hp.image",mode="summary",hdkey="",hdvalue="",verbose=False)
ImageMetaData  
ImageMetaData Image name      : twhya_n2hp.image
ImageMetaData Object name      : TW Hya
ImageMetaData Image type      : PagedImage
ImageMetaData Image quantity  : Intensity
ImageMetaData Pixel mask(s)    : mask0
ImageMetaData Region(s)        : None
ImageMetaData Image units      : Jy/beam
ImageMetaData Restoring Beam  : 0.750766 arcsec, 0.598023 arcsec, -59.397 deg
ImageMetaData  
ImageMetaData Direction reference : J2000
ImageMetaData Spectral  reference : LSRK
ImageMetaData Velocity  type      : RADIO
ImageMetaData Rest frequency      : 3.72672e+11 Hz
ImageMetaData Pointing center    :  11:01:51.796000  -34.42.17.366000
ImageMetaData Telescope          : ALMA
ImageMetaData Observer            : cqi
ImageMetaData Date observation    : 2012/11/19/07:56:27
ImageMetaData Telescope position: [2.22514e+06m, -5.44031e+06m, -2.48103e+06m] (ITRF)
ImageMetaData  
ImageMetaData Axis Coord Type      Name            Proj Shape Tile  Coord value at pixel      Coord incr Units
ImageMetaData --------------------------------------------------------------------------------------------------
ImageMetaData 0    0    Direction Right Ascension  SIN  250  125  11:01:51.796  125.00  -8.000000e-02 arcsec
ImageMetaData 1    0    Direction Declination      SIN  250  50 -34.42.17.366  125.00    8.000000e-02 arcsec
ImageMetaData 2    1    Stokes    Stokes                    1    1            I
ImageMetaData 3    2    Spectral  Frequency                15    5  3.726725e+11    0.00 -6.21550810e+05 Hz
ImageMetaData                     Velocity                                    0    0.00    5.000000e-01 km/s
imhead ##### End Task: imhead              #####
imhead ##########################################
</pre>
 
== Statistics ==
 
Often we'd like to measure basic image
statistics and fluxes. Recall you can do these operations pretty easily in the viewer, by
dragging out a box and then right-double clicking inside it to get statistics for that region. You can also get statistics using the {{imstat}} task.
 
<source lang="python">
# In CASA
my_stats = imstat("twhya_n2hp.image", chans="0~4")
</source>
 
The {{imstat}} task returns a python dictionary, and you can examine the contents as follows:
 
<source lang="python">
# In CASA
print(my_stats)
print(my_stats['rms'])
</source>
 
It's useful to see that the RMS is about 26 mJy/beam in the
line cube.
 
You can also just print the dictionary to the terminal screen:
 
<source lang="python">
# In CASA
imstat("twhya_n2hp.image", chans="0~4")
</source>
 
This task is useful for measuring basic source characteristics. For example,
calculate the statistics for a box encompassing the disk - you see that the
integrated flux is about 1.8 Jy.
 
<source lang="python">
# In CASA
imstat("sis14_twhya_cont.image", box="100,100,150,150")
my_stats = imstat("sis14_twhya_cont.image", box="100,100,150,150")
print(my_stats['flux'])
</source>
 
Alternatively, a box defined off the disk will give image noise statistics.
 
<source lang="python">
# In CASA
imstat("sis14_twhya_cont.image", box="25,150,225,200")
</source>
 
== Moments ==
 
For the spectral line cube, it's very useful to "collapse" the cube in
various ways to analyze the emission. The {{immoments}} task lets you do
this.
 
The first few "moments" are defined as follows:
<pre>
moments = -1 : mean value of the spectral coordinate
moments =  0 : integrated value of the spectral coordinate
moments =  1 : intensity weighted coordinate; traditionally used to get 'velocity fields'
moments =  2 : intensity weighted dispersion of the coordinate; traditionally used to get "velocity dispersion"
</pre>
 
A more complete description of image moments is given [http://casa.nrao.edu/docs/CasaRef/image.moments.html here in the CASA reference manual].
 
<figure id="Imaging-tutorial-analysis-mom0_5.7.png">
[[File:Imaging-tutorial-analysis-mom0_5.7.png|thumb|<caption>The zero moment of the N2H+ image.</caption>]]
</figure>
 
Let's make a moment 0 image clipped at ~1 sigma:
 
<source lang="python">
# In CASA
os.system("rm -rf sis14_twhya_n2hp.mom0")
immoments("twhya_n2hp.image",
          outfile="sis14_twhya_n2hp.mom0",
          includepix=[20e-3,100],
          chans="4~12",
          moments=0)
</source>
 
Make a moment 1 image clipped at ~2 sigma:
 
<source lang="python">
# In CASA
os.system("rm -rf sis14_twhya_n2hp.mom1")
immoments("sis14_twhya_n2hp.image",
          outfile="sis14_twhya_n2hp.mom1",
          includepix=[40e-3,100],
          chans="4~12",
          moments=1)
</source>
 
At this point we have a few really neat things to see: first the
line shows a hole in the middle. Overlay it on the dust (continuum)
disk using the viewer and you can notice that they align but with the N2H+
existing only outside a certain radius - in this case the "snow line".
Also have a look at the velocity field to see the rotation of the
disk.
 
The {{imview}} task has some command-line scripting capability.  For example, here we show how to overlay the line moment
0 (as a contour plot) on the continuum.
 
<figure id="Imaging-tutorial-analysis-overlay_5.7.png">
[[File:Imaging-tutorial-analysis-overlay_5.7.png|thumb|<caption>The N2H+ moment zero map as contours, overlaid on the continuum emission of TW Hydra.  The plot demonstrates the "show line".</caption>]]
</figure>
 
<source lang="python">
# In CASA
imview(raster={'file': 'sis14_twhya_cont.image',
'range': [-0.01,0.5]},
contour={'file': 'sis14_twhya_n2hp.mom0',
'levels': [0.5,0.6,0.7,0.8] })
</source>
 
Here, we are not using the primary beam corrected images and cubes when creating our moment maps. Remember that it's often very convenient to work in images before primary beam correction because the noise is the same across the field (e.g., this is a clean data set to search for signal) but it's very important to remember to apply this correction before calculating fluxes or intensities for science.
 
== Export FITS images ==
 
CASA is great (of course) but you will ultimately want to export
your data to a common format to analyze in other programs, share
with other astronomers, or archive. It's easy to export images from
CASA's image format to FITS images via the {{exportfits}} command.
 
<source lang="python">
# In CASA
exportfits(imagename="sis14_twhya_cont.image",
          fitsimage="twhya_cont.fits",
          overwrite=True)
</source>
 
For the cube we want to specify additionally that the frequency axis
will be written out as velocity.
 
<source lang="python">
# In CASA
exportfits(imagename="twhya_n2hp.image",
          fitsimage="twhya_n2hp.fits",
          velocity=True,
          overwrite=True)
</source>

Latest revision as of 21:18, 12 October 2023