Helper Functions

Functions for data I/O, masking and pre-processing that are useful when building custom workflows outside of the command-line interface.

Data I/O

bettermoments.io.load_cube(path, stokes=0)

Load a FITS data cube and return the data and velocity axis.

Parameters:

  • path (str) – Path to the FITS cube.

  • stokes (Optional[int]) – Stokes index to select if the cube has a Stokes axis. Defaults to 0.

Returns:

The data cube with non-finite values

replaced by zero, and the velocity axis in [m/s].

Return type:

data (ndarray), velax (ndarray)

bettermoments.io.save_to_FITS(moments, method, path, outname=None, overwrite=True)

Save the returned fits from collapse_{method_name} as FITS cubes. The filenames will replace the .fits extension with _{param}.fits.

Parameters:

  • moments (array) – Array of moment values from one of the collapse methods.

  • method (str) – Name of the collapse method used, e.g., 'zeroth' if collapse_zeroth was used.

  • path (str) – Path of the original data cube to grab header information.

  • outname (str) – Filename prefix for the saved images. Defaults to the path of the provided FITS file.

  • overwrite (Optional[bool]) – Whether to overwrite files.

Pre-processing & Masking

bettermoments.collapse_cube.estimate_RMS(data, N=5)

Estimate the RMS noise from the first and last N channels of the data cube, using the central 50% of the spatial extent.

Parameters:

  • data (ndarray) – 3D data cube with the spectral axis first.

  • N (Optional[int]) – Number of channels at each end to use for the noise estimate. Defaults to 5.

Returns:

Estimated RMS noise level.

Return type:

float

bettermoments.collapse_cube.smooth_data(data, smooth=0, polyorder=0)

Smooth the input data with a kernel of a width smooth. If polyorder is provided, will smooth with a Savitzky-Golay filter, while if polyorder=0, the default, then only a top-hat kernel will be used. From experimentation, smooth=5 with polyorder=3 provides a good result for noisy, but spectrally resolved data.

Warning

When smoothing low resolution data, this can substantially alter the line profile, so measurements must be taken with caution.

Parameters:

  • data (array) – Data to smooth.

  • smooth (optional[int]) – The width of the kernel for smooth in number of channels.

  • polyorder (optional[int]) – Polynomial order for the Savitzky-Golay filter. This must be smaller than smooth. If not provided, the smoothing will only be a top-hat filter.

Returns:

A smoothed copy of data.

Return type:

smoothed_data (array)

bettermoments.collapse_cube.get_channel_mask(data, firstchannel=0, lastchannel=-1, user_mask=None)

Returns the channel mask (a mask for the zeroth axis) based on a first and last channel. A chan_mask can also be provided for more complex masks, however be warned that the firstchannel and lastchannel will always take precedence over chan_mask.

Parameters:

  • data (array) – The data array to use for masking.

  • firstchannel (optional[int]) – The first channel to include. Defaults to the first channel.

  • lastchannel (optional[int]) – The last channel to include. Defaults to the last channel. This can be both a positive value, or a negative value following the normal indexing conventions, i.e. -1 describes the last channel.

  • user_mask (optional[array]) – A 1D array with size data.shape[0] detailing which channels to include in the moment map creation.

Returns:

A mask array the same shape as data.

Return type:

channel_mask (array)

bettermoments.collapse_cube.get_user_mask(data, user_mask_path=None)

Returns a mask based on a user-provided file. All positive values are included in the mask.

Parameters:

  • data (array) – The data array to mask.

  • user_mask_path (optional[str]) – Path to the FITS cube containing the user-defined mask.

Returns:

A mask array the same shape as data.

Return type:

user_mask (array)

bettermoments.collapse_cube.get_threshold_mask(data, clip=None, rms=None, smooth_threshold_mask=0, noise_channels=5)

Returns a mask based on a sigma-clip to the input data. The most standard approach would be to use clip=3 to mask out all pixels with intensities |I| \leq 3\sigma. If you wanted to specify an asymmetric criteria then you can provide a tuple, clip=(-2, 3) which would mask out all pixels where -2\sigma \leq I \leq 3\sigma.

[Some discussion on the smooth_threshold_mask coming…]

Parameters:

  • data (array) – The data array to mask.

  • clip (optional[float/tuple]) – The sigma clip to apply. If a single value is provided, this is taken to be a symmetric mask. If a tuple if provided, this is taking as a minimum and maximum clip value.

  • rms (optional[float]) – The RMS level to use for defining the noise. If not specified, will calculate it based on the standard deviation of the data.

  • smooth_threshold_mask (optional[float]) – Convolution kernel FWHM in pixels.

  • noise_channels (optional[int]) – Number of channels at the start and end of the velocity axis to use for estimating the noise.

