SpectralCubeReprojectMaskExample_v4p1.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "3898bc3f",
   "metadata": {},
   "source": [
    "# Cube Reprojection and Signal Masking Tutorial, version 4\n",
    "\n",
    "### Author: Alyssa Bulatek (https://github.com/abulatek)\n",
    "\n",
    "#### Last updated: May 9, 2022\n",
    "\n",
    "Adapted from a previous tutorial on the same subject by Adam Ginsburg, Eric Koch, Tom Robitaille, and Erik Rosolowsky. This tutorial shows how to take two spectral cubes observed toward the same part of the sky, but at different frequencies, and put them onto the same grid. Then, we use one spectral cube to mask the other, which is helpful for situations where the velocity structure of one cube is more complicated than the other."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d75448ba",
   "metadata": {},
   "source": [
    "## Index \n",
    "\n",
    " * [Step 0: Some imports and version checks](#Step-0:-Some-imports-and-version-checks)\n",
    " * [Step 1: Download](#Step-1:-Download-the-data)\n",
    " * [Step 2: Open files, collect metadata](#Step-2:-Load-the-cubes)\n",
    " * [Step 3: Convert to velocity](#Step-3:-Convert-cubes-from-frequency-to-velocity)\n",
    " * [Step 4: Spectral Interpolation](#Step-4:-Spectral-Interpolation)\n",
    " * [Step 5: Spatial Smoothing and Reprojection](#Step-5:-Spatial-Smoothing-and-Reprojection)\n",
    " * [Step 6: Signal Masking](#Step-6:-Signal-Masking)\n",
    " * [Step 7: Moment Maps](#Step-7:-Moment-Maps)\n",
    " * [Step 8: Masking Outliers](#Step-8.-Masking-Outliers)\n",
    " \n",
    "***NB:*** *In this example, we do spectral smoothing and interpolation (step 4) before spatial smoothing and interpolation (step 5), but if you have a varying-resolution cube (with a different beam size for each channel), you have to do spatial smoothing first. This is because smoothing two spectral channels together that have different beamsizes will blend unrelated objects on the sky.*"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e180368c",
   "metadata": {},
   "source": [
    "## Step 0: Some imports and version checks"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "303c0625",
   "metadata": {},
   "outputs": [],
   "source": [
   "%pip install -r https://github.com/radio-astro-tools/tutorials/raw/master/requirements.txt\n",
   "%pip install --pre -U astroquery"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b65a9e44",
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "# # Make plots have dark backgrounds; run this if your notebook is in dark mode\n",
    "# plt.style.use('dark_background')\n",
    "\n",
    "# Make sure spectral_cube is up-to-date\n",
    "from spectral_cube import SpectralCube\n",
    "import spectral_cube\n",
    "print(f'spectral_cube version: {spectral_cube.__version__}')\n",
    "print(f'path to spectral_cube installation: {spectral_cube.__path__}')\n",
    "\n",
    "# Make sure astropy and reproject are somewhat up-to-date\n",
    "# Used to require development version of astropy; 4.3.dev1788+ga3263b6 worked for me, but now so does 4.3.1\n",
    "import astropy\n",
    "import reproject\n",
    "print(f'astropy version: {astropy.__version__}')\n",
    "print(f'reproject version: {reproject.__version__}')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1e71967c",
   "metadata": {},
   "source": [
    "## Step 1: Download the data\n",
    "\n",
    "The first step is to pick out and download the cubes you will use for the reprojection. Here, we choose to download two cubes from the [ALMA Archive](https://almascience.nrao.edu/aq/), which is a publicly-available database where data taken by ALMA are stored after they exit their proprietary period. If you already have the cubes you will use for this tutorial locally, you can skip this step.\n",
    "\n",
    "The two cubes we download here both cover the same field at different wavelengths: a water maser in the Galactic Center molecular cloud G0.253+0.016, also known as \"The Brick.\" `filename_1` points to a cube containing the H2CS 4(1,3)-3(1,2) line (139.483699 GHz), and `filename_2` points to a cube containing the SiO v=5-4 line (217.104984 GHz)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bcb5dd48",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from astropy.utils.data import download_file"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6876df20",
   "metadata": {},
   "outputs": [],
   "source": [
    "filename_1 = download_file(\"https://almascience.nrao.edu/dataPortal/member.uid___A001_X1465_X3a33.BrickMaser_sci.spw71.cube.I.manual.image.pbcor.fits\",\n",
    "                           cache=True)\n",
    "filename_1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bdaccd0a",
   "metadata": {},
   "outputs": [],
   "source": [
    "filename_2 = download_file(\"https://almascience.nrao.edu/dataPortal/member.uid___A001_X87d_X141.a_sma1_sci.spw27.cube.I.pbcor.fits\",\n",
    "                          cache=True)\n",
    "filename_2"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "e9cc1664",
   "metadata": {},
   "source": [
    "## Step 2: Load the cubes\n",
    "\n",
    "Now, we will load the cubes into memory using the package `spectral_cube`, which we imported in Step 0. As previously mentioned, the two cubes are at different frequencies (139 GHz and 217 GHz). Once we've loaded the cubes in, we'll use the method `find_lines` to confirm the presence of the line we expect to find in each of the cubes (H2CS in `cube1`, SiO in `cube2`)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "596b65cc",
   "metadata": {},
   "outputs": [],
   "source": [
    "cube1 = SpectralCube.read(filename_1)\n",
    "cube1.allow_huge_operations = True\n",
    "cube1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7d19a001",
   "metadata": {},
   "outputs": [],
   "source": [
    "cube2 = SpectralCube.read(filename_2)\n",
    "cube2.allow_huge_operations = True\n",
    "cube2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ffa01793-a147-4c15-8d2c-421c141cf624",
   "metadata": {},
   "outputs": [],
   "source": [
    "# DEBUG\n",
    "import astroquery\n",
    "print(astroquery.__version__)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "765e0051",
   "metadata": {},
   "outputs": [],
   "source": [
    "cube1.find_lines(chemical_name=' H2CS ')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1fcc119b",
   "metadata": {},
   "outputs": [],
   "source": [
    "cube2.find_lines(chemical_name='SiO')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3a025ea9",
   "metadata": {},
   "source": [
    "We can see four listings for the H2CS 4(1,3)-3(1,2) line (139.483699 GHz) in `cube1`, each from different linelists. We can also see a few different SiO v=5-4 lines (217.104984 GHz), some from the same linelist (CDMS) and some from the other linelists. These printouts are just for verification purposes; we'll use the center frequencies we have written down in Step 1 for each of these transitions instead of the frequencies given in the `find_lines` function."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b6f595d7",
   "metadata": {},
   "source": [
    "### Convert from Jy/beam to kelvin\n",
    "\n",
    "At this point, we'll convert the cubes to be in units of kelvin instead of Jy/beam. This is because we'll be smoothing the cubes spatially in a future step, and the data values will only be preserved if the cubes are in units of kelvin instead of Jy/beam (since the latter units are dependent on the beamsize)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "deb201b6",
   "metadata": {},
   "outputs": [],
   "source": [
    "from astropy import units as u\n",
    "cube1 = cube1.to(u.K)\n",
    "cube2 = cube2.to(u.K)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9f11afbc",
   "metadata": {},
   "source": [
    "## Step 3: Convert cubes from frequency to velocity\n",
    "\n",
    "Now, we will convert these cubes from frequency units to velocity units using the rest frequencies of the transitions of interest. This is the first step to getting the cubes on a commmon velocity grid: we match up their rest frequencies at 0 km/s. We convert these cubes from frequency to velocity using the method `with_spectral_unit`, the documentation for which is available [here](https://spectral-cube.readthedocs.io/en/latest/api/spectral_cube.SpectralCube.html#spectral_cube.SpectralCube.with_spectral_unit) in the `spectral_cube` API."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "053eddbe",
   "metadata": {},
   "outputs": [],
   "source": [
    "cube1vel = cube1.with_spectral_unit(u.km/u.s, velocity_convention='radio', rest_value=139.483699*u.GHz)\n",
    "cube1vel"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "84f9d272",
   "metadata": {},
   "outputs": [],
   "source": [
    "cube2vel = cube2.with_spectral_unit(u.km/u.s, velocity_convention='radio', rest_value=217.104984*u.GHz)\n",
    "cube2vel"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "399cf7a8",
   "metadata": {},
   "source": [
    "From the shape of the cube, we can see the H2CS cube is narrower in velocity, so we'll use that as the target spectral reprojection (otherwise, we would have empty velocity channels on either side of the smaller cube that would be getting reprojected).  However, the SiO cube is the smaller footprint on the sky."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3e01349a",
   "metadata": {},
   "source": [
    "### Do some quicklooks of the peak intensity to see what we expect to see\n",
    "\n",
    "The cubes we are working with are three-dimensional: they have two spatial axes (on the plane of the sky) and one spectral axis. It is hard to visualize them in 3D easily, so we can collapse them into two dimensions in a variety of different ways. Here, we will collapse them along their spectral axis (`axis = 0`; `spectral_cube` usually makes the spectral axis the first axis as a default). In particular, for any given spatial coordinate in our cube, we will take the value of the brightest pixel along the spectral axis (the **peak flux density**). This summary image shows us where the bright emission is (spatially) in the cube, and how bright the emission in a certain region is in relation to the brightness of the emission in another region."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "31e042bf",
   "metadata": {},
   "outputs": [],
   "source": [
    "peak_fluxdensity1 = cube1.max(axis = 0)\n",
    "peak_fluxdensity1.quicklook()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1f850037",
   "metadata": {},
   "outputs": [],
   "source": [
    "peak_fluxdensity2 = cube2.max(axis = 0)\n",
    "peak_fluxdensity2.quicklook()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7809d54e",
   "metadata": {},
   "source": [
    "## Step 4: Spectral Interpolation\n",
    "\n",
    "We can safely assume that any two cubes from ALMA are not perfectly matched on a spatial and spectral grid. One cube will have a higher spatial resolution than the other, and one cube will have a finer spectral resolution than the other. We must match these two cubes in two different steps: matching them on a spatial grid (2D), and matching them along their spectral axes (1D). We can treat these as roughly independent steps, and could do either first. *However, if your cubes have a varying spatial resolution (different beam size at each spectral channel), you will need to perform spatial smoothing first—you don't want to blend unrelated components of your source by smoothing spectrally before all channels are smoothed to the same spatial resolution.* Because each of our cubes have a constant spatial resolution, we will perform spectral smoothing first.\n",
    "\n",
    "We need to match the spectral resolution of both cubes to each other. We'll start by printing out the spectral resolutions of both cubes."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e4e399f9",
   "metadata": {},
   "outputs": [],
   "source": [
    "velocity_res_1 = np.diff(cube1vel.spectral_axis)[0]\n",
    "velocity_res_2 = np.diff(cube2vel.spectral_axis)[0]\n",
    "velocity_res_1, velocity_res_2"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "46c42867",
   "metadata": {},
   "source": [
    "`cube1`, our H2CS cube, has broader channels.  We want to smooth `cube2` to that channel width before doing the spatial reprojection. To do this, we'll smooth the spectral axis of `cube2` using a Gaussian with an appropriate FWHM:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e91fdbf5",
   "metadata": {},
   "outputs": [],
   "source": [
    "fwhm_gaussian = (velocity_res_1**2 - velocity_res_2**2)**0.5\n",
    "fwhm_gaussian"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4c4364c9",
   "metadata": {},
   "source": [
    "Here is the [page on smoothing](https://spectral-cube.readthedocs.io/en/latest/smoothing.html#spectral-smoothing) in the `spectral_cube` documentation. Importantly, here, we're using a Gaussian to smooth the spectral channels, and then we say that the spectral resolution of the resulting cubes is equal to the FWHM of that Gaussian. However, spectral resolution is defined differently by the ALMA correlator: Adam thinks it's half of the channel width smoothed by a Hanning function and then sampled by the channel width, but the answer lies somewhere in the ALMA documentation. *The takeaway: we are not necessarily using the \"correct\" definition of resolution, but we could figure out what that is, if we wanted to.*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8dbfa36c",
   "metadata": {},
   "outputs": [],
   "source": [
    "from astropy.convolution import Gaussian1DKernel\n",
    "fwhm_to_sigma = np.sqrt(8*np.log(2))\n",
    "# We want the kernel in pixel units, so we force to km/s and take the value\n",
    "spectral_smoothing_kernel = Gaussian1DKernel(stddev = fwhm_gaussian.to(u.km/u.s).value / fwhm_to_sigma)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "49e56298",
   "metadata": {},
   "source": [
    "Then, let's cut down to the same spectral coverage. This is not technically necessary to do, but it is more computationally efficient to smooth over a smaller spectral range. We will do this in a general way, determining the smallest coverage and cutting both cubes based on those values, so no matter which cube has a larger spectral coverage, this code should work."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2ef227da",
   "metadata": {},
   "outputs": [],
   "source": [
    "print(cube1vel.spectral_extrema)\n",
    "print(cube2vel.spectral_extrema)\n",
    "\n",
    "lower_limit = np.max([cube1vel.spectral_extrema[0].value, cube2vel.spectral_extrema[0].value])*(cube1vel.spectral_axis.unit)\n",
    "upper_limit = np.min([cube1vel.spectral_extrema[1].value, cube2vel.spectral_extrema[1].value])*(cube1vel.spectral_axis.unit)\n",
    "\n",
    "print(lower_limit)\n",
    "print(upper_limit)\n",
    "\n",
    "cube1vel_cutout = cube1vel.spectral_slab(lower_limit, upper_limit)\n",
    "cube2vel_cutout = cube2vel.spectral_slab(lower_limit, upper_limit)\n",
    "cube1vel_cutout, cube2vel_cutout"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8b9dc868",
   "metadata": {},
   "source": [
    "We then smooth with the kernel. Note that this is doing 420x420 = 176400 smoothing operations on a length-221 spectrum: **it will take several minutes to run.**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "84d9cd8f",
   "metadata": {},
   "outputs": [],
   "source": [
    "cube2vel_smooth = cube2vel_cutout.spectral_smooth(spectral_smoothing_kernel)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "68c3d687",
   "metadata": {},
   "source": [
    "Now that we've done spectral smoothing, we can resample the spectrum. We include the parameter `suppress_smooth_warning = True` in this step, because otherwise we'll see a warning that \"the data should be smoothed prior to resampling.\" Since we have already smoothed the data, this warning can be safely ignored."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "40868f73",
   "metadata": {},
   "outputs": [],
   "source": [
    "cube2vel_spectralresample = cube2vel_smooth.spectral_interpolate(cube1vel_cutout.spectral_axis, suppress_smooth_warning = True)\n",
    "cube2vel_spectralresample"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f9262c35",
   "metadata": {},
   "source": [
    "## Step 5: Spatial Smoothing and Reprojection\n",
    "\n",
    "Recall that our two cubes started with different beam sizes. Our goal is to match the beam sizes of the two cubes, and then to reproject the cubes onto the same spatial pixel grid. We can start to do this by printing out the beam sizes of the two cubes as they currently stand."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "426c22b0",
   "metadata": {},
   "outputs": [],
   "source": [
    "cube1vel_cutout.beam, cube2vel_spectralresample.beam"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4b323381",
   "metadata": {},
   "source": [
    "`cube1` again has the larger beam, so we'll smooth `cube2` to its resolution."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "87c13bd2",
   "metadata": {},
   "source": [
    "### Aside: mixed beams \n",
    "\n",
    "If cube1 and cube2 had different sized beams, but neither was clearly larger, we would have to convolve _both_ to a [common beam](https://radio-beam.readthedocs.io/en/latest/commonbeam.html#finding-the-smallest-common-beam). In this case, it's redundant and we could have just used cube1's beam, but this is the more general approach:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "eb9bdafc",
   "metadata": {},
   "outputs": [],
   "source": [
    "import radio_beam\n",
    "common_beam = cube1vel_cutout.beam.commonbeam_with(cube2vel_spectralresample.beam)\n",
    "common_beam"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "bb3132c8",
   "metadata": {},
   "outputs": [],
   "source": [
    "# # This commented-out approach works in older versions of radio-beam:\n",
    "# import radio_beam\n",
    "# common_beam = radio_beam.commonbeam.common_2beams(radio_beam.Beams(beams=[cube1vel_cutout.beam, cube2vel_spectralresample.beam]))\n",
    "# common_beam"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ce348003",
   "metadata": {},
   "outputs": [],
   "source": [
    "cube2vel_spatialspectralsmooth = cube2vel_spectralresample.convolve_to(common_beam)\n",
    "cube2vel_spatialspectralsmooth"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "6f2f2964",
   "metadata": {},
   "source": [
    "Now, we can reproject `cube2` onto the same spatial grid as `cube1`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "89ee0def",
   "metadata": {},
   "outputs": [],
   "source": [
    "# these next two lines are a hack to prevent the WCS from trying to convert to frequency\n",
    "tgt_header = cube1vel_cutout.header\n",
    "tgt_header['RESTFRQ'] = cube2vel_spatialspectralsmooth.header['RESTFRQ']\n",
    "\n",
    "cube2vel_reproj = cube2vel_spatialspectralsmooth.reproject(tgt_header)\n",
    "cube2vel_reproj "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "20464757",
   "metadata": {},
   "source": [
    "These two cubes are now on an identical grid, and can be directly compared:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "38863b71",
   "metadata": {},
   "outputs": [],
   "source": [
    "cube2vel_reproj, cube1vel_cutout"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "16576500",
   "metadata": {},
   "source": [
    "These spectra can now be overplotted as they are in the same unit with the same beam. Using the quicklook maps of the peak intensity from step 2, we know that there is some signal in the central region of our pointing, so we can plot the spectrum for the pixel there. We can approximate the pixel coordinate using the peak intensity map with fewer pixels, since we know both cubes got smoothed to that (coarser) spatial resolution. The center is at about (125, 125)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "116b9ea8",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.plot(cube1vel_cutout[:,125,125].spectral_axis, cube1vel_cutout[:,125,125], linestyle = '-', drawstyle = 'steps-mid', label = 'Cube 1')\n",
    "plt.plot(cube2vel_reproj[:,125,125].spectral_axis, cube2vel_reproj[:,125,125], linestyle = '-', drawstyle = 'steps-mid', label = 'Cube 2')\n",
    "plt.xlabel(cube1vel_cutout.spectral_axis.unit)\n",
    "plt.ylabel(cube1vel_cutout.unit)\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d4031ef2",
   "metadata": {},
   "source": [
    "### Do a quick continuum subtraction\n",
    "\n",
    "If we looked at the spectra we calculated in the previous step and saw that the baseline level of the spectra were not at zero, we would want to perform continuum subtraction to subtract off that non-zero baseline (since we're only working with the spectral lines for our purposes). Though our above spectra look like they are likely already continuum-subtracted, we can perform the following steps anyway. First, we'll calculate the median values of each of the cubes along the spectral axes. This gives us the approximate continuum value in each pixel. The median is more resistant to the presence of the spectral lines than the mean would be."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f773644a",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Calculate the median value along the spectral axis\n",
    "med1 = cube1vel_cutout.median(axis=0)\n",
    "med1.quicklook()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "37b08c88",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Calculate the median value along the spectral axis\n",
    "med2 = cube2vel_reproj.median(axis=0)\n",
    "med2.quicklook()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "77e7ba93",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Perform background subtraction on cubes by subtracting the median spectral value for each pixel\n",
    "cube1vel_cutout_bgsubtract = cube1vel_cutout - med1\n",
    "cube2vel_reproj_bgsubtract = cube2vel_reproj - med2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "89389c86",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot continuum-subtracted spectra for the central pixel\n",
    "plt.plot(cube1vel_cutout_bgsubtract[:,125,125].spectral_axis, cube1vel_cutout_bgsubtract[:,125,125], linestyle = '-', drawstyle = 'steps-mid', label = 'Cube 1')\n",
    "plt.plot(cube2vel_reproj_bgsubtract[:,125,125].spectral_axis, cube2vel_reproj_bgsubtract[:,125,125], linestyle = '-', drawstyle = 'steps-mid', label = 'Cube 2')\n",
    "plt.xlabel(cube1vel_cutout.spectral_axis.unit)\n",
    "plt.ylabel(cube1vel_cutout.unit)\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d003359f",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Plot difference between non-cont-sub and cont-sub spectra for the central pixel\n",
    "diff1 = cube1vel_cutout[:,125,125] - cube1vel_cutout_bgsubtract[:,125,125]\n",
    "diff2 = cube2vel_reproj[:,125,125] - cube2vel_reproj_bgsubtract[:,125,125]\n",
    "\n",
    "plt.plot(cube1vel_cutout_bgsubtract[:,125,125].spectral_axis, diff1, linestyle = '-', drawstyle = 'steps-mid', label = 'Cube 1')\n",
    "plt.plot(cube2vel_reproj_bgsubtract[:,125,125].spectral_axis, diff2, linestyle = '-', drawstyle = 'steps-mid', label = 'Cube 2')\n",
    "plt.xlabel(cube1vel_cutout.spectral_axis.unit)\n",
    "plt.ylabel(cube1vel_cutout.unit)\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a5be621f",
   "metadata": {},
   "source": [
    "For our spectra, since they are already continuum-subtracted, there is almost no visible change."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4f5c4698",
   "metadata": {},
   "source": [
    "## Step 6: Signal Masking\n",
    "\n",
    "Now, we will use one of the molecules, H2CS, to mask the other one, SiO.  This can be useful if SiO has complicated multi-component velocity structure, but H2CS has a similar morphology on the sky and we want to use that to grab the matching component of the SiO velocity structure.\n",
    "\n",
    "First, we'll estimate the noise throughout the cube by using a signal-resistant measurement of the standard deviation (called the median absolute deviation (MAD) standard deviation). We'll sigma-clip the cube first at 3 sigma and plot the spectrum to look for where the noise is least affected by signal."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "09a9664c",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Make signal mask out of template molecule cube, in noise-free area (should be flat)\n",
    "cube_sclip = cube1vel_cutout_bgsubtract.sigma_clip_spectrally(3)\n",
    "mad_std_spectrum_sclip = cube_sclip.mad_std(axis=(1, 2))\n",
    "plt.plot(mad_std_spectrum_sclip.spectral_axis.value, mad_std_spectrum_sclip.value, \n",
    "         drawstyle='steps-mid')\n",
    "plt.xlabel('Velocity (km/s)')\n",
    "plt.ylabel(r' Noise standard deviation $\\sigma$ (K)')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8f0d4b3e",
   "metadata": {},
   "source": [
    "We can see that between -15 and 45 km/s, the noise estimate is affected by signal. However, above 60 km/s, the noise is flatter. We'll use this part of the spectrum to estimate the noise map."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "52d5ce54",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Calculate noise map based on flat part of spectrum\n",
    "cube_sclip_cut = cube_sclip.spectral_slab(60*u.km/u.s, 100*u.km/u.s)\n",
    "mad_std_spectrum_sclip_cut = cube_sclip_cut.mad_std(axis=(1, 2))\n",
    "plt.plot(mad_std_spectrum_sclip_cut.spectral_axis.value, mad_std_spectrum_sclip_cut.value, \n",
    "         drawstyle='steps-mid')\n",
    "plt.xlabel('Velocity (km/s)')\n",
    "plt.ylabel(r' Noise standard deviation $\\sigma$ (K)')\n",
    "mad_std_map_sclip = cube_sclip_cut.mad_std(axis=0) # Calculate sigma along the spectral dimension"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ad6d11d6",
   "metadata": {},
   "source": [
    "Now, we'll create a \"plain mask\" where we mask the template cube and only keep pixels that are at or above 3 sigma (where sigma is given on a pixel-by-pixel basis by the noise map we just created). This procedure was originally illustrated in this [signal masking tutorial](https://mybinder.org/v2/gh/radio-astro-tools/tutorials/HEAD?filepath=masking_and_moments%2Fmasking_and_moments_guide.ipynb)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b3a089ce",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Calculate and apply plain mask\n",
    "plain_mask = cube1vel_cutout_bgsubtract >= 3 * mad_std_map_sclip # Get plain 3sigma mask\n",
    "plain_masked_slab = cube1vel_cutout_bgsubtract.with_mask(plain_mask) # Mask the template molecule cube"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "45382c82",
   "metadata": {},
   "outputs": [],
   "source": [
    "plain_masked_slab.shape"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "642b20a5",
   "metadata": {},
   "source": [
    "Then, we'll apply that mask to the template molecule cube and create a more sophisticated \"signal mask\" from that cube, which takes into account the fact that connected regions of signal are more likely to be real than an individual high-sigma pixel."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "da5bd982",
   "metadata": {},
   "outputs": [],
   "source": [
    "import scipy.ndimage as nd\n",
    "# Make a low and high mask\n",
    "low_snr_mask = (plain_masked_slab > 3 * mad_std_map_sclip).include()\n",
    "high_snr_mask = (plain_masked_slab > 10 * mad_std_map_sclip).include()\n",
    "# low_snr_mask = low_snr_mask.compute() # Don't need this for this tutorial\n",
    "# high_snr_mask = high_snr_mask.compute() # Don't need this for this tutorial\n",
    "# Find connected structures\n",
    "structure = np.ones((3, 3, 3), dtype=bool)\n",
    "low_snr_mask_labels, num_labels = nd.label(low_snr_mask, structure=structure)\n",
    "print(f\"Initial number of regions found: {num_labels}\")\n",
    "# From the labels, count the number of pixels within each label.\n",
    "num_pixels_in_high_snr_mask = nd.sum(high_snr_mask,\n",
    "                                     labels=low_snr_mask_labels,\n",
    "                                     index=range(1, num_labels + 1)) # +1 offset for mask labels\n",
    "# Repeat for the high signal mask.\n",
    "num_pixels_in_low_snr_mask = nd.sum(low_snr_mask,\n",
    "                                    labels=low_snr_mask_labels,\n",
    "                                    index=range(1, num_labels + 1)) # +1 offset for mask labels\n",
    "# To preserve the low_snr_mask, we will create a new signal mask where we will remove \n",
    "# regions that do not pass the criteria.\n",
    "signal_mask = low_snr_mask\n",
    "low_min_pixels = 40\n",
    "high_min_pixels = 10\n",
    "for num, (high_pix_num, low_pix_num) in enumerate(zip(num_pixels_in_high_snr_mask, \n",
    "                                                      num_pixels_in_low_snr_mask)):\n",
    "    if high_pix_num >= high_min_pixels and low_pix_num >= low_min_pixels:\n",
    "        # This region passes the criteria. Keep it in the mask.\n",
    "        continue\n",
    "    # Remove regions that do not pass the criteria.\n",
    "    # NOTE: enumerate will start with 0, but the mask labels start at 1\n",
    "    # We apply a +1 offset to `num` to account for this.\n",
    "    signal_mask[low_snr_mask_labels == num + 1] = False\n",
    "signal_mask_labels, num_labels = nd.label(signal_mask,\n",
    "                                          structure=structure)\n",
    "print(f\"Final number of regions found: {num_labels}\")"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7908389d",
   "metadata": {},
   "source": [
    "### Binary dilation\n",
    "\n",
    "The process of binary dilation grows the mask by one pixel outwards in every direction. This is meant to smooth the rough edges of the mask that may be generated when we perform these cuts, and it could help us to recover real contiguous signal that we may have accidentally removed in the masking process. However, it can also pick up noise. *User beware: compare the number of outliers in your final moment maps (this is done in a future step) with and without binary dilation to determine if you wish to use it in your processing.*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "40b7614b",
   "metadata": {},
   "outputs": [],
   "source": [
    "# signal_mask = nd.binary_dilation(signal_mask, structure=structure, iterations=1)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f9189b0b",
   "metadata": {},
   "source": [
    "### Mask application"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b2e820a4",
   "metadata": {},
   "source": [
    "Finally, we'll apply the signal mask we just created using the template molecule cube (which has a mask at each velocity step) to the target molecule cube."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3de87353",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Apply signal mask from template molecule to target molecule\n",
    "masked_cube = cube2vel_reproj_bgsubtract.with_mask(signal_mask)\n",
    "\n",
    "# # We can write the masked cube to a file and look at it in DS9, but this is optional:\n",
    "# masked_cube.write('masked_target_cube.fits', overwrite=True)\n",
    "# masked_cube = SpectralCube.read('masked_target_cube.fits')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7ab21ad2",
   "metadata": {},
   "source": [
    "Visualizing the masks:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "7db5505e",
   "metadata": {},
   "outputs": [],
   "source": [
    "plain_mask_comp = plain_mask.include() # Need to compute this mask before we visualize it\n",
    "plt.imshow(plain_mask_comp[30]) # Pick a channel\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "9027ad47",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.imshow(signal_mask[30]) # Pick that same channel\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "054d93a3",
   "metadata": {},
   "source": [
    "## Step 7: Moment Maps\n",
    "\n",
    "First, we'll create some initial moment maps of the target molecule using the cube we just masked. The negative values in these moment maps are due to us using a mask from one cube to mask signal in another cube. We can pick up noise that way, and since the noise can be negative, the sum along the spectral axis (which is what we are plotting in a moment 0 map) can also be negative."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f4b020f9",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Make moment maps from the target cube\n",
    "masked_moment0 = masked_cube.moment0()\n",
    "\n",
    "ax = plt.subplot(projection = masked_moment0.wcs)\n",
    "im = ax.imshow(masked_moment0.value, origin='lower', cmap='inferno')\n",
    "cbar = plt.colorbar(im)\n",
    "cbar.set_label('Integrated Intensity (K km/s)')\n",
    "ax.set_ylabel('Declination')\n",
    "ax.set_xlabel('Right Ascension')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5ae3fd7d",
   "metadata": {},
   "outputs": [],
   "source": [
    "masked_moment1 = masked_cube.moment1()\n",
    "\n",
    "ax = plt.subplot(projection = masked_moment1.wcs)\n",
    "im = ax.imshow(masked_moment1.value, origin='lower', cmap='coolwarm')\n",
    "cbar = plt.colorbar(im)\n",
    "cbar.set_label('Centroid (km/s)')\n",
    "\n",
    "ax.set_ylabel('Declination')\n",
    "ax.set_xlabel('Right Ascension')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f97344a9",
   "metadata": {},
   "source": [
    "The color bar on the moment 1 map is skewed due to the presence of at least one outlier pixel with a large negative velocity."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f58f1f2b",
   "metadata": {},
   "source": [
    "## Step 8: Masking Outliers\n",
    "\n",
    "At the moment, to get a better idea of the velocity structure of the target molecule, we can mask the outliers by masking pixels with a moment 0 value less than some threshold, which can cause moment 1 values to blow up since the moment 0 value appears in the denominator of the moment 1 value."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "84228efd",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Look for outliers\n",
    "v_thresh = 1000\n",
    "masked_moment1 = masked_cube.moment1()\n",
    "masked_moment1_outliers = (masked_moment1 > v_thresh*u.km/u.s)|(masked_moment1 < -v_thresh*u.km/u.s)\n",
    "plt.imshow(masked_moment1_outliers, origin='lower') \n",
    "# Clumps of outliers might mean they're real, just outside of vel range\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "0a643beb",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Look at spectra of outliers individually, can check neighboring spectra here too\n",
    "max_vel_coord = np.unravel_index(np.nanargmin(masked_moment1), masked_moment1.shape)\n",
    "spectrum = masked_cube[:, max_vel_coord[0], max_vel_coord[1]]\n",
    "print(masked_moment1[max_vel_coord[0], max_vel_coord[1]])\n",
    "plt.plot(spectrum.spectral_axis, spectrum.value, drawstyle='steps-mid')\n",
    "plt.xlabel('Velocity (km/s)')\n",
    "plt.ylabel('Intensity (K)')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "63fab9c6",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Remove outliers\n",
    "mom0 = masked_cube.moment0()\n",
    "mom0_mask = mom0 > 1.*u.K*u.km/u.s # Mask pixels with mom0 less than threshold\n",
    "print(f\"Found {mom0_mask.sum()} good pixels\")\n",
    "masked_cube_no_outliers = masked_cube.with_mask(mom0_mask)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "058bede7",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Make moment maps from the target cube with outliers removed\n",
    "masked_moment0 = masked_cube_no_outliers.moment0()\n",
    "\n",
    "ax = plt.subplot(projection = masked_moment0.wcs)\n",
    "im = ax.imshow(masked_moment0.value, origin='lower', cmap='inferno')\n",
    "cbar = plt.colorbar(im)\n",
    "cbar.set_label('Integrated Intensity (K km/s)')\n",
    "ax.set_ylabel('Declination')\n",
    "ax.set_xlabel('Right Ascension')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "e087f5c5",
   "metadata": {},
   "outputs": [],
   "source": [
    "masked_moment1 = masked_cube_no_outliers.moment1()\n",
    "\n",
    "ax = plt.subplot(projection = masked_moment1.wcs)\n",
    "im = ax.imshow(masked_moment1.value, origin='lower', cmap='coolwarm')\n",
    "cbar = plt.colorbar(im)\n",
    "cbar.set_label('Centroid (km/s)')\n",
    "\n",
    "ax.set_ylabel('Declination')\n",
    "ax.set_xlabel('Right Ascension')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "d3830e95",
   "metadata": {},
   "source": [
    "The velocity structure is more visible here with the outliers removed."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "db47bd10",
   "metadata": {},
   "outputs": [],
   "source": [
    "# masked_moment1.write('moment1.fits')"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "13856eee",
   "metadata": {},
   "source": [
    "### Addendum: collapse masked cube along spatial axes to make PV projections"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3a1d4da2",
   "metadata": {},
   "outputs": [],
   "source": [
    "masked_cube_no_outliers.shape # (75, 250, 250)\n",
    "collapse1 = masked_cube_no_outliers.sum(axis = 1) # Collapse along one spatial axis\n",
    "collapse2 = masked_cube_no_outliers.sum(axis = 2) # Collapse along the other spatial axis"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "8292acde",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.imshow(collapse1.value, origin='lower')\n",
    "plt.xlabel('Spatial axis (pixels)')\n",
    "plt.ylabel(f'Velocity axis [{masked_cube_no_outliers.spectral_axis.unit}]')\n",
    "plt.colorbar(label = '# of valid pixels along other spatial axis')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "fe870421",
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.imshow(collapse2.value, origin='lower')\n",
    "plt.xlabel('Spatial axis (pixels)')\n",
    "plt.ylabel(f'Velocity axis [{masked_cube_no_outliers.spectral_axis.unit}]')\n",
    "plt.colorbar(label = '# of valid pixels along other spatial axis')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7acf8ad1",
   "metadata": {},
   "source": [
    "This is another visualization, but it includes outliers."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "77d6dc62",
   "metadata": {},
   "outputs": [],
   "source": [
    "# Imports and matplotlib settings\n",
    "from astropy.visualization import simple_norm\n",
    "from mpl_toolkits.axes_grid1 import make_axes_locatable\n",
    "# plt.rcParams['figure.facecolor'] = 'w'\n",
    "\n",
    "# Collapse the mask in three dimensions\n",
    "mask = signal_mask.sum(axis=0)\n",
    "collapse1 = signal_mask.sum(axis = 1)\n",
    "collapse2 = signal_mask.sum(axis = 2)\n",
    "\n",
    "# Set up maximum value for colorbar\n",
    "vmax = np.max([mask.max(), collapse1.max(), collapse2.max()])\n",
    "norm = simple_norm(mask, stretch='linear', max_cut=vmax, min_cut=0.1)\n",
    "cm = plt.matplotlib.cm.viridis.copy()\n",
    "cm.set_under('w') # Make sure the \"zero\" color is white\n",
    "\n",
    "fig = plt.figure(figsize = (8.5, 10)) # constrained_layout=True\n",
    "\n",
    "# Get WCS coordinates from an (arbitrary?) masked cube\n",
    "masked_moment0 = masked_cube.moment0()\n",
    "wcs = masked_moment0.wcs\n",
    "\n",
    "# Collapse along axis 0\n",
    "ax1 = plt.subplot(323, projection = wcs, aspect = 1)\n",
    "ax1.set_aspect(1)\n",
    "im1 = ax1.imshow(mask, origin = 'lower', cmap=cm, vmin=0.1)\n",
    "ax1.tick_params(direction = 'in')\n",
    "ax1.coords[0].set_major_formatter('hh:mm:ss')\n",
    "ax1.coords[1].set_major_formatter('d.dd')\n",
    "ax1.set_xlabel('Right ascension'), ax1.set_ylabel('Declination')\n",
    "cbar1 = plt.colorbar(mappable = im1, ax = ax1)\n",
    "cbar1.set_label('# of valid pixels along axis')\n",
    "\n",
    "# This may not be right:\n",
    "# overlay = ax1.get_coords_overlay('galactic')\n",
    "# overlay.grid(color='white', ls='dotted')\n",
    "# overlay[0].set_axislabel('Galactic longitude')\n",
    "# overlay[1].set_axislabel('Galactic latitude')\n",
    "\n",
    "# Collapse along axis 1\n",
    "ax2 = plt.subplot(321, sharex = ax1) # , adjustable='box'\n",
    "im2 = ax2.imshow(collapse1, origin = 'lower', norm=norm, cmap=cm)\n",
    "ax2.xaxis.set_tick_params(labelbottom = False) # Remove redundant tick labels\n",
    "ax2.tick_params(direction = 'in')\n",
    "ax2.set_aspect(collapse1.shape[1]/collapse1.shape[0])\n",
    "ax2.set_ylabel('Velocity [km s$^{{-1}}$]')\n",
    "cbar2 = plt.colorbar(mappable=im2, ax=ax2)\n",
    "cbar2.set_label('# of valid pixels along axis')\n",
    "\n",
    "# Collapse along axis 2\n",
    "ax3 = plt.subplot(324, sharey = ax1)\n",
    "im3 = ax3.imshow(collapse2.T, origin = 'lower', norm=norm, cmap=cm)\n",
    "ax3.yaxis.set_tick_params(labelleft = False) # Remove redundant tick labels\n",
    "ax3.tick_params(direction = 'in')\n",
    "ax3.set_aspect(collapse2.shape[0]/collapse2.shape[1])\n",
    "ax3.set_xlabel('Velocity [km s$^{{-1}}$]')\n",
    "cbar3 = plt.colorbar(mappable=im3, ax=ax3)\n",
    "cbar3.set_label('# of valid pixels along axis')\n",
    "\n",
    "# plt.savefig(\"figures/mask.pdf\", dpi = 300, facecolor='w', edgecolor='w', bbox_inches='tight')\n",
    "plt.show()"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.10.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}