{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Scripting `bettermoments` \n", "\n", "In this Notebook, we will step through how to integrate the moment map making process (in this case, a zeroth moment map, or integrated intensity map), into your workflow. This should elucidate the steps that are taken automatically when using the [command line interface](https://bettermoments.readthedocs.io/en/latest/user/command_line.html).\n", "\n", "### Standard Imports" ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "import bettermoments as bm" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Load Up the Data\n", "\n", "Here the `load_cube` function will return a 3D array for the data and a 1D array of the velocity axis (this should automatically convert any frequency axis to a velocity axis). Note that as we are collapsing along the velocity axis, we have no need for spatial axes, so we do not bother creating them." ] }, { "cell_type": "code", "execution_count": 2, "metadata": {}, "outputs": [], "source": [ "path = '../../../gofish/docs/user/TWHya_CS_32.fits'\n", "data, velax = bm.load_cube(path)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Spectrally Smooth the Data\n", "\n", "If you have relatively noisy data, a low level of smoothing along the spectral axis can be useful. `bettermoments` allows for two different methods: a convolution with a simple top-hat function, or the use of a [Savitzky-Golay filter](https://en.wikipedia.org/wiki/Savitzky-Golay_filter). For a top-hat convolution, you need only specify `smooth`, which describes the kernel size in the number of channels. For a Savitzky-Golay filter, you must also provide `polyorder` which describes the polynomial order which is used for the fitting. Note that `smooth` must always be larger than `polyorder`.\n", "\n", "It is important to remember that while a small level of smoothing can help with certain aspects of moment map creation, it also distorts the line profile (for example broadening the line in the case of a simple top-hat convolution). Such systematic effects must be considered when analysing the resulting moment maps.\n", "\n", "Here we just consider a smoothing with a top-hat kernel that is 3 channels wide." ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "smoothed_data = bm.smooth_data(data=data, smooth=3, polyorder=0)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Estimate the Noise of the Data\n", "\n", "We require an estimate of the noise of the data for two reasons:\n", "\n", "1. For the estimation of the uncertainties of the moment maps.\n", "2. For applying anything threshold clipping.\n", "\n", "To make this estimate, we assume that the noise in the image is constant both spatially (such that the primary beam correction is minimal) and spectrally. To avoid line emission, we consider the RMS of the line-free channels, defined as the first `N` and last `N` channels in the data cube." ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [], "source": [ "rms = bm.estimate_RMS(data=data, N=5)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Note that the noise estimated this way will differ whether you use the `smoothed_data` or the original `data` array. When using the command line interface for `bettermoments`, the RMS will be estimated on the _smoothed_ data." ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "RMS = 2.8 mJy/beam (original)\n", "RMS = 2.2 mJy/beam (smoothed)\n" ] } ], "source": [ "rms_smoothed = bm.estimate_RMS(data=smoothed_data, N=5)\n", "\n", "print('RMS = {:.1f} mJy/beam (original)'.format(rms * 1e3))\n", "print('RMS = {:.1f} mJy/beam (smoothed)'.format(rms_smoothed * 1e3))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### User-Defined Mask\n", "\n", "Sometimes you will want to mask particular regions within your PPV cube in order to disentangle various components, for example if you have multiple hyperfine components that you want to distinguish. Often the easiest way to do this is to define a mask elsewhere and apply it to the data you are collapsing (see for example the [keplerian_mask.py](https://github.com/richteague/keplerian_mask) routine to generate a Keplerian mask).\n", "\n", "Through the `get_user_mask` function, you can load a mask (a 3D array of 1s and 0s) saved as a FITS file, and apply that to the data. If no `user_mask_path` is provided, then this simply returns an array with the same shape as `data` filled with 1s.\n", "\n", "Note that the user-defined mask _must_ share the same pixel and channel resolution, and be the same shape as the data. No aligning or reshaping is done internally with `bettermoments`." ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [], "source": [ "user_mask = bm.get_user_mask(data=data, user_mask_path=None)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Threshold Mask\n", "\n", "A threshold mask, or a 'sigma-clip', is one of the most common approaches to masking used in moment map creation. The `get_threshold_mask` provides several features which will help you optimize your threshold masking.\n", "\n", "The `clip` argument takes a tuple of values, `clip=(-3.0, 3.0)` describing the minimum and maximum SNR of the pixels that will be removed (this is very similar to the `excludepix` argument in [CASA's immoments task](https://casa.nrao.edu/casadocs/casa-6.1.0/global-task-list/task_immoments/about), but with values given in units of sigma, the noise, rather than flux units). `clip` also accepts just a single value, and will convert that to a symmetric clip as above, for example `clip=(-2.0, 2.0)` and `clip=2.0` are equivalent. The option to provide a tuple allows the options to have asymmetric clip ranges, for example, `clip=(-np.inf, 3.0)`, to remove all pixels below 3 sigma, including high significance but negative pixels.\n", "\n", "It has been found that threshold masks can lead to large artifacts in the resulting moment map if there are large intensity gradients in low SNR regions of the PPV cube. To combate this, users have the option to first smooth the data (only temporarily to generate the threshold mask) which will allow for more conservative contours in the threshold mask. This can be achived by providing the FWHM of the Gaussian kernel used for this spatial smooth as `smooth_threshold_mask` in number of pixels. Note that because the data is smoothed, the effective RMS will drop and so the RMS is re-estimated interally on the smoothed image.\n", "\n", "Here we mask all pixels with a SNR less than 2 sigma, i.e., $|I \\, / \\, \\sigma| < 2$." ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [], "source": [ "threshold_mask = bm.get_threshold_mask(data=data,\n", " clip=2.0,\n", " smooth_threshold_mask=0.0)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Channel Mask\n", "\n", "For many PPV cubes, the line emission of interest only spans a small range of velocity axis. This region can be easily selected using the `firstchannel` and `lastchannel` arguments in `get_channel_mask`. Note that the `lastchannel` argument also accepts negative values, following the standard Python indexing convention, i.e., `lastchannel=-1` results in the final channel being the last.\n", "\n", "`get_channel_mask` also accepts a `user_mask` argument, which is an array the same size as the velocity axis of the data, specifying which channels to include. This may be useful if you want to integrate over several hyperfine components while excluding the line-free regions between them." ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [], "source": [ "channel_mask = bm.get_channel_mask(data=data,\n", " firstchannel=0,\n", " lastchannel=-1)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Mask Combination\n", "\n", "All the masks can be easily combined, either with `AND` or `OR`, with the `get_combined_mask` function. This can then be applied to the data used for the moment map creation through a simple multiplication. Note for all collapse functions, pixels with a 0 value will be ignored." ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [], "source": [ "mask = bm.get_combined_mask(user_mask=user_mask,\n", " threshold_mask=threshold_mask,\n", " channel_mask=channel_mask,\n", " combine='and')\n", "masked_data = smoothed_data * mask" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Collapse the Data\n", "\n", "Now that we have a smoothed and masked dataset, we can collapse it along the velocity axis through several different methods. (https://bettermoments.readthedocs.io/en/latest/user/collapse_cube.html#). In general, most functions require the velocity axis, `velax`, the masked data data, `data`, and the RMS of the data, `rms`. The available functions can be checked through the `available_collapse_methods` function such that the desired function is `collapse_{methodname}`." ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Available methods are:\n", "\n", "\t zeroth (integrated intensity)\n", "\t first (intensity weighted average velocity)\n", "\t second (intensity weighted velocity dispersion)\n", "\t eighth (peak intensity)\n", "\t ninth (velocity channel of peak intensity)\n", "\t maximum (both collapse_eighth and collapse_ninth)\n", "\t quadratic (quadratic fit to peak intensity)\n", "\t width (effective width for a Gaussian profile)\n", "\t gaussian (gaussian fit)\n", "\t gaussthick (gaussian with optically thick core fit)\n", "\t gausshermite (gaussian-hermite expansion fit)\n", "\n", "Call the function with `collapse_{method_name}`.\n" ] } ], "source": [ "bm.available_collapse_methods()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Each function will return `moments`, an `(N, Y, X)` shaped array, where `(Y, X)` is the shape of a single channel of the data and `N` is twice the number of statistics (with the uncertainty of each value interleaved). To see which parameters are returned for each `collapse_method`, we can use the `collapse_method_products` function. For the `'zeroth'` method:" ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [ { "data": { "text/plain": [ "'M0, dM0'" ] }, "execution_count": 11, "metadata": {}, "output_type": "execute_result" } ], "source": [ "bm.collapse_method_products('zeroth')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "So we have the zeroth moment, `M0`, and it's associated uncertainty `dM0`.\n", "\n", "Here we will collapse the cube to a zeroth (integrated intensity) map." ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [], "source": [ "moments = bm.collapse_zeroth(velax=velax, data=masked_data, rms=rms)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Save the Data to FITS\n", "\n", "It's possible to work with the data directly, however it's often useful to save these for later. The `save_to_FITS` function will split up the `moments` array and save each one as a new FITS file, replacing the `.fits` exention with `_{moment_name}.fits` for easy identification. The header will be copied from the original file." ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [], "source": [ "bm.save_to_FITS(moments=moments, method='zeroth', path=path)" ] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.7.4" } }, "nbformat": 4, "nbformat_minor": 2 }