Smoothing¶
There are two types of smoothing routine available in spectral_cube
:
spectral and spatial.
Spatial Smoothing¶
The convolve_to
method will convolve each plane
of the cube to a common resolution, assuming the cube’s resolution is known
in advanced and stored in the cube’s beam
or beams
attribute.
A simple example:
import radio_beam
from spectral_cube import SpectralCube
from astropy import units as u
cube = SpectralCube.read('file.fits')
beam = radio_beam.Beam(major=1*u.arcsec, minor=1*u.arcsec, pa=0*u.deg)
new_cube = cube.convolve_to(beam)
Note that the convolve_to()
method will work
for both VaryingResolutionSpectralCube
instances and
single-resolution SpectralCube
instances, but for a
VaryingResolutionSpectralCube
, the convolution kernel
will be different for each slice.
Spectral Smoothing¶
Only SpectralCube
instances with a consistent beam can
be spectrally smoothed, so if you have a
VaryingResolutionSpectralCube
, you must convolve each
slice in it to a common resolution before spectrally smoothing.
spectral_smooth()
will apply a convolution
kernel to each spectrum in turn. As of July 2016, a parallelized version is
partly written but incomplete.
Example:
import radio_beam
from spectral_cube import SpectralCube
from astropy import units as u
from astropy.convolution import Gaussian1DKernel
cube = SpectralCube.read('file.fits')
kernel = Gaussian1DKernel(2.5)
new_cube = cube.spectral_smooth(kernel)
This can be useful if you want to interpolate onto a coarser grid but maintain
Nyquist sampling. You can then use the
spectral_interpolate
method to regrid your
smoothed spectrum onto a new grid.
Say, for example, you have a cube with 0.5 km/s resolution, but you want to resample it onto a 2 km/s grid. You might then choose to smooth by a factor of 4, then downsample by the same factor:
# cube.spectral_axis is np.arange(0,10,0.5) for this example
new_axis = np.arange(0,10,2)*u.km/u.s
fwhm_factor = np.sqrt(8*np.log(2))
smcube = cube.spectral_smooth(Gaussian1DKernel(4/fwhm_factor))
interp_Cube = smcube.spectral_interpolate(new_axis,
suppress_smooth_warning=True)
We include the suppress_smooth_warning
override because there is no way for
SpectralCube
to know if you’ve done the appropriate smoothing (i.e., making
sure that your new grid nyquist samples the data) prior to the interpolation
step. If you don’t specify this, it will still work, but you’ll be warned that
you should preserve Nyquist sampling.
If you have a cube with 0.1 km/s resolution (where we assume resolution corresponds to the fwhm of a gaussian), and you want to smooth it to 0.25 km/s resolution, you can smooth the cube with a Gaussian Kernel that has a width of (0.25^2 - 0.1^2)^0.5 = 0.229 km/s. For simplicity, it can be done in the unit of pixel. In our example, each channel is 0.1 km/s wide:
import numpy as np
from astropy import units as u
from spectral_cube import SpectralCube
from astropy.convolution import Gaussian1DKernel
cube = SpectralCube.read('file.fits')
fwhm_factor = np.sqrt(8*np.log(2))
current_resolution = 0.1 * u.km/u.s
target_resolution = 0.25 * u.km/u.s
pixel_scale = 0.1 * u.km/u.s
gaussian_width = ((target_resolution**2 - current_resolution**2)**0.5 /
pixel_scale / fwhm_factor)
kernel = Gaussian1DKernel(gaussian_width)
new_cube = cube.spectral_smooth(kernel)
new_cube.write('newfile.fits')