diff --git a/tutorials/notebooks/specutils-fit-line/requirements.txt b/tutorials/notebooks/specutils-fit-line/requirements.txt new file mode 100644 index 00000000..80bc50ef --- /dev/null +++ b/tutorials/notebooks/specutils-fit-line/requirements.txt @@ -0,0 +1,6 @@ +astropy=4.2.1 +specutils=1.1 +numpy=1.20.2 +matplotlib=3.3.4 +astroquery=0.4.3 +ipython=7.22.0 diff --git a/tutorials/notebooks/specutils-fit-line/specutils-fit-line1.ipynb b/tutorials/notebooks/specutils-fit-line/specutils-fit-line1.ipynb new file mode 100644 index 00000000..d71de2eb --- /dev/null +++ b/tutorials/notebooks/specutils-fit-line/specutils-fit-line1.ipynb @@ -0,0 +1,476 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": { + "id": "uuznfl0QFMd_" + }, + "source": [ + "# Spectra 1: Load, Plot, and Smooth an Optical Spectrum from SDSS\n", + "\n", + "## Authors\n", + "Lia Corrales, Camilla Pacifici, Kelle Cruz\n", + "\n", + "## Learning Goals\n", + "* Download an SDSS spectrum using `astroquery`\n", + "* Load the spectrum with `specutils`\n", + "* Write a function to display the spectrum in `matplotlib`\n", + "* Smooth the spectrum with `specutils`\n", + "\n", + "## Keywords\n", + "\n", + "spectroscopy, specutils, astroquery\n", + "\n", + "## Companion Content\n", + "\n", + "This tutorial is Part 1 of a series:\n", + "* **Spectra 1: Load, Plot, and Smooth an Optical Spectrum from SDSS**\n", + "* Spectra 2: Measure Properties of a Spectral Emission Line\n", + "\n", + "These tutorials are based on the [MOS spectrum example](https://spacetelescope.github.io/jdat_notebooks/pages/mos-spectroscopy/mos-spectroscopy.html) by Camilla Pacifici.\n", + "\n", + "More information about optical SDSS spectra can be found on the [Getting Started with the Optical Data](https://www.sdss.org/dr16/spectro/spectro_basics/) SDSS help page.\n", + "\n", + "\n", + "## Summary\n", + "\n", + "In this tutorial, we will download the spectrum of a galaxy from the Sloan Digital Sky Survey (SDSS), load the spectrum as a Spectrum1D object using `specutils`, and then write a function to display the spectrum using `matplotlib`. We will use `specutils` to smooth the spectrum and plot it with the original." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "colab": { + "base_uri": "https://localhost:8080/", + "height": 358 + }, + "id": "3PhFPkSUFMeB", + "outputId": "2c03200d-0e80-429e-b61c-46d8e57d3a07" + }, + "outputs": [], + "source": [ + "import numpy as np\n", + "import matplotlib.pyplot as plt\n", + "\n", + "import astropy.units as u\n", + "\n", + "from astroquery.sdss import SDSS\n", + "from astropy.coordinates import SkyCoord\n", + "\n", + "from specutils import Spectrum1D\n", + "import specutils.manipulation" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "id": "S6AToKbJB4wM" + }, + "source": [ + "## Part 1: Download an SDSS spectrum" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "id": "Ooz_Hs2fG_4I" + }, + "source": [ + "Here, we will follow the `astroquery.sdss` [documentation](https://astroquery.readthedocs.io/en/latest/sdss/sdss.html) in identifying an example SDSS spectrum for download.\n", + "\n", + "First, we will query the SDSS archive for a quasar (QSO) spectrum at the following coordinates.\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "colab": { + "base_uri": "https://localhost:8080/", + "height": 153 + }, + "id": "326itk9JCRH9", + "outputId": "1dd5718e-b8e0-4865-a22e-486f6b70e6c6" + }, + "outputs": [], + "source": [ + "object_coords = SkyCoord(\"0h8m05.63s +14d50m23.3s\", frame='icrs')\n", + "\n", + "xid = SDSS.query_region(object_coords, spectro=True)\n", + "xid" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "id": "mNf1EQ6QE3YT" + }, + "source": [ + "The `SDSS.query_region` function returns an Astropy Table object that stores database identification information. This object, which we've named `xid`, can now be passed to `SDSS.get_spectra` to download the spectrum from the archive." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "id": "KTiZ_zrOCsVW" + }, + "outputs": [], + "source": [ + "sp_hdus = SDSS.get_spectra(matches=xid)" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "id": "JbEs3seNFj0e" + }, + "source": [ + "The `SDSS.get_spectra` method returns a list of [HDUList](https://docs.astropy.org/en/stable/io/fits/api/hdulists.html) objects, so that each item of the list contains the contents of FITS file. Because our query only returns one object, there is one HDUList returned by `SDSS.get_spectra`.\n", + "\n", + "We can treat the contents of `sp_hdus` as if we loaded a FITS file. Check out the [Viewing and manipulating data from FITS tables tutorial](https://learn.astropy.org/rst-tutorials/FITS-tables.html?highlight=fits) to learn more about working with FITS files. For example, we can print the information about the content contained in the HDUList:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "colab": { + "base_uri": "https://localhost:8080/" + }, + "id": "KtjoDpalC04j", + "outputId": "190532bb-0f96-46bb-871d-341c3f8107b5" + }, + "outputs": [], + "source": [ + "print(\"Number of HDUList objects:\", len(sp_hdus))\n", + "\n", + "print(sp_hdus[0].info())" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Part 2: Load the SDSS spectrum into a `Spectrum1D` object" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "id": "NWb7IoCwLYch" + }, + "source": [ + "We will load the spectral information into a `specutils.Spectrum1D` object, which is the simplest case for using specutils." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Method 1: Using the `specutils` built-in loader for SDSS spectrum files" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "id": "5wKmayM7TtHk" + }, + "source": [ + "Specutils has [many loaders](https://specutils.readthedocs.io/en/stable/spectrum1d.html#list-of-loaders) for common spectra formats, including one for SDSS spectra. If you have downloaded the SDSS spectrum FITS files, you can provide the filename to `Spectrum1D.read` and set the keyword argument `format=\"SDSS-III/IV spec\"`." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "qso0008_spectrum = Spectrum1D.read(sp_hdus[0], format=\"SDSS-III/IV spec\")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Method 2: Creating a `Spectrum1D` object manually" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "id": "QWbIferhHjFy" + }, + "source": [ + "Based on the [SDSS data model for spectra](https://data.sdss.org/datamodel/files/BOSS_SPECTRO_REDUX/RUN2D/spectra/PLATE4/spec.html), we suspect that the \"COADD\" extension holds the spectral information.\n", + "\n", + "To make sure this is true, we will extract the table contents to a new variable and examine the column names.\n", + "\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "colab": { + "base_uri": "https://localhost:8080/" + }, + "id": "mmEZuQA3IVn7", + "outputId": "4dfbf13b-6c00-4f30-d7da-1ac6e74ce5fc" + }, + "outputs": [], + "source": [ + "qso0008_coadd_data = sp_hdus[0][\"COADD\"].data\n", + "\n", + "print(qso0008_coadd_data.columns)" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "id": "y2BaxjeWJNwC" + }, + "source": [ + "Let's plot the columns that look like they might contain the spectrum information we want (wavelength and flux). We assume the 'flux' column contains the flux. Since the column name 'loglam' indicates that it contains the logarithim of lambda, a common name choice for wavelength, we'll plot that column after solving for lambda:\n", + "\n", + "$$ \\lambda = 10^{\\log \\lambda} $$" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "colab": { + "base_uri": "https://localhost:8080/", + "height": 283 + }, + "id": "rzOXuMqzH7JT", + "outputId": "cece376e-cc98-4c05-d01c-3ce1bbe1f971" + }, + "outputs": [], + "source": [ + "lambda = np.power(10.0, qso0008_coadd_data['loglam'])\n", + "plt.plot(lambda, qso0008_coadd_data['flux'])" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "id": "4NooTts_Jd2n" + }, + "source": [ + "Excellent, this looks like a QSO spectrum! \n", + "\n", + "Based on the [Getting Started with the Optical Data](https://www.sdss.org/dr16/spectro/spectro_basics/) SDSS help page, the SDSS spectrograph covers the wavelength range 3800 - 9200 Angstrom. So the wavelength units for the 'loglam' column must be Angstrom. The website also states that the flux units are $10^{-17}$erg cm$^{-2}$ s$^{-1}$ Angstrom$^{-1}$." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "colab": { + "base_uri": "https://localhost:8080/", + "height": 181 + }, + "id": "2zDm8iezKrcD", + "outputId": "ec9decb2-43a2-493d-d40f-9a36f87c855a" + }, + "outputs": [], + "source": [ + "qso0008_spectrum = Spectrum1D(flux=qso0008_coadd_data['flux'] * 1.e-17 * u.Unit('erg cm^-2 s^-1 AA^-1'), \n", + " spectral_axis=np.power(10.0, qso0008_coadd_data['loglam']) * u.angstrom)" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "id": "GrklgIDjLRax" + }, + "source": [ + "## Part 3: Write a function to plot the spectrum" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "id": "t9qoGTcDJdIp" + }, + "source": [ + "While it's possible to plot the spectrum directly by passing information stored in the `Spectrum1D` object to `matplotlib.pyplot.plot`, we will write a function for this process so that we can re-plot any `Spectrum1D` with ease." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "id": "pDsfwhbjLAXD" + }, + "outputs": [], + "source": [ + "def plot_spectrum1d(ax, spec, **kwargs):\n", + " \"\"\"\n", + " Plot a Spectrum1D object\n", + "\n", + " ax : matplotlib.pyplot.axis object on which to plot the data\n", + "\n", + " spec : specutils.Spectrum1D object to plot\n", + "\n", + " **kwargs will be passed to the matplotlib.pyplot.plot command\n", + " \"\"\"\n", + " ax.plot(spec.spectral_axis, spec.flux, **kwargs)\n", + " ax.set_xlabel(spec.spectral_axis.unit)\n", + " ax.set_ylabel(spec.flux.unit)\n", + " return" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "id": "_Oky3XBMLBof" + }, + "source": [ + "Next, set up the matplotlib axis and pass it to our newly defined function. We will then change the x-axis range in order to zoom in on some of the sharp emission features." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "id": "uEHoPDaSK_FD" + }, + "outputs": [], + "source": [ + "ax = plt.subplot(111)\n", + "\n", + "plot_spectrum1d(ax, qso0008_spectrum)\n", + "plt.xlim(5100, 5300)" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "id": "Ws1Bls9aNGo8" + }, + "source": [ + "Adding plotting functionality to specutils is currently under discussion and you can [read the details and contribute to the conversation on Github](https://github.com/astropy/specutils/pull/650)." + ] + }, + { + "cell_type": "markdown", + "metadata": { + "id": "6gwBkRTfMgO4" + }, + "source": [ + "## Part 4: Smooth the spectrum with `specutils` and re-plot" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "id": "_V-o2sQhNhum" + }, + "source": [ + "We can use a Gaussian kernel to smooth the spectrum data with the `specutils.manipulation.gaussian_smooth` function.\n", + "\n", + "The standard deviation of the Gaussian function is defined in units of pixels (i.e., number of bins), with the `stddev` keyword." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "id": "2o7GWoINLkZT" + }, + "outputs": [], + "source": [ + "qso0008_smoothed_spectrum = specutils.manipulation.gaussian_smooth(qso0008_spectrum, stddev=5)" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "id": "0BOWHVWUMyms" + }, + "source": [ + "Now, using the function we defined above, let's plot our two spectra together so that we can see the difference between the original and smoothed spectrum." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "id": "JdCKwqGjMu0Z" + }, + "outputs": [], + "source": [ + "ax = plt.subplot(111)\n", + "\n", + "plot_spectrum1d(ax, qso0008_spectrum, color='k', label='Original')\n", + "plot_spectrum1d(ax, qso0008_smoothed_spectrum, color='r', alpha=0.8, label='Smoothed')\n", + "\n", + "plt.legend(loc='upper right')\n", + "plt.xlim(5100, 5300)" + ] + }, + { + "cell_type": "markdown", + "metadata": { + "id": "Y_kYkFK3OXlS" + }, + "source": [ + "## Exercise\n", + "\n", + "Modify the plotting function so you can plot the spectrum with a different x-axis unit, such as frequency. Your function definition might look like:\n", + "\n", + "`def plot_spectrum1d(ax, spec, xunit='Angstrom'):`\n", + "\n", + "**Hint:** Check out the [Astropy units equivalencies module](https://docs.astropy.org/en/stable/units/equivalencies.html#spectral-units) to learn how to manage conversions for spectral units." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "id": "izWtOHpTNHT_" + }, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "id": "NC6LhuwMNdkk" + }, + "outputs": [], + "source": [] + } + ], + "metadata": { + "colab": { + "collapsed_sections": [], + "name": "specutils-fit-line1.ipynb", + "provenance": [] + }, + "kernelspec": { + "display_name": "astropy-tutorials", + "language": "python", + "name": "astropy-tutorials" + }, + "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.10" + } + }, + "nbformat": 4, + "nbformat_minor": 1 +}