09b-Specutils/Specutils_overview_solutions.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# ![Specutils: An Astropy Package for Spectroscopy](data/specutils_logo.png)\n",
    "\n",
    "\n",
    "This notebook provides an overview of the Specutils Astropy coordinated package.  While this notebook is intended as an interactive introduction to specutils at the time of its writing, the canonical source of information for the package is the latest version's documentation: \n",
    "\n",
    "https://specutils.readthedocs.io"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Installation\n",
    "\n",
    " \n",
    "\n",
    "<div class=\"alert alert-block alert-warning\">\n",
    "If you followed the workshop instructions, specutils should be already installed. But depending on when you installed it, you might need to follow <a href=\"../../edit/00-Install_and_Setup/UPDATING.md\">the UPDATING.md instructions</a> at the base of this repo.\n",
    "</div>\n",
    "\n",
    "\n",
    "If it is not already installed, installing `specutils` should be as straightforward as either:\n",
    "```\n",
    "$ conda install -c conda-forge specutils\n",
    "```\n",
    "\n",
    "or\n",
    "\n",
    "```\n",
    "$ pip install specutils\n",
    "```\n",
    "\n",
    "although for futher information see [the installation section of the docs](https://specutils.readthedocs.io/en/latest/installation.html)."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Once these are installed, fundamental imports necessary for this notebook are possible:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "\n",
    "import astropy.units as u\n",
    "from astropy.coordinates import SpectralCoord\n",
    "\n",
    "import specutils\n",
    "from specutils import Spectrum1D\n",
    "specutils.__version__"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# for plotting:\n",
    "%matplotlib inline\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "\n",
    "# for showing quantity units on axes automatically:\n",
    "from astropy.visualization import quantity_support\n",
    "quantity_support();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Background/Spectroscopic ecosystem\n",
    "\n",
    "The large-scale plan for spectroscopy support in the Astropy project is outlined in  the [Astropy Proposal For Enhancement 13](https://github.com/astropy/astropy-APEs/blob/main/APE13.rst).  In summary, this APE13 lays out three broad packages:\n",
    "\n",
    "* `specutils` - a Python package containing the basic data structures for representing spectroscopic data sets, as well as a suite of fundamental spectroscopic analysis tools to work with these data structures.\n",
    "* `specreduce` - a general Python package to reduce raw astronomical spectral images to 1d spectra (represented as `specutils` objects).\n",
    "* `specviz` - a Python package (or possibly suite of packages) for visualization of astronomical spectra.\n",
    "\n",
    "\n",
    "While all are still in development, the first of these is furthest along, and is the subject of this notebook, as it contains the core data structures and concepts required for the others."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Fundamentals of specutils"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Objects for representing spectra\n",
    "\n",
    "The most fundamental purpose of `specutils` is to contain the shared Python-level data structures for storing astronomical spectra.  It is important to recognize that this is not the same as the *on-disk* representation.  As desecribed later specutils provides loaders and writers for various on-disk representations, with the intent that they all load to a common set of in-memory/Python interfaces.  Those intefaces (implemented as Python classes) are described in detail in the [relevant section of the documentation](https://specutils.readthedocs.io/en/latest/types_of_spectra.html), which contains this diagram:\n",
    "\n",
    "![Specutils Classes](data/specutils_classes_diagrams.png)\n",
    "\n",
    "The core principal is that all of these representations contain a `spectral_axis` attribute as well as a `flux` attribute (as well as optional matching `uncertainty`).  The former is often wavelength for OIR spectra, but might be frequency or energy for e.g. Radio or X-ray spectra.  Regardless of which spectral axis is used, the class attempts to interpret it appropriately, using the features of `astropy.Quantity` to distinguish different types of axes.  Similarly, `flux` may or may not be a traditional astronomical `flux` unit (e.g. Jy or  erg sec$^{-1}$ cm$^{-2}$ angstrom$^{-1}$), but is treated as the portion of the spectrum that acts in that manner.  The various classes are then distinguished by whether these attributes are one-dimensional or not, and how to map the `spectral_axis` dimensionality onto the `flux`.  The simplest case (and the one primarily considered here) is the scalar `Spectrum1D` case, which is a single spectrum with a matched-size `flux` and `spectral_axis`."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Basics of creating Spectrum1D Objects"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "While often a spectrum will be loaded from a file (see below), if the format is not compatible, or particular customization is required, spectra can be greated directly from arrays and astropy quantities.  E.g.:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "flux = (np.random.randn(100)*0.2 + 2) * u.erg / u.s / (u.cm**2) / u.AA \n",
    "wavelength = np.linspace(3000, 8000, 100)*u.AA\n",
    "spec1d = Spectrum1D(spectral_axis=wavelength, flux=flux)\n",
    "\n",
    "plt.step(spec1d.spectral_axis, spec1d.flux)\n",
    "\n",
    "spec1d"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "These can be straightforwardly transformed to other spectral units if desired:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "jyspec1d = spec1d.new_flux_unit(u.Jy)\n",
    "plt.plot(jyspec1d.frequency, jyspec1d.flux)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Working with Spectral Axes \n",
    "\n",
    "We created `Spectrum1D` just as Quantity arrays, so they can be treated just as `Quantity` objects when convenient with unit conversions and the like:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "spec1d.spectral_axis"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "spec1d.spectral_axis.to(u.nm)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "spec1d.spectral_axis.to(u.THz, u.spectral())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "But under the hood this are are fully-featured WCS following the [Astropy APE14](https://github.com/astropy/astropy-APEs/blob/main/APE14.rst) WCS interface along with the [GWCS](https://gwcs.readthedocs.io/) package. So you can use that to do conversions to and from spectral to pixel axes:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "spec1d.wcs.pixel_to_world([10, 10.5])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "spec1d.wcs.world_to_pixel(SpectralCoord([400, 450]*u.nm))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Uncertainties\n",
    "\n",
    "Currently the most compatible way to use uncertainties is the machinery built for the `astropy.nddata` object (although in many cases simply passing in an uncertainty array will also work):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from astropy.nddata import StdDevUncertainty\n",
    "\n",
    "unc = StdDevUncertainty(0.2 * np.ones(100))\n",
    "spec1d_unc = Spectrum1D(spectral_axis=wavelength, flux=flux, uncertainty=unc)\n",
    "\n",
    "# the errorbar matplotlib function does not natively support quantities\n",
    "plt.errorbar(spec1d_unc.wavelength.value, \n",
    "             spec1d_unc.flux.value, \n",
    "             spec1d_unc.uncertainty.array,\n",
    "             fmt='-',\n",
    "             drawstyle='steps-mid',\n",
    "             ecolor='k')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise\n",
    "\n",
    "Create a `Spectrum1D` object for an ideal 5800 K blackbody and plot it. Try the same, but with (random) noise added and stored as the uncertainty.\n",
    "\n",
    "Hint: while you can do this manually if you know the Planck function, there is an Astropy function to help you with this - you can find it via appropriate searches in the [astropy docs](http://docs.astropy.org)."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from astropy.modeling import models"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "wl = np.linspace(4000, 12000, 1024)*u.angstrom\n",
    "flux = models.BlackBody(temperature=5800*u.K)(wl)\n",
    "Spectrum1D(spectral_axis=wl, flux=flux)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "noise = np.random.randn(flux.size)*1e-6 * flux.unit\n",
    "bbspec = Spectrum1D(spectral_axis=wl, flux=flux+noise, uncertainty=StdDevUncertainty(1e-6*np.ones_like(flux)))\n",
    "bbspec"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.step(bbspec.spectral_axis, bbspec.flux)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Loading spectra from files\n",
    "\n",
    "Specutils also comes with readers for a variety of spectral data formats (including loaders for future JWST instruments). While support for specific formats depends primarily from users (like you!) providing readers, you may find that one has already been implemented for your favorite spectrum format.  As an example, we consider a spectrum form the [Sloan Digital Sky Survey](http://skyserver.sdss.org/):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from urllib.request import urlretrieve\n",
    "\n",
    "url = 'https://data.sdss.org/sas/dr16/sdss/spectro/redux/26/spectra/1323/spec-1323-52797-0012.fits'\n",
    "urlretrieve(url, 'data/sdss_spectrum.fits')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sdss_spec = Spectrum1D.read('data/sdss_spectrum.fits', format='SDSS-III/IV spec')\n",
    "plt.plot(sdss_spec.wavelength, sdss_spec.flux)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "To see the full list of formats readable in your current version of specutils, see the table at the bottom of the `Spectrum1d.read` method:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "help(Spectrum1D.read)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise\n",
    "\n",
    "If you have your own spectroscopic data, try loading a file here using either one of the built-in loaders, or the `Spectrum1D` interface, and plotting it.  If you don't have your own data on-hand, try downloading something of interest via a public archive (e.g., public HST data using MAST), and load it."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Note it is not really critical to use astroquery, one could do this from the MAST web site\n",
    "from astroquery.mast import Observations\n",
    "\n",
    "prods = Observations.get_product_list(Observations.query_criteria(proposal_id=1364))\n",
    "dls = Observations.download_products(prods[(prods['productSubGroupDescription'] == 'X1D') & (prods['type'] == 'S')][:1])\n",
    "spec_fn = dls['Local Path'][0]\n",
    "spec_fn"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "spec = Spectrum1D.read(spec_fn)\n",
    "spec"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.step(spec.spectral_axis, spec.flux)\n",
    "plt.ylim(0, 1e3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Manipulating spectra\n",
    "\n",
    "In addition to the analysis tools described in more detail in [the next notebook](Specutils_analysis.ipynb), Specutils also provides functionality for manipulating spectra.  In general these follow the pattern of creating *new* specutils objects with the results of the operation instead of in-place operations.\n",
    "\n",
    "The most straightforward of operations are arithmetic manipulations.  In general these follow expected patterns. E.g.:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "newspec1d = spec1d - 2 * u.erg / u.s / (u.cm**2) / u.AA\n",
    "plt.step(newspec1d.wavelength, newspec1d.flux)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "However, when there is ambiguity in your intent - for example, two spectra with different units where it is not clear what the desired output is - errors are generally produced instead of the code attempting to guess:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true,
    "tags": [
     "raises-exception"
    ]
   },
   "outputs": [],
   "source": [
    "# This raises an error:\n",
    "newspec1d - jyspec1d "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Resolving this requires explicit conversion:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "spec1d.new_flux_unit(u.Jy) - jyspec1d"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise\n",
    "\n",
    "Take the blackbody spectrum you generated above, and create a \"spectral feature\" by adding a Gaussian absorption or emission line to it using the arithmetic operators demonstrated above.  (Hint: [astropy.modeling](http://docs.astropy.org/en/stable/modeling/) contains [an implementation of the Gaussian line profile](http://docs.astropy.org/en/stable/api/astropy.modeling.functional_models.Gaussian1D.html#astropy.modeling.functional_models.Gaussian1D) which you may find useful.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gflux = models.Gaussian1D(amplitude=bbspec.flux.mean(), mean=6000, stddev=50)(bbspec.spectral_axis.value)\n",
    "\n",
    "bbspec_w_gaussian = bbspec + Spectrum1D(spectral_axis=bbspec.spectral_axis, flux=gflux)\n",
    "plt.step(bbspec_w_gaussian.spectral_axis, bbspec_w_gaussian.flux)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Smoothing, etc.\n",
    "\n",
    "More complex manipulation tools exist, outlined in the [relevant doc section](https://specutils.readthedocs.io/en/latest/manipulation.html). As a final example of this, we smooth our example spectrum using Gaussian smoothing:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from specutils.manipulation import gaussian_smooth\n",
    "\n",
    "smoothed_spec = gaussian_smooth(spec1d, 0.75)  # 0.75 pixel gaussian kernel\n",
    "plt.step(smoothed_spec.wavelength, smoothed_spec.flux)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercise\n",
    "\n",
    "Try smoothing your spectrum loaded in the above example (or the SDSS spectrum).  Compare all available kernel types and decide which seems most appropriate for your spectrum."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# we use the fake blackbody just to make sure this is clearly reproducible\n",
    "smoothed_spec = gaussian_smooth(bbspec_w_gaussian, 1)\n",
    "plt.step(smoothed_spec.spectral_axis, smoothed_spec.flux)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from specutils import manipulation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# we use the fake blackbody just to make sure this is clearly reproducible\n",
    "smoothed_spec = gaussian_smooth(bbspec_w_gaussian, 1)\n",
    "plt.step(smoothed_spec.spectral_axis, smoothed_spec.flux)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "smoothed_spec = manipulation.smoothing.median_smooth(bbspec_w_gaussian, 3)\n",
    "plt.step(smoothed_spec.spectral_axis, smoothed_spec.flux)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from astropy import convolution\n",
    "\n",
    "k = convolution.Trapezoid1DKernel(3)\n",
    "smoothed_spec = manipulation.smoothing.convolution_smooth(bbspec_w_gaussian, k)\n",
    "plt.step(smoothed_spec.spectral_axis, smoothed_spec.flux)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "k = convolution.kernels.CustomKernel([1,1,1,5,4,3,2])\n",
    "smoothed_spec = manipulation.smoothing.convolution_smooth(bbspec_w_gaussian, k)\n",
    "plt.step(smoothed_spec.spectral_axis, smoothed_spec.flux)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Which one to use is a bit of a matter of judgement, but you can decide from the plots above!"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Next Steps\n",
    "\n",
    "While the above focuses on loading and viewing spectra, the next notebook is aimed at quantitative scientific analysis of these spectra. Go to the [Specutils analysis notebook](Specutils_analysis.ipynb) for more on this."
   ]
  }
 ],
 "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.11.11"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}