Masked moments of radio data cubes, using a dilated masking approach. The basic idea is that the mask is produced from the cube itself, by taking a high significance contour (e.g. 4-sigma) and expanding it to a surrounding lower-significance contour (e.g. 2-sigma). User has the option to smooth the data (spatially and/or spectrally) before defining these contours, and to pad the mask by a certain number of pixels in all directions. An estimate of the noise at each location in the cube is needed, which can be provided from a different input cube, or can be calculated by the script from the data, optionally using a gain image.
Based on idl_mommaps: https://github.com/tonywong94/idl_mommaps
Required packages: spectral_cube, radio_beam, astropy, scipy
How to use:
from maskmoment import maskmoment
maskmoment(img_fits='file.fits', *other parameters*)
img_fits : FITS file name, required
The image cube, this should be in units of K, Jy/beam, or equivalent.
gain_fits : FITS file name, optional
The gain cube, e.g. pb cube from CASA. This should have a value between
0 and 1, with 0 near the edges and 1 near the center of the image, and
have the same dimensions as (or broadcastable to) the image cube.
NOTE: The gain cube is ignored if a noise cube (rms_fits) is given.
rms_fits : FITS file name, optional
The noise cube, providing an estimate of the rms noise at each pixel.
This should have the same dimensions and units as the image cube.
If it is a 2D array it will be replicated along the velocity axis.
NOTE: If rms_fits is not given, a noise cube is generated from the
image cube, after removing any gain variation using the gain cube.
rms : float, optional
Global estimate of the rms noise, to be used instead of trying
to estimate it from the data. It should have the units of the input cube.
Default is to use rms_fits, or to calculate the rms from the data.
mask_fits : FITS file name, optional
External mask cube to use. This cube should have 1's for valid pixels
and 0's for excluded pixels. If this is provided then the mask generation
is skipped and the program goes straight to calculating the moments.
outdir : string, optional
Directory to write the output files.
Default: Write to the directory where img_fits resides.
outname : string, optional
Basename for output files. For instance, outname='foo' produces files
'foo.mom0.fits.gz', etc.
Default: Based on root name of img_fits.
snr_hi : float, optional
The high significance sigma threshold from which to begin mask dilation.
Use a very low number (such as -10) to produce unmasked moment maps.
Default: 4
snr_lo : float, optional
The low significance sigma threshold at which to end mask dilation.
Default: 2
snr_hi_minch : int, optional
High significance mask is required to span at least this many channels
at all pixels.
Default: 1
snr_lo_minch : int, optional
Low significance mask is required to span at least this many channels
at all pixels.
Default: 1
min_tot_ch : int, optional
Each contiguous mask region must span at least this many channels (but
it's not necessary that every pixel in the region span this many channels).
Default: 2
minbeam : float, optional
Minimum velocity-integrated area of a mask region in units of the beam size.
Default: 1
nguard : tuple of two ints, optional
Expand the final mask by nguard[0] pixels in the sky directions and
nguard[1] channels in velocity. Currently these values must be equal
if both are non-zero.
If nguard[0] = 0 then no expansion is done in sky coordinates.
If nguard[1] = 0 then no expansion is done in velocity.
Default: [0,0]
edgech : int, optional
Number of channels at each end of vel axis to use for rms estimation.
Default: 5
fwhm : float or :class:`~astropy.units.Quantity`, optional
Spatial resolution to smooth to before generating the dilated mask.
If value is not an astropy quantity, assumed to be given in arcsec.
Default: No spatial smoothing is applied.
vsm : float or :class:`~astropy.units.Quantity`, optional
Full width of the spectral smoothing kernel (or FWHM for gaussian).
If given as astropy quantity, should be given in velocity units.
If not given as astropy quantity, interpreted as number of channels.
Default: No spectral smoothing is applied.
vsm_type : string, optional
What type of spectral smoothing to employ. Currently three options:
(1) 'boxcar' - 1D boxcar smoothing, vsm rounded to integer # of chans.
(2) 'gauss' - 1D gaussian smoothing, vsm is the convolving gaussian FWHM.
(3) 'gaussfinal' - 1D gaussian smoothing, vsm is the gaussian FWHM
after convolution, assuming FWHM before convolution is 1 channel.
Default: 'gauss'
perpixel : boolean, optional
Whether to calculate the rms per XY pixel instead of over whole image.
Set to True if you know there is a sensitivity variation across the image
but you don't have a gain cube - requires rms_fits and gain_fits unset.
Default: False
output_peak : boolean, optional
Output the peak brightness, velocity at peak brightness, and effective
line width (mom0/peak) in addition to the moment maps. The peak brightness
has no mask applied, the other two are derived from the masked cube. The
effective line width is normalized to match mom-2 for a pure Gaussian profile.
Default: False
output_snr_cube : boolean, optional
Output the cube in SNR units in addition to the moment maps.
Default: False
output_snr_peak : boolean, optional
Output the peak SNR image in addition to the moment maps.
Default: False
output_2d_mask : boolean, optional
Output the projected 2-D mask as well as the newly generated mask.
The projected mask at a given pixel is valid for all channels as
long as the parent mask is valid for any channel.
Default: False
to_kelvin : boolean, optional
Output the moment maps in K units if the cube is in Jy/beam units.
Default: True
vmin : float or :class:`~astropy.units.Quantity`, optional
Minimum channel number or velocity to use for all calculations. Note
that channels are discarded before any rms estimation using edgech.
If given as astropy quantity, should be provided in velocity units.
If NOT given as astropy quantity, interpreted as channel number.
Default: Do not impose a velocity cut.
vmax : float or :class:`~astropy.units.Quantity`, optional
Maximum channel number or velocity to use for all calculations. Note
that channels are discarded before any rms estimation using edgech.
If given as astropy quantity, should be provided in velocity units.
If NOT given as astropy quantity, interpreted as channel number.
Default: Do not impose a velocity cut.
Base code was developed by Tony Wong in 2019-2020. Hailin Wang assisted with scripting, testing, and debugging the code during development.