From 7acc36d431845062a6451d5a7f0cd14946b51173 Mon Sep 17 00:00:00 2001 From: Erik Tollerud Date: Tue, 11 May 2021 18:55:20 -0400 Subject: [PATCH] Merge pull request #23 from eteq/specutils-dry-run Specutils based on dry run feedback --- jdat_session/Specutils_Webbinar_live.ipynb | 344 ++++++++++----- .../Specutils_Webbinar_solutions.ipynb | 400 ++++++++++++------ 2 files changed, 506 insertions(+), 238 deletions(-) diff --git a/jdat_session/Specutils_Webbinar_live.ipynb b/jdat_session/Specutils_Webbinar_live.ipynb index 747ade8..a64c7e1 100644 --- a/jdat_session/Specutils_Webbinar_live.ipynb +++ b/jdat_session/Specutils_Webbinar_live.ipynb @@ -32,9 +32,10 @@ "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\n", + "from specutils import Spectrum1D, SpectralRegion, analysis, manipulation, fitting\n", "specutils.__version__" ] }, @@ -94,14 +95,29 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "## Basics of creating Spectrum1D Objects" + "## 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": [ - "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." + "## 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": [ + "from astropy.utils.data import download_file\n", + "\n", + "spec_fn = download_file('https://stsci.box.com/shared/static/b22b1fzhimtdqfp8597m4bg67kovvauu.fits', cache=True)" ] }, { @@ -115,9 +131,7 @@ "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:" + "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:" ] }, { @@ -125,13 +139,17 @@ "execution_count": null, "metadata": {}, "outputs": [], - "source": [] + "source": [ + "help(Spectrum1D.read)" + ] }, { "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:" + "### 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, or [an sdss galaxy](https://dr14.sdss.org/optical/spectrum/view/data/format=fits/spec=lite?plateid=1323&mjd=52797&fiberid=12)), and load it." ] }, { @@ -145,9 +163,14 @@ "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):" + "# 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: " ] }, { @@ -161,11 +184,9 @@ "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", + "## Working with Units and Spectral Axes \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)." + "We created `Spectrum1D` just as Quantity arrays, so they can be treated just as `Quantity` objects when convenient with unit conversions and the like:" ] }, { @@ -179,9 +200,7 @@ "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 simulated high-redshift (z > 1) galaxy like that you might see from NIRSpec:" + "There's even fast-accessors to make some of this more convenient:" ] }, { @@ -195,9 +214,7 @@ "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, or [an sdss galaxy](https://dr14.sdss.org/optical/spectrum/view/data/format=fits/spec=lite?plateid=1323&mjd=52797&fiberid=12)), and load it." + "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:" ] }, { @@ -211,11 +228,16 @@ "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.:" + "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" ] }, { @@ -229,18 +251,20 @@ "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:" + "But specutils makes this much easier!" ] }, { "cell_type": "code", "execution_count": null, - "metadata": { - "scrolled": true, - "tags": [ - "raises-exception" - ] - }, + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, "outputs": [], "source": [] }, @@ -248,9 +272,9 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "### Exercise\n", + "## Uncertainties\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.)" + "Currently the most compatible way to use uncertainties is the machinery built for the `astropy.nddata` object (although in some cases simply passing in an uncertainty array will also work):" ] }, { @@ -264,9 +288,11 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "## Smoothing, etc.\n", + "# Arithmetic on Spectra\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:" + "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.:" ] }, { @@ -280,9 +306,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "### Exercise\n", - "\n", - "Try smoothing your spectrum loaded in the above example (or the JWST spectrum). Compare all available kernel types and decide which seems most appropriate for your spectrum." + "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:" ] }, { @@ -296,9 +320,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "# Making Measurements on specutils Spectra\n", - "\n", - "The remainder of this Webbinar foucuses on analysis steps that lead to measurements of spectra - e.g. line flux, line width, model-fitting, and parameter extraction. Most of this functionality is in `specutils.analysis` and `specutils.manipulation`:" + "Resolving this requires explicit conversion:" ] }, { @@ -306,17 +328,33 @@ "execution_count": null, "metadata": {}, "outputs": [], + "source": [] + }, + { + "cell_type": "markdown", + "metadata": {}, "source": [ - "from specutils import analysis, manipulation" + "# 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": [] + }, { "cell_type": "markdown", "metadata": {}, "source": [ - "## Sample Spectrum and SNR\n", + "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", - "As a reminder, here is the simulated JWST spectrum we loaded above:" + "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$" ] }, { @@ -325,15 +363,21 @@ "metadata": {}, "outputs": [], "source": [ - "plt.step(jwst_spec.spectral_axis, jwst_spec.flux)\n", - "plt.ylim(0, plt.ylim()[-1])" + "halpha_lines_region = SpectralRegion(16380*u.angstrom, 16440*u.angstrom)" ] }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + }, { "cell_type": "markdown", "metadata": {}, "source": [ - "Lets start with 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 straightofrward algorithm detailed in the literature." + "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:" ] }, { @@ -347,7 +391,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "Ideally we'd use an uncertainty *we* understand to derive the S/N:" + "Now we simply plug this into the redshift formula for H$\\alpha$'s rest wavelength:" ] }, { @@ -361,26 +405,25 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "But we need to define the uncertainty first! To do that we need to specify the region to compute the uncertainty:" + "We now have strong evidence this is a manufactured spectrum! (it's suspiciously close to *exactly* 1.5 ...)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "# Spectral Regions\n", + "# Additional Materials\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." + "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": "code", - "execution_count": null, + "cell_type": "markdown", "metadata": {}, - "outputs": [], "source": [ - "from specutils import SpectralRegion\n", - "from specutils.manipulation import extract_region" + "## Sample Spectrum and SNR\n", + "\n", + "Lets start with 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 straightofrward algorithm detailed in the literature." ] }, { @@ -388,16 +431,13 @@ "execution_count": null, "metadata": {}, "outputs": [], - "source": [ - "ha_region = SpectralRegion((6563-50)*u.AA, (6563+50)*u.AA)\n", - "ha_region" - ] + "source": [] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "Regions can also be raw pixel values (although of course this is more applicable to a specific spectrum):" + "Ideally we'd use an uncertainty *we* understand to derive the S/N:" ] }, { @@ -405,16 +445,22 @@ "execution_count": null, "metadata": {}, "outputs": [], + "source": [] + }, + { + "cell_type": "markdown", + "metadata": {}, "source": [ - "ha_pixel_region = SpectralRegion(2650*u.pixel, 2850*u.pixel)\n", - "ha_pixel_region" + "But we need to define the uncertainty first! To do that we need to specify the region to compute the uncertainty:" ] }, { "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:" + "## 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." ] }, { @@ -422,16 +468,13 @@ "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" - ] + "source": [] }, { "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:" + "Regions can also be raw pixel values (although of course this is more applicable to a specific spectrum):" ] }, { @@ -445,7 +488,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "With that uncertainty present, we can now compute the S/N:" + "Additionally, *multiple* regions can be in the same `SpectralRegion` object. This is useful for e.g. measuring multiple spectral features in one call:" ] }, { @@ -459,7 +502,7 @@ "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):" + "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:" ] }, { @@ -473,16 +516,21 @@ "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." + "With that uncertainty present, we can now compute the S/N:" ] }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + }, { "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." + "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):" ] }, { @@ -496,7 +544,16 @@ "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. " + "## 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." ] }, { @@ -510,7 +567,9 @@ "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." + "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." ] }, { @@ -524,7 +583,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "Note this simple centroid can also give us a reasonable estimate for the redshift:" + "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." ] }, { @@ -622,7 +681,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "# Continuum Subtraction\n", + "## 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." ] @@ -641,8 +700,15 @@ "execution_count": null, "metadata": {}, "outputs": [], + "source": [] + }, + { + "cell_type": "markdown", + "metadata": {}, "source": [ - "from specutils.fitting import fit_generic_continuum" + "(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:" ] }, { @@ -656,9 +722,9 @@ "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", + "### The customizable way\n", "\n", - "With this model in hand, continuum-subtracted or continuum-normalized spectra can be produced using basic spectral manipulations:" + "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)." ] }, { @@ -672,9 +738,7 @@ "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)." + "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`:" ] }, { @@ -682,16 +746,13 @@ "execution_count": null, "metadata": {}, "outputs": [], - "source": [ - "from specutils.fitting import fit_continuum\n", - "from astropy.modeling import models" - ] + "source": [] }, { "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`:" + "This then provides total flexibility. For example, you can also try other polynomials like higher-degree Hermite polynomials:" ] }, { @@ -705,7 +766,14 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "## Exercise\n", + "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?" ] @@ -721,8 +789,6 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "## Exercise\n", - "\n", "Using the spectrum from the previous exercise, first subtract a continuum, then re-do your measurement. Is it better?" ] }, @@ -737,7 +803,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "# Line-Fitting\n", + "## 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." ] @@ -747,9 +813,7 @@ "execution_count": null, "metadata": {}, "outputs": [], - "source": [ - "from specutils import fitting" - ] + "source": [] }, { "cell_type": "markdown", @@ -817,11 +881,91 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "## Exercise\n", + "### 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": [] + }, + { + "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, diff --git a/jdat_session/Specutils_Webbinar_solutions.ipynb b/jdat_session/Specutils_Webbinar_solutions.ipynb index 34d406d..5389b82 100644 --- a/jdat_session/Specutils_Webbinar_solutions.ipynb +++ b/jdat_session/Specutils_Webbinar_solutions.ipynb @@ -32,9 +32,10 @@ "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\n", + "from specutils import Spectrum1D, SpectralRegion, analysis, manipulation, fitting\n", "specutils.__version__" ] }, @@ -94,14 +95,18 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "## Basics of creating Spectrum1D Objects" + "## 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": [ - "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.:" + "## 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:" ] }, { @@ -110,20 +115,26 @@ "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", + "from astropy.utils.data import download_file\n", "\n", - "spec1d" + "spec_fn = download_file('https://stsci.box.com/shared/static/b22b1fzhimtdqfp8597m4bg67kovvauu.fits', 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": [ - "These can be straightforwardly transformed to other spectral units if desired:" + "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:" ] }, { @@ -132,17 +143,16 @@ "metadata": {}, "outputs": [], "source": [ - "jyspec1d = spec1d.new_flux_unit(u.Jy)\n", - "plt.plot(jyspec1d.frequency, jyspec1d.flux)" + "help(Spectrum1D.read)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "## Working with Spectral Axes \n", + "### Exercise\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:" + "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, or [an sdss galaxy](https://dr14.sdss.org/optical/spectrum/view/data/format=fits/spec=lite?plateid=1323&mjd=52797&fiberid=12)), and load it." ] }, { @@ -150,8 +160,20 @@ "execution_count": null, "metadata": {}, "outputs": [], + "source": [] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Creating a Spectrum \"by-hand\"" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, "source": [ - "spec1d.spectral_axis" + "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: " ] }, { @@ -160,7 +182,8 @@ "metadata": {}, "outputs": [], "source": [ - "spec1d.spectral_axis.to(u.nm)" + "fitsfile = fits.open(spec_fn)\n", + "fitsfile.info()" ] }, { @@ -169,14 +192,17 @@ "metadata": {}, "outputs": [], "source": [ - "spec1d.spectral_axis.to(u.THz, u.spectral())" + " fitsfile[1].header" ] }, { - "cell_type": "markdown", + "cell_type": "code", + "execution_count": null, "metadata": {}, + "outputs": [], "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:" + "table_data = fitsfile[1].data\n", + "table_data" ] }, { @@ -185,7 +211,23 @@ "metadata": {}, "outputs": [], "source": [ - "spec1d.wcs.pixel_to_world([10, 10.5])" + "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": [ + "## 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:" ] }, { @@ -194,16 +236,16 @@ "metadata": {}, "outputs": [], "source": [ - "spec1d.wcs.world_to_pixel([400, 450]*u.nm)" + "jwst_spec.spectral_axis" ] }, { - "cell_type": "markdown", + "cell_type": "code", + "execution_count": null, "metadata": {}, + "outputs": [], "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):" + "jwst_spec.spectral_axis.to(u.angstrom)" ] }, { @@ -212,29 +254,14 @@ "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')" + "jwst_spec.spectral_axis.to(u.THz, u.spectral())" ] }, { "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)." + "There's even fast-accessors to make some of this more convenient:" ] }, { @@ -242,15 +269,15 @@ "execution_count": null, "metadata": {}, "outputs": [], - "source": [] + "source": [ + "plt.step(jwst_spec.frequency, jwst_spec.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 simulated high-redshift (z > 1) galaxy like that you might see from NIRSpec:" + "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:" ] }, { @@ -259,9 +286,7 @@ "metadata": {}, "outputs": [], "source": [ - "from astropy.utils.data import download_file\n", - "\n", - "spec_fn = download_file('https://stsci.box.com/shared/static/b22b1fzhimtdqfp8597m4bg67kovvauu.fits', cache=True)" + "jwst_spec.wcs.pixel_to_world([10, 10.5])" ] }, { @@ -270,15 +295,14 @@ "metadata": {}, "outputs": [], "source": [ - "jwst_spec = Spectrum1D.read(spec_fn)\n", - "plt.plot(jwst_spec.spectral_axis, jwst_spec.flux);" + "jwst_spec.wcs.world_to_pixel([1.2, 1.3]*u.um)" ] }, { "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:" + "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:" ] }, { @@ -287,16 +311,44 @@ "metadata": {}, "outputs": [], "source": [ - "help(Spectrum1D.read)" + "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": [ - "### Exercise\n", + "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": [ + "## Uncertainties\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, or [an sdss galaxy](https://dr14.sdss.org/optical/spectrum/view/data/format=fits/spec=lite?plateid=1323&mjd=52797&fiberid=12)), and load it." + "Currently the most compatible way to use uncertainties is the machinery built for the `astropy.nddata` object (although in some cases simply passing in an uncertainty array will also work):" ] }, { @@ -304,17 +356,28 @@ "execution_count": null, "metadata": {}, "outputs": [], - "source": [] + "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": [ - "# Manipulating spectra\n", + "# Arithmetic on 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", + "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 expected patterns. E.g.:" + "The most straightforward of operations are arithmetic manipulations. In general these follow patterns that are based on fundametal arithmetic. E.g.:" ] }, { @@ -323,7 +386,7 @@ "metadata": {}, "outputs": [], "source": [ - "newspec1d = spec1d - 2 * u.erg / u.s / (u.cm**2) / u.AA\n", + "newspec1d = jwst_spec - 0.5*u.uJy\n", "plt.step(newspec1d.wavelength, newspec1d.flux)" ] }, @@ -346,7 +409,7 @@ "outputs": [], "source": [ "# This raises an error:\n", - "newspec1d - jyspec1d " + "newspec1d - jwst_spec_flamb " ] }, { @@ -362,16 +425,18 @@ "metadata": {}, "outputs": [], "source": [ - "spec1d.new_flux_unit(u.Jy) - jyspec1d" + "newspec1d.new_flux_unit(f_lamb_units) - jwst_spec_flamb" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "### Exercise\n", + "# Sample Analysis: Line Center for Redshift\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.)" + "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." ] }, { @@ -379,15 +444,23 @@ "execution_count": null, "metadata": {}, "outputs": [], - "source": [] + "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": [ - "## Smoothing, etc.\n", + "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", - "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:" + "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$" ] }, { @@ -396,19 +469,22 @@ "metadata": {}, "outputs": [], "source": [ - "from specutils.manipulation import gaussian_smooth\n", + "halpha_lines_region = SpectralRegion(16380*u.angstrom, 16440*u.angstrom)\n", "\n", - "smoothed_spec = gaussian_smooth(spec1d, 0.75) # 0.75 pixel gaussian kernel\n", - "plt.step(smoothed_spec.wavelength, smoothed_spec.flux)" + "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": [ - "### Exercise\n", - "\n", - "Try smoothing your spectrum loaded in the above example (or the JWST spectrum). Compare all available kernel types and decide which seems most appropriate for your spectrum." + "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:" ] }, { @@ -416,15 +492,16 @@ "execution_count": null, "metadata": {}, "outputs": [], - "source": [] + "source": [ + "cen = analysis.centroid(jwst_halpha_contsub, halpha_lines_region)\n", + "cen" + ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "# Making Measurements on specutils Spectra\n", - "\n", - "The remainder of this Webbinar foucuses on analysis steps that lead to measurements of spectra - e.g. line flux, line width, model-fitting, and parameter extraction. Most of this functionality is in `specutils.analysis` and `specutils.manipulation`:" + "Now we simply plug this into the redshift formula for H$\\alpha$'s rest wavelength:" ] }, { @@ -433,32 +510,31 @@ "metadata": {}, "outputs": [], "source": [ - "from specutils import analysis, manipulation" + "(cen/(6563*u.angstrom)) - 1" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "## Sample Spectrum and SNR\n", - "\n", - "As a reminder, here is the simulated JWST spectrum we loaded above:" + "We now have strong evidence this is a manufactured spectrum! (it's suspiciously close to *exactly* 1.5 ...)" ] }, { - "cell_type": "code", - "execution_count": null, + "cell_type": "markdown", "metadata": {}, - "outputs": [], "source": [ - "plt.step(jwst_spec.spectral_axis, jwst_spec.flux)\n", - "plt.ylim(0, plt.ylim()[-1])" + "# 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 Spectrum and SNR\n", + "\n", "Lets start with 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 straightofrward algorithm detailed in the literature." ] }, @@ -498,21 +574,11 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "# Spectral Regions\n", + "## 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": [ - "from specutils import SpectralRegion\n", - "from specutils.manipulation import extract_region" - ] - }, { "cell_type": "code", "execution_count": null, @@ -615,7 +681,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "# Line Measurements\n", + "## 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." ] @@ -634,13 +700,12 @@ "outputs": [], "source": [ "# estimate a reasonable continuum-level estimate for the h-alpha area of the spectrum\n", - "jwst_continuum = 560*u.nJy # note this is more convenient units than the spectrum itself\n", + "jwst_continuum = 550*u.nJy # note this is more convenient units than the spectrum itself\n", "\n", - "jwst_halpha_contsub = extract_region(jwst_spec, ha_pixel_region) - jwst_continuum\n", + "jwst_halpha_contsub = manipulation.extract_region(jwst_spec, ha_pixel_region) - jwst_continuum\n", "\n", "plt.axhline(0, c='k', ls=':')\n", - "plt.step(jwst_halpha_contsub.spectral_axis, jwst_halpha_contsub.flux)\n", - "#plt.ylim(-50, 50)" + "plt.step(jwst_halpha_contsub.spectral_axis, jwst_halpha_contsub.flux)" ] }, { @@ -658,8 +723,8 @@ "metadata": {}, "outputs": [], "source": [ - "LOWER = 16380 * u.angstrom\n", - "UPPER = 16440 * u.angstrom\n", + "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", @@ -667,7 +732,8 @@ "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)" + "plt.ylim(yl1, yl2)\n", + "plt.xlim(1.6, 1.7);" ] }, { @@ -687,22 +753,6 @@ "cen" ] }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Note this simple centroid can also give us a reasonable estimate for the redshift:" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "(cen/(6563*u.angstrom)) - 1" - ] - }, { "cell_type": "markdown", "metadata": {}, @@ -800,7 +850,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "# Continuum Subtraction\n", + "## 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." ] @@ -942,7 +992,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "## Exercise\n", + "### 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?" ] @@ -958,8 +1008,6 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "## Exercise\n", - "\n", "Using the spectrum from the previous exercise, first subtract a continuum, then re-do your measurement. Is it better?" ] }, @@ -974,20 +1022,11 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "# Line-Fitting\n", + "## 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": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "from specutils import fitting" - ] - }, { "cell_type": "markdown", "metadata": {}, @@ -1030,7 +1069,7 @@ "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 = extract_region(jwst_ha_contsub, line_region)\n", + " line_spectrum = manipulation.extract_region(jwst_ha_contsub, line_region)\n", " line_spectrum\n", " \n", " # here's the workaround from above again\n", @@ -1071,7 +1110,7 @@ "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 = extract_region(jwst_ha_contsub, line_region)\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", @@ -1112,11 +1151,96 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "## Exercise\n", + "### 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,