How to rotate and slice a cube for pV diagrams: Difference between revisions

From CASA Guides
Jump to navigationJump to search
No edit summary
No edit summary
 
(2 intermediate revisions by the same user not shown)
Line 3: Line 3:
Here's a recipe:
Here's a recipe:


1) rotate the input cube by a given position angle (to align the slice in the vertical direction):
1) find out the axes of your cube.
 
<source lang="python">
imhead('original.image')
</source>
 
and in the example here, the logger spits out:
 
<pre style="background-color: #fffacd;">
Axis Coord Type      Name            Proj Shape Tile  Coord value at pixel      Coord incr Units
--------------------------------------------------------------------------------------------------
0    0    Direction Right Ascension  SIN  300  60  09:47:57.382  150.00  -4.000000e-01 arcsec
1    0    Direction Declination      SIN  300  60 +13.16.40.660  150.00    4.000000e-01 arcsec
2    1    Stokes    Stokes                    1    1            I
3    2    Spectral  Frequency                54    9  3.63945e+10    0.00 -1.25014600e+05 Hz
                    Velocity                              -17.7893    0.00    1.029965e+00 km/s
</pre>
 
Note that axis 0 is RA, axis 1 ins Dec, Stokes is axis 2, and frequency/velocity is the last axis 3.
 
2) rotate the input cube by a given position angle (to align the slice in the vertical or horizontal direction):


<source lang="python">
<source lang="python">
Line 10: Line 30:


#rotate by the position angle and write out a new, rotated cube
#rotate by the position angle and write out a new, rotated cube
ia.rotate(pa='45deg',outfile='rotated.image')
ia.rotate(pa='40deg',outfile='rotated.image')


# close the image
# close the image
Line 19: Line 39:
You now have a rotated cube 'rotated.image'
You now have a rotated cube 'rotated.image'


2) now create a box region by finding the bottom left corner (blc) and top right corner (trc) pixel coordinates. Attach this to a variable, here 'mybox'
Load and inspect with the viewer
 
<source lang="python">
viewer 'rotated.image'
</source>
 
 
 
3) now create a box region by finding the bottom left corner (blc) and top right corner (trc) pixel coordinates. Attach this to a variable, here 'mybox'. Note that the last two coords are Stokes (we don't touch it, so running from 0 to 0) and the velocity which we select to be channels 3 through 50:


<source lang="python">
<source lang="python">
mybox=rg.box(blc=[65,152],trc=[215,187])
mybox=rg.box(blc=[99,147,0,3],trc=[205,179,0,50])
</source>
</source>


3) bin along one of the axes within that region, drop the degenerate axis, write to a new file 'slice.image'
4) bin along one of the axes within that region (here the y axis - the other axes we bin with factor of 1 which corresponds to no binning at all), drop the degenerate axis, and write to a new file 'slice.image'


<source lang="python">
<source lang="python">
ia.open('rotated.image')
ia.open('rotated.image')
ia.rebin(region=mybox,outfile='slice.image',bin=[1,35,1,1],dropdeg=True)
ia.rebin(region=mybox,outfile='slice.image',bin=[1,33,1,1],dropdeg=True)
ia.close('rotated.image')
ia.close('rotated.image')
</source>
</source>


4) load the slice into the viewer - don't forget that the axes need to be swapped in the control panel (in 'Display Axes' e.g. Declination vs Frequency). Unfortunately, the coordinates are now a bit scrambled as they mix RA and Dec. Safest is to display the x-axis as pixels via 'Axis Label Properties' -> 'World or Pixel coordinates' -> 'Pixel'
4) load the slice into the viewer. Unfortunately, the coordinates are now a bit scrambled as they mix RA and Dec. Safest is to display the x-axis as relative coordinates via 'Axis Label Properties' -> 'Absolute or Relative' -> 'Relative'
 
<source lang="python">
viewer 'slice.image'
</source>

Latest revision as of 18:16, 9 November 2010

A position-velocity diagram is typically along slice along a direction in the RA/DEC plane with a given position angle. The y-axis consist of the velocities along that slice. How to do this? Currently there's no task but one can use the tools in CASA.

Here's a recipe:

1) find out the axes of your cube.

imhead('original.image')

and in the example here, the logger spits out:

Axis Coord Type      Name             Proj Shape Tile   Coord value at pixel      Coord incr Units
-------------------------------------------------------------------------------------------------- 
	0    0     Direction Right Ascension   SIN   300   60  09:47:57.382   150.00   -4.000000e-01 arcsec
	1    0     Direction Declination       SIN   300   60 +13.16.40.660   150.00    4.000000e-01 arcsec
	2    1     Stokes    Stokes                    1    1             I
	3    2     Spectral  Frequency                54    9   3.63945e+10     0.00 -1.25014600e+05 Hz
	                     Velocity                              -17.7893     0.00    1.029965e+00 km/s

Note that axis 0 is RA, axis 1 ins Dec, Stokes is axis 2, and frequency/velocity is the last axis 3.

2) rotate the input cube by a given position angle (to align the slice in the vertical or horizontal direction):

#open the image
ia.open('original.image')

#rotate by the position angle and write out a new, rotated cube
ia.rotate(pa='40deg',outfile='rotated.image')

# close the image
ia.close('original.image')

You now have a rotated cube 'rotated.image'

Load and inspect with the viewer

viewer 'rotated.image'


3) now create a box region by finding the bottom left corner (blc) and top right corner (trc) pixel coordinates. Attach this to a variable, here 'mybox'. Note that the last two coords are Stokes (we don't touch it, so running from 0 to 0) and the velocity which we select to be channels 3 through 50:

mybox=rg.box(blc=[99,147,0,3],trc=[205,179,0,50])

4) bin along one of the axes within that region (here the y axis - the other axes we bin with factor of 1 which corresponds to no binning at all), drop the degenerate axis, and write to a new file 'slice.image'

ia.open('rotated.image')
ia.rebin(region=mybox,outfile='slice.image',bin=[1,33,1,1],dropdeg=True)
ia.close('rotated.image')

4) load the slice into the viewer. Unfortunately, the coordinates are now a bit scrambled as they mix RA and Dec. Safest is to display the x-axis as relative coordinates via 'Axis Label Properties' -> 'Absolute or Relative' -> 'Relative'

viewer 'slice.image'