jdat_session/Cubeviz_Webbinar_solutions.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Cubeviz visualization tool"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Introduction\n",
    "\n",
    "This notebook introduces the basic functionality of the Cubeviz tool:\n",
    "\n",
    "* Launching the tool\n",
    "* Loading data from a notebook cell\n",
    "* Adjusting display parameters\n",
    "* Selecting spatial and spectral regions and passing them to the notebook\n",
    "* Extracting a spectrum with Cubeviz and passing it to the notebook\n",
    "* Using the Collapse plugin to extract an image from the cube\n",
    "* Displaying contours on an image\n",
    "* Extracting a spectrum from a notebook cell.\n",
    "\n",
    "The other data analysis plugins will be covered in a future Jwebbinar on IFU and MRS spectroscopy.\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Imports "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "* _numpy_ for array math\n",
    "* _scipy_ for ndcube gaussian smoothing\n",
    "* _specutils_ for Spectrum1D data model and cube manipulation\n",
    "* _jdaviz_ : Cubeviz data visualization tool\n",
    "* _photutils_ to define circular apertures\n",
    "* _astropy.io_ for reading and writing FITS cubes and images\n",
    "* _astropy.wcs, units, coordinates_ for defining and reading WCS\n",
    "* _astropy.stats_ for sigma_clipping\n",
    "* _astropy.utils_ for downloading files from URLs\n",
    "* _matplotlib_ for plotting spectra and images"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Resize notebook to full width\n",
    "from IPython.core.display import display, HTML\n",
    "display(HTML(\"<style>.container { width:100% !important; }</style>\"))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import time\n",
    "\n",
    "import numpy as np\n",
    "import scipy\n",
    "\n",
    "import specutils\n",
    "from specutils import Spectrum1D, SpectralRegion\n",
    "from specutils.manipulation import extract_region, spectral_slab\n",
    "\n",
    "from jdaviz import CubeViz\n",
    "\n",
    "from photutils import CircularAperture, SkyCircularAperture, aperture_photometry\n",
    "from photutils.detection import DAOStarFinder\n",
    "\n",
    "from regions import PixCoord, CirclePixelRegion\n",
    "\n",
    "from astropy.io import fits\n",
    "from astropy.stats import sigma_clip\n",
    "from astropy.stats import sigma_clipped_stats\n",
    "from astropy.utils.data import download_file\n",
    "import astropy.units as u\n",
    "\n",
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "from matplotlib.colors import LogNorm\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Simulated NIRSpec IFU Cube\n",
    "\n",
    "A point source (quasar) was simulated using the NIRSpec Instrument Performance Simulator (IPS), then run through the JWST Spec2 pipeline. We will use this for our science dataset."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Download the cube data from Box\n",
    "filename = \"https://stsci.box.com/shared/static/ff0bj31acot1272x5qq2clbmto3s6v4f.fits\"\n",
    "cube_file = download_file(filename, cache=True)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Visualize Science Data with Cubeviz"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "cubeviz = CubeViz()\n",
    "cubeviz.app"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### UI Instructions:\n",
    "#### Task 1:  Load the cube and view it\n",
    "* Load science datacube into Cubeviz using the next code cell below\n",
    "* Go to Tools icon: Gear icon: Layer in the leftmost image viewer \n",
    "* In that tab, change the stretch to Logarithmic, 99 percentile to see the target PSF wings.\n",
    "* Scrub through the cube using the Slice slider.\n",
    "\n",
    "#### Task 2:  Extract a spectrum from a region (subset) of the cube\n",
    "* Select a circular subset region centered on the source, using the region selection tool under the Tools icon. \n",
    "* Note that the region is pixelated and doesn't include fractional pixels.\n",
    "* Change the collapse method to \"Sum\" in spectrum viewer: Tools icon: Gear icon : Viewer. (Note that the collapse method is for the erroneously labeled 2nd \"x axis\")\n",
    "* Adjust the size and location of the region until you are happy with the extraction that shows up in the spectrum viewer.\n",
    "\n",
    "#### Task 3: Select a spectral region for further analysis\n",
    "* In the spectrum viewer, use the vertical and horizontal pan and zoom to focus on the broad H-alpha line (redshifted to 1.64 microns).  To pan: drag and drop spectrum with mouse or trackpad. To zoom: scroll with mouse or trackpad.\n",
    "* Select 'No selection (create new)' in the Subsets menu in the Cubeviz top banner.\n",
    "* Select a spectral subset region in spectrum viewer, using the region selection tool under the Tools icon, centered on a spectral feature of interest.\n",
    "\n",
    "#### Task 4: Extract an image by collapsing the cube\n",
    "* Find the Collapse data analysis Plugin by clicking the 'Lego' icon at upper right.\n",
    "* Select a dataset to collapse with the Data dropdown (e.g. SCI extension of our cube).\n",
    "* Enter the cube Axis to collapse (0 for the spectral axis).\n",
    "* Pick the 'Sum' collapse method.\n",
    "* Choose a spectral region to collapse in the Spectral Region dropdown (e.g. Subset 2).\n",
    "* Press 'Apply' \n",
    "* Use the 2nd cube/image viewer to display the collapsed data cube image. Go to Tools icon: Gear icon : Data dropdown, deselect the 'ERR' extension checkbox, and select the 'Collapsed' dataset that was generated by the previous step. \n",
    "* NB: you may need to resize the viewer(s) in order to see the full tools menu.\n",
    "* Adjust the Stretch in the viewer as needed.\n",
    "\n",
    "#### Task 5: Display contours on an image or cube\n",
    "* Go to Tools icon: Gear icon : Layer in the 2nd viewer\n",
    "* Make sure the 'SCI' extension is selected.\n",
    "* Change the stretch in of the collapsed cube image to Logarithmic 99%.\n",
    "* Click on the Contour 'eye with slash' icon to remove the slash and display the autogenerated contours.\n",
    "* Click on the wrench icon to enter custom contour levels (try 1e+12, 1e+13, 1e+14 for the collapsed data set). Note: There is a bug that autocorrects scientific notation. Try entering a junk character like at the beginning of the list (e.g. #1e+12, 1e+13, 1e+14) then delete the '#', as a temporary work-around."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Load Cube into Cubeviz"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "time.sleep(2)  # Sleep to avoid glue-jupyter timing issue\n",
    "cubeviz.app.load_data(cube_file)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Export Region from Cubeviz\n",
    "Export the region defined by the user in Cubeviz as an astropy CirclePixel Region, which has units of pixels.  We need \"try:, except:\" to catch the case where there is no Subset 1 selected in Cubeviz. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "cubeviz_data = cubeviz.app.data_collection[0]\n",
    "try:\n",
    "    region1 = cubeviz.app.get_subsets_from_viewer('spectrum-viewer')['Subset 1']\n",
    "    print(region1)\n",
    "    region1_exists = True\n",
    "    center1_xy = [region1.center.x, region1.center.y]  \n",
    "    r_pix = region1.radius\n",
    "\n",
    "except Exception:\n",
    "    print(\"There are no regions selected in the cube viewer.\")\n",
    "    region1_exists = False\n",
    "    center1_xy = [17.1, 20.]\n",
    "    r_pix = 6.0\n",
    "      "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Extract Subset Spectrum in Cubeviz Spectrum Viewer\n",
    "Retrieve the spectrum (Subset1) of the user-defined region from the Spectrum Viewer as a Spectrum1D object. Trim to remove bad wavelength ranges."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "wave_trim =[1.0*u.um,1.88*u.um]\n",
    "trim_region = SpectralRegion(wave_trim[0], wave_trim[1])\n",
    "try:\n",
    "    spectrum_subset1_untrimmed = cubeviz.app.get_data_from_viewer('spectrum-viewer')['Subset 1']\n",
    "    print(spectrum_subset1_untrimmed)\n",
    "    \n",
    "    # Trim the extracted spectrum\n",
    "    spectrum_subset1 = extract_region(spectrum_subset1_untrimmed, trim_region)\n",
    "    \n",
    "    print()\n",
    "    print('Wavelength:', spectrum_subset1_untrimmed.spectral_axis)\n",
    "    print('Trimmed:', spectrum_subset1.spectral_axis)\n",
    "    print()\n",
    "\n",
    "except Exception:\n",
    "    print(\"There are no subsets selected in the spectrum viewer.\")\n",
    "    \n",
    "\n",
    "#spectrum_subset1_untrimmed = spectrum_subset1_untrimmed.with_spectral_unit(u.um)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise: Extract Spectrum in Linearly Expanding Circular Aperture (Cone)\n",
    "This method is appropriate for a point source PSF with width proportional to wavelength"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Open and inspect the file with astropy.fits.open\n",
    "with fits.open(cube_file, memmap=False) as hdulist:\n",
    "    hdulist.info()\n",
    "    sci = hdulist['SCI'].data\n",
    "    err = hdulist['ERR'].data  \n",
    "    \n",
    "# Load original (untrimmed version) with Spectrum1D    \n",
    "spec1d_untrimmed = Spectrum1D.read(cube_file, format='JWST s3d')\n",
    "wavelength_untrimmed = spec1d_untrimmed.spectral_axis\n",
    "print()\n",
    "print(\"Wavelength: \", wavelength_untrimmed)\n",
    "\n",
    "# Trim the data cube and adjust region location\n",
    "x_trim = [2,-1]\n",
    "y_trim = [5, -4]\n",
    "spec1d = spectral_slab(spec1d_untrimmed, wave_trim[0], wave_trim[1])[x_trim[0]:x_trim[1],y_trim[0]:y_trim[1],:]\n",
    "wavelength = spec1d.spectral_axis\n",
    "spec1d_len = len(wavelength.value)\n",
    "\n",
    "#Adjust region location in trimmed cube\n",
    "center_trim = PixCoord(x=center1_xy[0]-x_trim[0], y=center1_xy[1]-y_trim[0])\n",
    "region_trim = CirclePixelRegion(center=center_trim, radius=r_pix)\n",
    "print(region_trim)\n",
    "\n",
    "center_x = center_trim.x\n",
    "center_y = center_trim.y\n",
    "\n",
    "#Use photutils.aperture_photometry to measure flux in expanding aperture\n",
    "cone_sum = []\n",
    "for idx in range(spec1d_len):\n",
    "    r_cone = r_pix * wavelength.value[idx]/ wavelength.value[0]\n",
    "    aperture_cone = CircularAperture((center_x,center_y), r=r_cone)\n",
    "    phot_table = aperture_photometry(spec1d.flux.value[:, :, idx], aperture_cone)\n",
    "    cone_sum.append(phot_table['aperture_sum'][0])\n",
    "    \n",
    "cone = Spectrum1D(flux=np.array(cone_sum)*u.MJy/u.sr, spectral_axis=spec1d.spectral_axis)\n",
    "print(cone)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot and compare spectra extracted in cone to Cubeviz subset:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "print(len(wavelength))\n",
    "print(len(spectrum_subset1.flux.value))\n",
    "f, (ax1) = plt.subplots(1, 1, figsize=(15, 5)) \n",
    "\n",
    "ax1.set_title(\"Spectral extractions\")\n",
    "ax1.set_xlabel(\"Observed Wavelength (microns)\")  \n",
    "ax1.set_ylabel(\"Flux Density\")\n",
    "ax1.set_xlim(0.99, 1.9)\n",
    "#ax1.set_ylim(0, 0.6)\n",
    "\n",
    "ax1.plot(wavelength, cone.flux.value, label=\"Cone\", c='darkorange', alpha=0.5)\n",
    "try:\n",
    "    ax1.plot(wavelength, spectrum_subset1.flux.value, c='r', label=\"Subset1\", alpha=0.4)\n",
    "except Exception:\n",
    "    print(\"There is no Cubeviz Subset1 spectrum to plot.\")\n",
    "ax1.legend()\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Comparison of the conical and Cubeviz subset spectral extractions. \n",
    "The conical extraction captures slightly more flux at long wavelengths.\n",
    "Red-shifted Broad Balmer and narrow [O III] lines  are visible in the quasar spectra. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<img style=\"float: right;\" src=\"https://raw.githubusercontent.com/spacetelescope/notebooks/master/assets/stsci_pri_combo_mark_horizonal_white_bkgd.png\" alt=\"Space Telescope Logo\" width=\"200px\"/>"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Notebook created by Patrick Ogle."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Jdaviz",
   "language": "python",
   "name": "jdaviz"
  },
  "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.9.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}