Returns:

A mask array the same shape as data.

Return type:

threshold_mask (array)

bettermoments.collapse_cube.get_combined_mask(user_mask, threshold_mask, channel_mask, combine='and')

Return the combined user, threshold and channel masks, user_mask, threshold_mask and velo_mask, respectively. The user and threshold masks can be combined either through AND or OR, which is controlled through the combine argument. This defaults to AND, such that all mask requirements are met.

Parameters:

  • user_mask (array) – User-defined mask from get_user_mask.

  • threshold_mask (array) – Threshold mask from get_threshold_mask.

  • channel_mask (array) – Channel mask from get_channel_mask.

Returns:

A combined mask.

Return type:

combined_mask (array)

Spectral Noise Correlation

Tools for handling spectrally correlated noise (e.g.from Hanning smoothing in the imaging pipeline). These underpin the --acf command-line flag and the acf= keyword on the collapse_* functions, and can also be used directly when building custom workflows.

bettermoments.collapse_cube.estimate_spectral_acf(data, N=5, max_lag=None, rms=None, threshold=2.0)

Estimate the normalised spectral autocorrelation function (ACF) of the noise from off-source spatial pixels, using the full spectral axis.

Off-source pixels are identified as those whose peak absolute intensity across all channels is below threshold * rms. Each such full-length spectrum is mean-subtracted and used to form the biased ACF estimator,

\hat{\rho}(\tau) = \frac{\sum_{i} n_i n_{i+\tau}}{\sum_{i} n_i^2},

averaged across all selected pixels. The returned ACF is truncated at the first lag at which it falls within the white-noise band \pm 2/\sqrt{N_{\rm chan}}, where N_{\rm chan} is the number of channels per spectrum. acf[0] is always 1.

Note that with threshold = 2 a fraction of true noise-only sightlines will be rejected (the maximum of N_chan Gaussian draws routinely exceeds 2 sigma), so the selection is conservative: it favours including only clearly off-source pixels at the cost of sample size, which is fine for a first-order ACF estimate.

Parameters:

  • data (ndarray) – 3D data cube with the spectral axis first.

  • N (Optional[int]) – Number of channels at each end of the cube used to estimate rms when one is not supplied. Defaults to 5.

  • max_lag (Optional[int]) – Hard cap on the maximum lag returned. If None, defaults to N_chan - 1.

  • rms (Optional[float]) – Per-channel RMS used to define the off-source threshold. If None, estimated via estimate_RMS().

  • threshold (Optional[float]) – Pixels with peak |intensity| < threshold * rms across all channels are treated as off-source. Defaults to 2.0.

Returns:

1D array of normalised ACF values, acf[0] = 1, length

between 1 and max_lag + 1.

Return type:

ndarray

bettermoments.collapse_cube.build_spectral_covariance(rms, acf, nchan, eps=1e-06)

Build the Toeplitz spectral noise covariance matrix,

\mathbf{C}_{ij} = \sigma^2 \rho(|i - j|),

where \rho is the normalised autocorrelation function returned by estimate_spectral_acf() and \sigma is the per-channel RMS. Lags beyond len(acf) - 1 are taken to be zero.

Truncated empirical ACFs need not yield a positive-definite Toeplitz matrix (the implied spectral density can go negative at high frequency). To guarantee a valid covariance suitable for Cholesky decomposition, the eigenvalues are clipped at a small fraction eps of the spectral radius before reconstruction.

Parameters:

  • rms (float) – Per-channel RMS noise.

  • acf (ndarray) – 1D normalised ACF with acf[0] = 1.

  • nchan (int) – Number of channels (size of the returned matrix).

  • eps (Optional[float]) – Eigenvalue floor as a fraction of the largest eigenvalue. Defaults to 1e-6.

Returns:

(nchan, nchan) symmetric positive-definite covariance.

Return type:

ndarray

Utilities

bettermoments.methods.check_finite_errors(moments)

Check if the errors are finite, and if not, replace them with nearest-neighbor interpolated values.

Parameters:

moments (tuple) – Tuple of arrays in (value, error, value, error, ...) order, as returned by one of the collapse_* functions.

Returns:

Same structure as moments with any

non-finite error values replaced by interpolated finite errors.

Return type:

verified_moments (tuple)