Smoothing¶
There are two types of smoothing routine available in spectral_cube
:
spatial and spectral.
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.
Common Beam selection¶
You may want to convolve your cube to the smallest beam that is still larger
than all contained beams. To do this, you can use the
common_beam
tool. For example:
common_beam = cube.beams.common_beam()
new_cube = cube.convolve_to(common_beam)
Sometimes, you’ll encounter the error “Could not find common beam to deconvolve all beams.” This is a real issue, as the algorithms we have in hand so far do not always converge on a common containing beam. There are two ways to get the algorithm to converge to a valid common beam:
1. Changing the tolerance. - You can try to change the tolerance used in the
getMinVolEllipse
code by
passing tolerance=1e-5
to the common beam function:
cube.beams.common_beam(tolerance=1e-5)
Convergence may be met by either increasing or decreasing the tolerance; it depends on having the algorithm not step within the minimum enclosing ellipse, leading to the error. Note that decreasing the tolerance by an order of magnitude will require an order of magnitude more iterations for the algorithm to converge and will take longer to run.
2. Changing epsilon - A second parameter epsilon
controls the fraction
to overestimate the beam size, ensuring that solutions that are marginally
smaller than the common beam will not be found by the algorithm:
cube.beams.common_beam(epsilon=1e-3)
The default value of epsilon=1e-3
will sample points 0.1% larger than the
edge of each beam in the set. Increasing epsilon
ensures that a valid common
beam can be found, avoiding the tolerance issue, but will result in
overestimating the common beam area. For most radio data sets, where the beam
is oversampled by \(\sim 3 \mbox{-5}\) pixels, moderate increases in
epsilon
will increase the common beam area far less than a pixel area, making
the overestimation negligible.
We recommend testing different values of tolerance to find convergence, and if the error persists, to then slowly increase epsilon until a valid common beam is found. More information can be found in the radio-beam documentation.
Alternative approach to spatial smoothing¶
There is an alternative way to spatially smooth a data cube, which is using the
spatial_smooth()
method. This is an example
of how you can do this:
from spectral_cube import SpectralCube
from astropy.io import fits
from astropy.convolution import Gaussian2DKernel
cube = SpectralCube.read('/some_path/some_file.fits')
kernel = Gaussian2DKernel(x_stddev=1)
new_cube = cube.spatial_smooth(kernel)
new_cube.write('/some_path/some_other_file.fits')
x_stddev
specifies the width of the Gaussian kernel.
Any astropy convolution
is acceptable.
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.value)
new_cube = cube.spectral_smooth(kernel)
new_cube.write('newfile.fits')
gaussian_width
is in pixel units but is defined as a unitless Quantity
.
By using gaussian_width.value
, we convert the pixel width into a float.