{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Data Analysis Tools JWebbinar: Specutils\n", "\n", "![Specutils: An Astropy Package for Spectroscopy](specutils_logo.png)\n", "\n", "\n", "This notebook provides an overview of the Astropy coordinated package `specutils`. 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": [ "`specutils` should already be in your JWebbinar environment. If you wish to install locally, you can follow the instructions in [the installation section of the specutils docs](https://specutils.readthedocs.io/en/latest/installation.html).\n", "\n", "Once specutils is 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.io import fits\n", "\n", "import specutils\n", "from specutils import Spectrum1D, SpectralRegion, analysis, manipulation, fitting\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. While Astropy has not developed this, STScI has built [`jdaviz`](https://jdaviz.readthedocs.io/en/latest/) for this purpose (especially, but not exclusively, for JWST).\n", "\n", "\n", "The first is the subject of this notebook." ] }, { "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](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\n", "\n", "If your spectrum is in a format that specutils understands, loading it is very straightforward. There are times when you want a bit more control, though, so lets look at loading from a file and creating an object directly from arrays:" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Loading spectra from files\n", "\n", "Specutils 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 simulated high-redshift (z > 1) galaxy like that you might see from NIRSpec:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "url_to_download = 'https://stsci.box.com/shared/static/b22b1fzhimtdqfp8597m4bg67kovvauu.fits'" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from astropy.utils.data import download_file\n", "\n", "spec_fn = download_file(url_to_download, cache=True)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "jwst_spec = Spectrum1D.read(spec_fn)\n", "plt.plot(jwst_spec.spectral_axis, jwst_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": [ "# Creating a Spectrum \"by-hand\"" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "If you have a format that is not compatible, or you want to do some sort of customization of the loading process, spectra can be created directly from aastropy quantities (which are basically arrays with associated units). Here we'll show how you can do that, using data from the same file above: " ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "fitsfile = fits.open(spec_fn)\n", "fitsfile.info()" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "fitsfile[1].header" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "table_data = fitsfile[1].data\n", "table_data" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "flux = table_data['flux']*u.uJy\n", "wavelength = table_data['wavelength']*u.um\n", "\n", "spec1d_byhand = Spectrum1D(spectral_axis=wavelength, flux=flux)\n", "\n", "plt.step(spec1d_byhand.spectral_axis, spec1d_byhand.flux)\n", "\n", "spec1d_byhand" ] }, { "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, you can try downloading something of interest via a public archive (e.g., public HST data using MAST), or you can try with an [SDSS galaxy spectrum](https://dr14.sdss.org/optical/spectrum/view/data/format=fits/spec=lite?plateid=1323&mjd=52797&fiberid=12) included in this repo (look for the file `example_sdss_spectrum.fits`). Once you've got your dataset, try loading it into a `Spectrum1D` object whichever of the ways above makes most sense to you.\n", "\n", "If you have time and enough knowledge of the format, try loading your spectrum *both* ways - with a built-in loader and manually creating a `Spectrum1D`." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Working with Units and 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": [ "jwst_spec.spectral_axis" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "jwst_spec.spectral_axis.to(u.angstrom)" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "jwst_spec.spectral_axis.to(u.THz, u.spectral())" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "There's even fast-accessors to make some of this more convenient:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "plt.step(jwst_spec.frequency, jwst_spec.flux)" ] }, { "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": [ "jwst_spec.wcs.pixel_to_world([10, 10.5])" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "jwst_spec.wcs.world_to_pixel([1.2, 1.3]*u.um)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The other dimension is the flux - that requires slightly more complex transformation because it matters where in the spectrum you are to do the flux transformation:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "f_lamb_units = u.erg / u.s / (u.cm**2) / u.um\n", "\n", "jwst_spec.flux.to(f_lamb_units, u.spectral_density(jwst_spec.spectral_axis))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "But specutils makes this much easier!" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "jwst_spec_flamb = jwst_spec.new_flux_unit(f_lamb_units)\n", "jwst_spec_flamb" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "plt.plot(jwst_spec_flamb.spectral_axis, jwst_spec_flamb.flux)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Arithmetic on Spectra\n", "\n", "Specutils provides a lot of 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 patterns that are based on fundametal arithmetic. E.g.:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "newspec1d = jwst_spec - 0.5*u.uJy\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": "raw", "metadata": { "scrolled": true, "tags": [ "raises-exception" ] }, "source": [ "# This raises an error:\n", "newspec1d - jwst_spec_flamb " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Resolving this requires explicit conversion:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "newspec1d.new_flux_unit(f_lamb_units) - jwst_spec_flamb" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Exercises\n", "\n", "1. Try continuum *normalizing* using the above assumed continuum level, instead of subtracting. \n", "2. You can tell by-eye that the continuum isn't a flat level, but has a slope to it. Try subtracting a *sloping* continuum instead of a single value (it's fine to just estimate it by-eye, fancier continuum subtracting will come later)." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## SNR and Uncertainties\n", "\n", "Now lets try a simple calculation: the S/N of this spectrum. While pipeline-output JWST files will have uncertainties, for this example we are using a basic simulation without the uncertainities. Hence we start using an SNR estimate that follows a [straightforward algorithm detailed in the literature](https://www.stecf.org/software/ASTROsoft/DER_SNR/)." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "analysis.snr_derived(jwst_spec)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "That's straightforward enough, but does not use any uncertainty information beyond the spectrum itself. To mock up what real JWST might look like (or how you would provide your *own* uncertainties), we can create a new Spectrum1D object with by-hand added uncertainties:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from astropy.nddata import StdDevUncertainty\n", "\n", "unc = StdDevUncertainty(0.05 * np.ones_like(jwst_spec.flux))\n", "spec1d_unc = Spectrum1D(spectral_axis=jwst_spec.spectral_axis, \n", " flux=jwst_spec.flux, \n", " uncertainty=unc)\n", "\n", "plt.fill_between(spec1d_unc.spectral_axis, \n", " spec1d_unc.flux - spec1d_unc.uncertainty.quantity, \n", " spec1d_unc.flux + spec1d_unc.uncertainty.quantity);" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "After that, it is a one-liner to compute the S/N directly, which will use the uncertainty already in the spectrum:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "analysis.snr(spec1d_unc)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Additional Materials\n", "\n", "From here on this notebook provides more in-depth details about the. While there is likely not enough time to work through this during the Webbinar, it is provided for you to work through whichever parts might appeal to your workflows. This is the intended mode of working with `specutils`: it is a toolbox for you to build your own science cases!\n" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "# Sample Analysis: Line Center for Redshift\n", "\n", "Specutils has a large set of analysis functions, more than we have time to cover here. So instead we focus on a specific concrete science case: determining the center of a line and using that to get a redshift estimate.\n", "\n", "If you are used to looking at galaxies, you probably recognized that our example spectrum has a characteristic emission line pattern with a strong H$\\alpha$ line. While we might be able to center up on a line that strong, to be sure we should always subtract the continuum first. Many other `specutils` analaysis function expect continuum-subtracted spectra as well." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# note this doesn't error because with a *scalar* it's unambiguous that the user wants the spectrum units\n", "jwst_continuum = 550*u.nJy \n", "\n", "jwst_halpha_contsub = jwst_spec - jwst_continuum\n", "\n", "plt.step(jwst_halpha_contsub.spectral_axis, jwst_halpha_contsub.flux)\n", "plt.axhline(0, c='k', ls=':')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "While there are techniques for identifying the line automatically (see the fitting section below), here we assume we are doing \"quick-look\" procedures where manual identification is possible. \n", "\n", "We do this by defining a `SpectralRegion` object spanning the area of the line. Further below is a worked example of how to get these values easily in a notebook, but for now we can take the numbers I've pre-eyeballed: $1.638 \\mu m$ to $1.644 \\mu m$" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "halpha_lines_region = SpectralRegion(16380*u.angstrom, 16440*u.angstrom)\n", "\n", "plt.step(jwst_halpha_contsub.spectral_axis, jwst_halpha_contsub.flux)\n", "\n", "yl1, yl2 = plt.ylim()\n", "plt.fill_between([halpha_lines_region.lower, halpha_lines_region.upper], \n", " -1, 10, alpha=.2)\n", "plt.xlim(1.6, 1.7);\n", "plt.ylim(-.2, 2.5)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "You can now call a variety of analysis functions on the continuum-subtracted spectrum to estimate various properties of the line (you can see the full list of relevant analysis functions [in the analysis part of the specutils docs](https://specutils.readthedocs.io/en/stable/analysis.html#functions)). Here we do just the center:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "cen = analysis.centroid(jwst_halpha_contsub, halpha_lines_region)\n", "cen" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now we simply plug this into the redshift formula for H$\\alpha$'s rest wavelength:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "(cen/(6563*u.angstrom)) - 1" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We now have strong evidence this is a manufactured spectrum! (it's suspiciously close to *exactly* 1.5 ...)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Exercise\n", "\n", "Look for other star forming galaxy emission lines. Get their center wavelengths as well, and use those multiple lines to improve your redshift estimate (and perhaps an uncertainty on your estimate)." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Spectral Regions\n", "\n", "Most analysis required on a spectrum requires specification of a part of the spectrum - e.g., a spectral line. Because such regions may have value independent of a particular spectrum, they are represented as objects distrinct from a given spectrum object. Below we outline a few ways such regions are specified." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "ha_region = SpectralRegion((6563-50)*u.AA, (6563+50)*u.AA)\n", "ha_region" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Regions can also be raw pixel values (although of course this is more applicable to a specific spectrum):" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "ha_pixel_region = SpectralRegion(2650*u.pixel, 2850*u.pixel)\n", "ha_pixel_region" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Additionally, *multiple* regions can be in the same `SpectralRegion` object. This is useful for e.g. measuring multiple spectral features in one call:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "HI_wings_region = SpectralRegion([(1.44*u.GHz, 1.43*u.GHz), (1.41*u.GHz, 1.4*u.GHz)])\n", "HI_wings_region" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "These regions are used for many other analysis steps, including, for example, specifying a \"featureless\" region to use to estimate a noise level to assume for the spectrum:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "noise_region = SpectralRegion(1.4*u.micron, 1.5*u.micron)\n", "jwst_spec_with_unc = manipulation.noise_region_uncertainty(jwst_spec, noise_region)\n", "jwst_spec_with_unc" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "With that uncertainty present, we can now compute the S/N:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "analysis.snr(jwst_spec_with_unc)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Regions can also be used to they can be used to extract sub-spectra from larger spectra, which can be used to do yet further operations (like in this case, per-pixel S/N):" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "jwst_spec_with_unc_ha = manipulation.extract_region(jwst_spec_with_unc, ha_pixel_region)\n", "plt.step(jwst_spec_with_unc_ha.spectral_axis, \n", " jwst_spec_with_unc_ha.flux/jwst_spec_with_unc_ha.uncertainty.quantity)\n", "\n", "analysis.snr(jwst_spec_with_unc_ha)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Line Measurements\n", "\n", "While line-fitting (detailed more below) is a good choice for high signal-to-noise spectra or when detailed kinematics are desired, more empirical measures are often used in the literature for noisier spectra or just simpler analysis procedures. Specutils provides a set of functions to provide these sorts of measurements, as well as similar summary statistics about spectral regions. The [analysis part of the specutils documentation](https://specutils.readthedocs.io/en/latest/analysis.html) provides a full list and detailed examples of these, but here we demonstrate some example cases." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Note: these line measurements generally assume your spectrum is continuum-subtracted or continuum-normalized. Some spectral pipelines do this for you, but often this is not the case. For our examples here we will do this step \"by-eye\", but for a more detailed discussion of continuum modeling, see the next section. Based on the plots above we estimate a continuum level for the area of the simulated JWST spectrum around the H-alpha emission line, and use basic math to construct the continuum-normalized and continuum-subtracted spectra." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "# estimate a reasonable continuum-level estimate for the h-alpha area of the spectrum\n", "jwst_continuum = 550*u.nJy # note this is more convenient units than the spectrum itself\n", "\n", "jwst_halpha_contsub = manipulation.extract_region(jwst_spec, ha_pixel_region) - jwst_continuum\n", "# THE TWO LINES BELOW ARE A WORK_AROUND FOR A BUG IN SPECUTILS 1.5.0 AND CAN BE REMOVED IN FAVOR OF THE ABOVE IF YOU'RE ON A VERSION WITHOUT THE BUG\n", "jwst_halpha_ex = manipulation.extract_region(jwst_spec, ha_pixel_region)\n", "jwst_halpha_contsub = Spectrum1D(spectral_axis=jwst_halpha_ex.spectral_axis , flux=jwst_halpha_ex.flux- jwst_continuum) \n", "\n", "plt.axhline(0, c='k', ls=':')\n", "plt.step(jwst_halpha_contsub.spectral_axis, jwst_halpha_contsub.flux)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "With the continuum level identified, we can now make some measurements of the spectral lines that are apparent by eye - in particular we will focus on the H-alpha emission line. While there are techniques for identifying the line automatically (see the fitting section below), here we assume we are doing \"quick-look\" procedures where manual identification is possible. \n", "\n", "In the cell below, change the values for `LOWER` and `UPPER` to make a spectral region that just encompasses the H-alpha line (the middle of the three lines). You may find it useful to change the values, re-run the cell, and change again to \"hone in\" on the right number." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "LOWER = 16000 * u.angstrom\n", "UPPER = 17000 * u.angstrom\n", "halpha_lines_region = SpectralRegion(LOWER, UPPER)\n", "\n", "plt.step(jwst_halpha_contsub.spectral_axis, jwst_halpha_contsub.flux)\n", "\n", "yl1, yl2 = plt.ylim()\n", "plt.fill_between([halpha_lines_region.lower, halpha_lines_region.upper], \n", " yl1, yl2, alpha=.2)\n", "plt.ylim(yl1, yl2)\n", "plt.xlim(1.6, 1.7);" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "You can now call a variety of analysis functions on the continuum-subtracted spectrum to estimate various properties of the line (you can see the full list of relevant analysis functions [in the analysis part of the specutils docs](https://specutils.readthedocs.io/en/stable/analysis.html#functions)). Here we highlight the line center, flux, and equivalent width." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "cen = analysis.centroid(jwst_halpha_contsub, halpha_lines_region)\n", "cen" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The line flux is equally straightforward, although deceptively complicated:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "analysis.line_flux(jwst_halpha_contsub, halpha_lines_region)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "While technically correct, those are not necessarily the kinds of units you would need to compare with, say, optical line fluxes. Fortunately, `Spectrum1D` objects provide straightforward unit conversion:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "f_lamb_unit = u.erg*u.m**-2*u.s**-1*u.micron**-1\n", "jwst_halpha_contsub_flamb = jwst_halpha_contsub.new_flux_unit(f_lamb_unit)\n", "\n", "analysis.line_flux(jwst_halpha_contsub_flamb, halpha_lines_region).to(u.erg/u.s/u.cm**2)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Equivalent width, being a continuum dependent property, can be computed directly from the spectrum if the continuum level is given:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "analysis.equivalent_width(jwst_spec, jwst_continuum, regions=halpha_lines_region).to(u.angstrom)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "As per the normal convention, it is negative for an emission line." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Exercises" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Suppose you wanted to compare some of the above line measurements to a relatively local galaxy that has known line Hα line fluxes - say, some galaxy you know to be at 100 Mpc. How can you do this comparison in just a few lines of code? (Hint: making use of the units in the ``line_flux(...)`` outputs and [astropy.cosmology](https://docs.astropy.org/en/stable/cosmology/index.html) are your friends here!) If you'd like to try this with some real data, you might compare to a ground-based spectrum of a local emission line galaxy, [like the SDSS](https://dr14.sdss.org/optical/spectrum/view/data/format=fits/spec=lite?plateid=1323&mjd=52797&fiberid=12)." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Load one of the spectrum datasets you made in the overview exercises into this notebook (i.e., your own dataset, a downloaded one, or the blackbody with an artificially added spectral feature). Make a flux or width measurement of a line in that spectrum directly (i.e. without subtracting a continuum). Is anything odd?" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Continuum Subtraction\n", "\n", "While continuum-fitting for spectra is sometimes thought of as an \"art\" as much as a science, specutils provides the tools to do a variety of approaches to continuum-fitting, without making a specific recommendation about what is \"best\" (since it is often very data-dependent). More details are available [in the relevant specutils doc section](https://specutils.readthedocs.io/en/latest/fitting.html#continuum-fitting), but here we outline the two basic options as it stands: an \"often good-enough\" function, and a more customizable tool that leans on the [`astropy.modeling`](http://docs.astropy.org/en/stable/modeling/index.html) models to provide its flexibility." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### The \"often good-enough\" way\n", "\n", "The `fit_generic_continuum` function provides a function that is often sufficient for reasonably well-behaved continuua, particular for \"quick-look\" or similar applications where high precision is not that critical. The function yields a continuum model, which can be evaluated at any spectral axis value:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from specutils.fitting import fit_generic_continuum" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "generic_continuum = fit_generic_continuum(jwst_spec)\n", "\n", "generic_continuum_evaluated = generic_continuum(jwst_spec.spectral_axis)\n", "\n", "plt.step(jwst_spec.spectral_axis, jwst_spec.flux)\n", "plt.plot(jwst_spec.spectral_axis, generic_continuum_evaluated)\n", "plt.ylim(0, 1.5);" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "(Note that in some versions of astropy/specutils you may see a warning that the \"Model is linear in parameters\" upon executing the above cell. This is not a problem unless performance is a serious concern, in which case more customization is required.)\n", "\n", "With this model in hand, continuum-subtracted or continuum-normalized spectra can be produced using basic spectral manipulations:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "jwst_gencont_sub = jwst_spec - generic_continuum(jwst_spec.spectral_axis)\n", "jwst_gencont_norm = jwst_spec / generic_continuum(jwst_spec.spectral_axis)\n", "\n", "ax1, ax2 = plt.subplots(2, 1)[1]\n", "\n", "ax1.step(jwst_gencont_sub.wavelength, jwst_gencont_sub.flux)\n", "ax1.set_ylim(-50, 50)\n", "ax1.axhline(0, color='k', ls=':') # continuum should be at flux=0\n", "\n", "ax2.step(jwst_gencont_norm.wavelength, jwst_gencont_norm.flux)\n", "ax2.set_ylim(0, 2)\n", "ax2.axhline(1, color='k', ls='--'); # continuum should be at flux=1" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### The customizable way\n", "\n", "The `fit_continuum` function operates similarly to `fit_generic_continuum`, but is meant for you to provide your favorite continuum model rather than being tailored to a specific continuum model. To see the list of models, see the [astropy.modeling documentation](http://docs.astropy.org/en/stable/modeling/index.html)." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "from specutils.fitting import fit_continuum\n", "from astropy.modeling import models" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "For example, suppose you want to use a 3rd-degree Chebyshev polynomial as your continuum model. You can use `fit_continuum` to get an object that behaves the same as for `fit_generic_continuum`:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "chebdeg3_continuum = fit_continuum(jwst_spec, models.Chebyshev1D(3))\n", "\n", "generic_continuum_evaluated = generic_continuum(jwst_spec.spectral_axis)\n", "\n", "plt.step(jwst_spec.spectral_axis, jwst_spec.flux)\n", "plt.plot(jwst_spec.spectral_axis, chebdeg3_continuum(jwst_spec.spectral_axis))\n", "plt.ylim(0.4, 1.1);" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This then provides total flexibility. For example, you can also try other polynomials like higher-degree Hermite polynomials:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "hermdeg7_continuum = fit_continuum(jwst_spec, models.Hermite1D(degree=7))\n", "hermdeg17_continuum = fit_continuum(jwst_spec, models.Hermite1D(degree=17))\n", "\n", "plt.step(jwst_spec.spectral_axis, jwst_spec.flux)\n", "plt.plot(jwst_spec.spectral_axis, chebdeg3_continuum(jwst_spec.spectral_axis))\n", "plt.plot(jwst_spec.spectral_axis, hermdeg7_continuum(jwst_spec.spectral_axis))\n", "plt.plot(jwst_spec.spectral_axis, hermdeg17_continuum(jwst_spec.spectral_axis))\n", "plt.ylim(0.4, 1.1);" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This immediately demonstrates the tradeoffs in polynomial fitting: while the high-degree polynomials capture the wiggles of the spectrum better than the low, they also *over*-fit near the strong emission lines." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Exercises\n", "\n", "Try combining the `SpectralRegion` and continuum-fitting functionality to only fit the parts of the spectrum that *are* continuum (i.e. not including emission lines). Can you do better?" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Using the spectrum from the previous exercise, first subtract a continuum, then re-do your measurement. Is it better?" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Line-Fitting\n", "\n", "In addition to the more empirical measurements described above, `specutils` provides tools for doing spectral line fitting. The approach is akin to that for continuum modeling: models from [astropy.modeling](http://docs.astropy.org/en/stable/modeling/index.html) are fit to the spectrum, and either those models can be used directly, or their parameters." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The fitting machinery must first be given guesses for line locations. This process can be automated using functions designed to identify lines (more detail on the options is [in the docs](https://specutils.readthedocs.io/en/latest/fitting.html#line-finding)). For data sets where these algorithms are not ideal, you may substitute your own (i.e., skip this step and start with line location guesses). \n", "\n", "Here we identify the three lines near the Halpha region in our spectrum, finding the lines above about a $\\sim 5 \\sigma$ flux threshold. They are then output as an astropy Table:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "jwst_ha_contsub = jwst_spec_with_unc_ha - jwst_continuum\n", "\n", "halpha_lines = fitting.find_lines_threshold(jwst_ha_contsub, 5)\n", "plt.step(jwst_spec_with_unc_ha.spectral_axis, jwst_spec_with_unc_ha.flux, where='mid')\n", "for line in halpha_lines:\n", " plt.axvline(line['line_center'], color='k', ls=':')\n", "\n", "halpha_lines" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now for each of these lines, we need to fit a model. Sometimes it is sufficient to simply create a model where the center is at the line and excise the appropriate area of the line to do a line estimate. This is not *too* sensitive to the size of the region, at least for well-separated lines like these. The result is a list of models that carry with them them the details of the fit:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "halpha_line_models = []\n", "for line in halpha_lines:\n", " line_region = SpectralRegion(line['line_center']-50*u.angstrom,\n", " line['line_center']+50*u.angstrom)\n", " line_spectrum = manipulation.extract_region(jwst_ha_contsub, line_region)\n", " line_spectrum\n", " \n", " # here's the workaround from above again\n", " line_spectrum = Spectrum1D(flux=line_spectrum.flux, \n", " spectral_axis=line_spectrum.spectral_axis, \n", " uncertainty=line_spectrum.uncertainty)\n", " line_estimate = models.Gaussian1D(mean=line['line_center'],stddev=.001*u.um)\n", " line_model = fitting.fit_lines(line_spectrum, line_estimate)\n", " \n", " halpha_line_models.append(line_model)\n", " \n", "plt.step(jwst_ha_contsub.spectral_axis, jwst_ha_contsub.flux, where='mid')\n", "for line_model in halpha_line_models:\n", " evaluated_model = line_model(jwst_ha_contsub.spectral_axis)\n", " plt.plot(jwst_ha_contsub.spectral_axis, evaluated_model) \n", " \n", "halpha_line_models" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "For more complicated models or fits it may be better to use the `estimate_line_parameters` function instead of manually creating e.g. a `Gaussian1D` model and setting the center. An example of this pattern is given below.\n", "\n", "Note that we provided a default `Gaussian1D` model to the `estimate_line_parameters` function above. This function makes reasonable guesses for `Gaussian1D`, `Voigt1D`, and `Lorentz1D`, the most common line profiles used for spectral lines, but may or may not work for other models. See [the relevant docs section](https://specutils.readthedocs.io/en/latest/fitting.html#parameter-estimation) for more details.\n", "\n", "In this example we also show an example of a *joint* fit of both lines at the same time. While the difference may seems subtle, in cases of blended lines this typically provides much better fits:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "halpha_line_estimates = []\n", "for line in halpha_lines:\n", " line_region = SpectralRegion(line['line_center']-15*u.angstrom,\n", " line['line_center']+15*u.angstrom)\n", " line_spectrum = manipulation.extract_region(jwst_ha_contsub, line_region)\n", " line_estimate = fitting.estimate_line_parameters(line_spectrum, models.Gaussian1D())\n", " \n", " halpha_line_estimates.append(line_estimate)\n", "\n", "# this could be done more flexibly with a for loop but we are explicit here for simplicity\n", "combined_model_estimate = halpha_line_estimates[0] + halpha_line_estimates[1]\n", "\n", "plt.step(jwst_ha_contsub.spectral_axis, jwst_ha_contsub.flux, where='mid')\n", "plt.plot(jwst_ha_contsub.spectral_axis, \n", " combined_model_estimate(jwst_ha_contsub.spectral_axis)) \n", "\n", "combined_model_estimate" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "While these estimates are already pretty decent, to complete the story (especially for blended lines) the last step is to then fit the combined model:" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "combined_model = fitting.fit_lines(jwst_ha_contsub, combined_model_estimate)\n", "\n", "plt.step(jwst_ha_contsub.spectral_axis, jwst_ha_contsub.flux, where='mid')\n", "plt.plot(jwst_ha_contsub.spectral_axis, \n", " combined_model(jwst_ha_contsub.spectral_axis)) \n", " \n", "combined_model" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "### Exercise\n", "\n", "Fit a spectral feature from your own spectrum using the fitting methods outlined above. Try the different line profile types (Gaussian, Lorentzian, or Voigt). If you are using the blackbody spectrum (where you know the \"true\" answer for the spectral line), compare your answer to the true answer." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Additional Manipulations\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(jwst_spec, 0.75) # 0.75 pixel gaussian kernel\n", "plt.step(smoothed_spec.wavelength, smoothed_spec.flux)" ] }, { "cell_type": "markdown", "metadata": { "tags": [] }, "source": [ "### Exercise\n", "\n", "Try smoothing your spectrum loaded in the first exercise, or the JWST spectrum. Compare all available kernel types (see the `convolution_smooth()` function) and decide which seems most appropriate for your spectrum." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Additional Exercises\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": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "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": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Now, try to make line measurements as described above on your model spectral feature. Try varying the amplitude and see how well you recover the line's properties as it sinks down into the noise." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "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.9.13" } }, "nbformat": 4, "nbformat_minor": 4 }