Manipulating cubes and extracting subcubes¶
Modifying the spectral axis¶
As mentioned in Accessing data, it is straightforward to find the
coordinates along the spectral axis using the
spectral_axis
attribute:
>>> cube.spectral_axis
[ -2.97198762e+03 -2.63992044e+03 -2.30785327e+03 -1.97578610e+03
-1.64371893e+03 -1.31165176e+03 -9.79584583e+02 -6.47517411e+02
...
3.15629983e+04 3.18950655e+04 3.22271326e+04 3.25591998e+04
3.28912670e+04 3.32233342e+04] m / s
The default units of a spectral axis are determined from the FITS header or
WCS object used to initialize the cube, but it is also possible to change the
spectral axis unit using with_spectral_unit()
:
>>> from astropy import units as u
>>> cube2 = cube.with_spectral_unit(u.km / u.s)
>>> cube2.spectral_axis
[ -2.97198762e+00 -2.63992044e+00 -2.30785327e+00 -1.97578610e+00
-1.64371893e+00 -1.31165176e+00 -9.79584583e-01 -6.47517411e-01
...
3.02347296e+01 3.05667968e+01 3.08988639e+01 3.12309311e+01
3.15629983e+01 3.18950655e+01 3.22271326e+01 3.25591998e+01
3.28912670e+01 3.32233342e+01] km / s
It is also possible to change from velocity to frequency for example, but this requires specifying the rest frequency or wavelength as well as a convention for the doppler shift calculation:
>>> cube3 = cube.with_spectral_unit(u.GHz, velocity_convention='radio',
... rest_value=200 * u.GHz)
[ 220.40086492 220.40062079 220.40037667 220.40013254 220.39988841
220.39964429 220.39940016 220.39915604 220.39891191 220.39866778
...
220.37645231 220.37620818 220.37596406 220.37571993 220.3754758
220.37523168 220.37498755 220.37474342 220.3744993 220.37425517] GHz
The new cubes will then preserve the new spectral units when computing moments for example (see Moment maps and statistics).
Extracting a spectral slab¶
Given a spectral cube, it is easy to extract a sub-cube covering only a subset
of the original range in the spectral axis. To do this, you can use the
spectral_slab()
method. This
method takes lower and upper bounds for the spectral axis, as well as an
optional rest frequency, and returns a new
SpectralCube
instance. The bounds can
be specified as a frequency, wavelength, or a velocity but the units have to
match the type of the spectral units in the cube (if they do not match, first
use with_spectral_unit()
to ensure that they
are in the same units). The bounds should be given as Astropy
Quantities
as follows:
>>> from astropy import units as u
>>> subcube = cube.spectral_slab(-50 * u.km / u.s, +50 * u.km / u.s)
The resulting cube subcube
(which is also a
SpectralCube
instance) then contains all channels
that overlap with the range -50 to 50 km/s relative to the rest frequency
assumed by the world coordinates, or the rest frequency specified by a prior
call to with_spectral_unit()
.
Extracting a sub-cube by indexing¶
It is also easy to extract a sub-cube from pixel coordinates using standard Numpy slicing notation:
>>> sub_cube = cube[:100, 10:50, 10:50]
This returns a new SpectralCube
object
with updated WCS information.
Extracting a subcube from a DS9/CRTF region¶
You can use DS9/CRTF regions to extract subcubes. The minimal enclosing subcube will be extracted with a two-dimensional mask corresponding to the DS9/CRTF region. Regions is required for region parsing. CRTF regions may also contain spectral cutout information.
This example shows extraction of a subcube from a ds9 region file file.reg
.
read
parses the ds9 file and converts it to a list of
Region
objects:
>>> import regions
>>> region_list = regions.Regions.read('file.reg')
>>> sub_cube = cube.subcube_from_regions(region_list)
This one shows extraction of a subcube from a CRTF region file file.crtf
,
parsed using read_crtf
:
>>> import regions
>>> region_list = regions.read_crtf('file.reg')
>>> sub_cube = cube.subcube_from_regions(region_list)
If you want to loop over individual regions with a single region file, you need to convert the individual regions to lists of that region:
>>> region_list = regions.Regions.read('file.reg')
>>> for region in region_list:
>>> sub_cube = cube.subcube_from_regions([region])
You can also directly use a ds9 region string. This example extracts a 0.1 degree circle around the Galactic Center:
>>> region_str = "galactic; circle(0, 0, 0.1)"
>>> sub_cube = cube.subcube_from_ds9region(region_str)
Similarly, you can also use a CRTF region string:
>>> region_str = "circle[[0deg, 0deg], 0.1deg], coord=galactic, range=[150km/s, 300km/s]"
>>> sub_cube = cube.subcube_from_crtfregion(region_str)
CRTF regions that specify a subset in the spectral dimension can be used to
produce full 3D cutouts. The meta
attribute of a regions.Region
object
contains the spectral information for that region in the three special keywords
range
, restfreq
, and veltype
:
>>> import regions
>>> from astropy import units as u
>>> regpix = regions.RectanglePixelRegion(regions.PixCoord(0.5, 1), width=4, height=2)
>>> regpix.meta['range'] = [150 * u.km/u.s, 300 * u.km/u.s] # spectral range
>>> regpix.meta['restfreq'] = [100 * u.GHz] # rest frequency
>>> regpix.meta['veltype'] = 'OPTICAL' # velocity convention
>>> subcube = cube.subcube_from_regions([regpix])
If range
is specified, but the other two keywords are not, the code will
likely crash.
Extract the minimal valid subcube¶
If you have a mask that masks out some of the cube edges, such that the resulting sub-cube might be smaller in memory, it can be useful to extract the minimal enclosing sub-cube:
>>> sub_cube = cube.minimal_subcube()
You can also shrink any cube by this mechanism:
>>> sub_cube = cube.with_mask(smaller_region).minimal_subcube()
If you have many cubes that are the same shape, and you want to cut them out in the same way (e.g., for CASA image, model, residual, and other cubes), you can get the slice to make the cutout and reuse it. It can also be helpful to cut only in the spatial dimensions:
>>> subcube_slice = cube.subcube_slices_from_mask(cube.mask, spatial_only=True)
>>> mod_subcube = modcube[subcube_slice]
>>> resid_subcube = residcube[subcube_slice]
>>> subcube = cube[subcube_slice]
Extract a spatial and spectral subcube¶
There is a generic subcube function that allows slices in the spatial and spectral axes simultaneously, as long as the spatial axes are aligned with the pixel axes. An arbitrary example looks like this:
>>> sub_cube = cube.subcube(xlo=5*u.deg, xhi=6*u.deg,
ylo=2*u.deg, yhi=2.1*u.deg,
zlo=50*u.GHz, zhi=51*u.GHz)