First Look at Image Analysis CASA 6: Difference between revisions
From CASA Guides
Jump to navigationJump to search
No edit summary |
No edit summary |
||
| (6 intermediate revisions by 3 users not shown) | |||
| Line 1: | Line 1: | ||
''This guide is written for CASA 6.2.1.7 and uses Python 3.'' | |||
The first few tutorials provided us with a spectral line data cube and a continuum image. Next we'd like to | The first few tutorials provided us with a spectral line data cube and a continuum image. Next we'd like to | ||
| Line 20: | Line 20: | ||
# In CASA | # In CASA | ||
imhead("sis14_twhya_cont.image") | imhead("sis14_twhya_cont.image") | ||
imhead(" | imhead("twhya_n2hp.image") | ||
</source> | </source> | ||
| Line 33: | Line 33: | ||
ImageMetaData Image type : PagedImage | ImageMetaData Image type : PagedImage | ||
ImageMetaData Image quantity : Intensity | ImageMetaData Image quantity : Intensity | ||
ImageMetaData Pixel mask(s) : | ImageMetaData Pixel mask(s) : mask0 | ||
ImageMetaData Region(s) : None | ImageMetaData Region(s) : None | ||
ImageMetaData Image units : Jy/beam | ImageMetaData Image units : Jy/beam | ||
ImageMetaData Restoring Beam : 0. | ImageMetaData Restoring Beam : 0.830341 arcsec, 0.683704 arcsec, -66.1199 deg | ||
ImageMetaData | ImageMetaData | ||
ImageMetaData Direction reference : J2000 | ImageMetaData Direction reference : J2000 | ||
ImageMetaData Spectral reference : LSRK | ImageMetaData Spectral reference : LSRK | ||
ImageMetaData Velocity type : RADIO | ImageMetaData Velocity type : RADIO | ||
ImageMetaData Rest frequency : 3.72637e+11 Hz | ImageMetaData Rest frequency : 3.72637e+11 Hz | ||
| Line 45: | Line 45: | ||
ImageMetaData Telescope : ALMA | ImageMetaData Telescope : ALMA | ||
ImageMetaData Observer : cqi | ImageMetaData Observer : cqi | ||
ImageMetaData Date observation : 2012/11/19/07: | ImageMetaData Date observation : 2012/11/19/07:56:27 | ||
ImageMetaData Telescope position: [2. | ImageMetaData Telescope position: [2.22514e+06m, -5.44031e+06m, -2.48103e+06m] (ITRF) | ||
ImageMetaData | ImageMetaData | ||
ImageMetaData Axis Coord Type Name Proj Shape Tile Coord value at pixel Coord incr Units | ImageMetaData Axis Coord Type Name Proj Shape Tile Coord value at pixel Coord incr Units | ||
ImageMetaData ---------------------------------------------------------------------------------------------------- | ImageMetaData ---------------------------------------------------------------------------------------------------- | ||
ImageMetaData 0 0 Direction Right Ascension SIN | 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 | 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 2 1 Stokes Stokes 1 1 I | ||
ImageMetaData 3 2 Spectral Frequency 1 1 3.72637e+11 0.00 2. | ImageMetaData 3 2 Spectral Frequency 1 1 3.72637e+11 0.00 2.34445114878e+08 Hz | ||
ImageMetaData Velocity 0 0.00 -1. | ImageMetaData Velocity 0 0.00 -1.886149e+02 km/s | ||
imhead ##### End Task: imhead ##### | imhead ##### End Task: imhead ##### | ||
imhead ########################################## | imhead ########################################## | ||
| Line 60: | Line 60: | ||
imhead ########################################## | imhead ########################################## | ||
imhead ##### Begin Task: imhead ##### | imhead ##### Begin Task: imhead ##### | ||
imhead imhead(imagename=" | imhead imhead(imagename="twhya_n2hp.image",mode="summary",hdkey="",hdvalue="",verbose=False) | ||
ImageMetaData | ImageMetaData | ||
ImageMetaData Image name : | ImageMetaData Image name : twhya_n2hp.image | ||
ImageMetaData Object name : TW Hya | ImageMetaData Object name : TW Hya | ||
ImageMetaData Image type : PagedImage | ImageMetaData Image type : PagedImage | ||
| Line 69: | Line 69: | ||
ImageMetaData Region(s) : None | ImageMetaData Region(s) : None | ||
ImageMetaData Image units : Jy/beam | ImageMetaData Image units : Jy/beam | ||
ImageMetaData Restoring Beam : 0. | ImageMetaData Restoring Beam : 0.750766 arcsec, 0.598023 arcsec, -59.397 deg | ||
ImageMetaData | ImageMetaData | ||
ImageMetaData Direction reference : J2000 | ImageMetaData Direction reference : J2000 | ||
| Line 79: | Line 79: | ||
ImageMetaData Observer : cqi | ImageMetaData Observer : cqi | ||
ImageMetaData Date observation : 2012/11/19/07:56:27 | ImageMetaData Date observation : 2012/11/19/07:56:27 | ||
ImageMetaData Telescope position: [2. | ImageMetaData Telescope position: [2.22514e+06m, -5.44031e+06m, -2.48103e+06m] (ITRF) | ||
ImageMetaData | ImageMetaData | ||
ImageMetaData Axis Coord Type Name Proj Shape Tile Coord value at pixel Coord incr Units | ImageMetaData Axis Coord Type Name Proj Shape Tile Coord value at pixel Coord incr Units | ||
| Line 86: | Line 86: | ||
ImageMetaData 1 0 Direction Declination SIN 250 50 -34.42.17.366 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 2 1 Stokes Stokes 1 1 I | ||
ImageMetaData 3 2 Spectral Frequency 15 5 | 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 | ImageMetaData Velocity 0 0.00 5.000000e-01 km/s | ||
imhead ##### End Task: imhead ##### | imhead ##### End Task: imhead ##### | ||
| Line 100: | Line 100: | ||
<source lang="python"> | <source lang="python"> | ||
# In CASA | # In CASA | ||
my_stats = imstat(" | my_stats = imstat("twhya_n2hp.image", chans="0~4") | ||
</source> | </source> | ||
| Line 111: | Line 111: | ||
</source> | </source> | ||
It's useful to see that the RMS is about | It's useful to see that the RMS is about 26 mJy/beam in the | ||
line cube. | line cube. | ||
| Line 118: | Line 118: | ||
<source lang="python"> | <source lang="python"> | ||
# In CASA | # In CASA | ||
imstat(" | imstat("twhya_n2hp.image", chans="0~4") | ||
</source> | </source> | ||
This task is useful for measuring basic source characteristics. For example, | This task is useful for measuring basic source characteristics. For example, | ||
calculate the statistics for a box encompassing the disk - you see that the | calculate the statistics for a box encompassing the disk - you see that the | ||
integrated flux is about 1. | integrated flux is about 1.8 Jy. | ||
<source lang="python"> | <source lang="python"> | ||
| Line 155: | Line 155: | ||
A more complete description of image moments is given [http://casa.nrao.edu/docs/CasaRef/image.moments.html here in the CASA reference manual]. | A more complete description of image moments is given [http://casa.nrao.edu/docs/CasaRef/image.moments.html here in the CASA reference manual]. | ||
[[File:Imaging-tutorial-analysis- | [[File:Imaging-tutorial-analysis-mom0_6.1.png|thumb|<caption>Figure 1: The zero moment of the N2H+ image.</caption>]] | ||
Let's make a moment 0 image clipped at ~1 sigma: | Let's make a moment 0 image clipped at ~1 sigma: | ||
| Line 164: | Line 164: | ||
# In CASA | # In CASA | ||
os.system("rm -rf sis14_twhya_n2hp.mom0") | os.system("rm -rf sis14_twhya_n2hp.mom0") | ||
immoments(" | immoments("twhya_n2hp.image", | ||
outfile="sis14_twhya_n2hp.mom0", | outfile="sis14_twhya_n2hp.mom0", | ||
includepix=[20e-3,100], | includepix=[20e-3,100], | ||
| Line 176: | Line 176: | ||
# In CASA | # In CASA | ||
os.system("rm -rf sis14_twhya_n2hp.mom1") | os.system("rm -rf sis14_twhya_n2hp.mom1") | ||
immoments(" | immoments("twhya_n2hp.image", | ||
outfile="sis14_twhya_n2hp.mom1", | outfile="sis14_twhya_n2hp.mom1", | ||
includepix=[40e-3,100], | includepix=[40e-3,100], | ||
| Line 193: | Line 193: | ||
0 (as a contour plot) on the continuum. | 0 (as a contour plot) on the continuum. | ||
[[File:Imaging-tutorial-analysis-overlay_6.1.png|thumb|<caption>Figure 2: The N2H+ moment zero map as contours, overlaid on the continuum emission of TW Hydra. The plot demonstrates the "show line".</caption>]] | |||
[[File:Imaging-tutorial-analysis- | |||
<source lang="python"> | <source lang="python"> | ||
| Line 226: | Line 225: | ||
<source lang="python"> | <source lang="python"> | ||
# In CASA | # In CASA | ||
exportfits(imagename=" | exportfits(imagename="twhya_n2hp.image", | ||
fitsimage="twhya_n2hp.fits", | fitsimage="twhya_n2hp.fits", | ||
velocity=True, | velocity=True, | ||
overwrite=True) | overwrite=True) | ||
</source> | </source> | ||
Latest revision as of 18:28, 6 April 2022
This guide is written for CASA 6.2.1.7 and uses Python 3.
The first few tutorials provided us with a spectral line data cube and a continuum image. Next we'd like to 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_data directory to our current directory:
# 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 .")
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.
# In CASA
imhead("sis14_twhya_cont.image")
imhead("twhya_n2hp.image")
CASA prints the image header information to the logger:
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 ##########################################
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.
# In CASA
my_stats = imstat("twhya_n2hp.image", chans="0~4")
The imstat task returns a python dictionary, and you can examine the contents as follows:
# In CASA
print(my_stats)
print(my_stats['rms'])
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:
# In CASA
imstat("twhya_n2hp.image", chans="0~4")
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.
# 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'])
Alternatively, a box defined off the disk will give image noise statistics.
# In CASA
imstat("sis14_twhya_cont.image", box="25,150,225,200")
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:
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"
A more complete description of image moments is given here in the CASA reference manual.
Let's make a moment 0 image clipped at ~1 sigma:
# 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)
Make a moment 1 image clipped at ~2 sigma:
# In CASA
os.system("rm -rf sis14_twhya_n2hp.mom1")
immoments("twhya_n2hp.image",
outfile="sis14_twhya_n2hp.mom1",
includepix=[40e-3,100],
chans="4~12",
moments=1)
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.
# 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] })
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.
# In CASA
exportfits(imagename="sis14_twhya_cont.image",
fitsimage="twhya_cont.fits",
overwrite=True)
For the cube we want to specify additionally that the frequency axis will be written out as velocity.
# In CASA
exportfits(imagename="twhya_n2hp.image",
fitsimage="twhya_n2hp.fits",
velocity=True,
overwrite=True)
