How to rotate and slice a cube for pV diagrams
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'