{
"cells": [
{
"cell_type": "markdown",
"id": "e6c242e1",
"metadata": {},
"source": [
"# JWebbinar 8: Redshift and Template Fitting\n",
"\n",
"This notebook covers basic examples on how a user can measure the redshift of a source using the visualization tool [Jdaviz](https://jdaviz.readthedocs.io/en/latest/) or programmatically with [Specutils](https://specutils.readthedocs.io/en/latest/).\n",
"\n",
"**Content**\n",
"- [Resources and documentation](#resources)\n",
"- [Installation](#installation)\n",
"- [Imports](#imports)\n",
"- [Fetch the example data and the mystery data](#data)\n",
"- [\"By eye\" redshift measurement with Specviz](#byeye)\n",
" - Exercise: measure the redshift on the mystery data\n",
"- [Redshift measurement with cross-correlation method](#crosscorr)\n",
" - Get a template and prepare it for use\n",
" - Subtract the continuum from the observed spectrum\n",
" - Exercise: fit the continuum using Specviz\n",
" - Run the cross correlation function\n",
"- [Template fitting using a library of spectra](#fitting)\n",
" - Get the library of spectra\n",
" - Run the template_match function\n",
"\n",
"**Author**: Camilla Pacifici (cpacifici@stsci.edu)\n",
"\n",
"**Review**: Kyle Conroy (kconroy@stsci.edu)"
]
},
{
"cell_type": "markdown",
"id": "9ae98afa",
"metadata": {},
"source": [
"\n",
"## Resources and documentation\n",
"\n",
"This notebook uses functionality from [Specutils](https://specutils.readthedocs.io/en/stable/) and [Jdaviz](https://jdaviz.readthedocs.io/en/latest/). Developers at the Space Telescope Science Institute are available to answer questions and resolve problems through the [JWST Help Desk](https://jwsthelp.stsci.edu/). If you wish to send feedback or report problems, you can also submit an issue directly on Github, both for [Specutils](https://github.com/astropy/specutils) and for [Jdaviz](https://github.com/spacetelescope/jdaviz)."
]
},
{
"cell_type": "markdown",
"id": "95a0cf36",
"metadata": {},
"source": [
"\n",
"## Installation\n",
"\n",
"If you want to create an environment for all the notebooks presented during JWebbinar 8, please follow the instructions in the README to use the requirement file. If you want to run this notebook only, you can install just the jdaviz package with the following instructions. \n",
"`conda create -n jdaviz python=3.8` \n",
"`conda activate jdaviz` \n",
"from the latest release \n",
"`pip install jdaviz` \n",
"or from git \n",
"`pip install git+https://github.com/spacetelescope/jdaviz.git`"
]
},
{
"cell_type": "markdown",
"id": "1f2020d3",
"metadata": {},
"source": [
"\n",
"## Imports"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "45a3ae5e",
"metadata": {},
"outputs": [],
"source": [
"# This ensures that our notebook is using the full width of the browser\n",
"from IPython.core.display import display, HTML\n",
"display(HTML(\"\"))\n",
"\n",
"# general os\n",
"import os\n",
"import zipfile\n",
"import urllib.request\n",
"\n",
"# numpy\n",
"import numpy as np\n",
"\n",
"# specviz\n",
"import jdaviz # this is needed to get the version number later\n",
"from jdaviz import Specviz\n",
"\n",
"# astropy\n",
"import astropy # again for the version number\n",
"import astropy.units as u\n",
"from astropy.io import fits, ascii\n",
"from astropy.utils.data import download_file\n",
"from astropy.modeling.models import Linear1D, Polynomial1D, Chebyshev1D\n",
"from astropy.nddata import StdDevUncertainty\n",
"from astropy.visualization import quantity_support\n",
"\n",
"# specutils\n",
"import specutils # again for the version number\n",
"from specutils import Spectrum1D, SpectralRegion\n",
"from specutils.fitting import fit_generic_continuum\n",
"from specutils.analysis import correlation\n",
"from specutils.analysis import template_comparison\n",
"from specutils.manipulation import LinearInterpolatedResampler\n",
"from specutils.manipulation import extract_region\n",
"\n",
"\n",
"# matplotlib (in case we need it)\n",
"from matplotlib import pyplot as plt\n",
"# customization of matplotlib style\n",
"plt.rcParams[\"figure.figsize\"] = (10, 5)\n",
"params = {'legend.fontsize': '18', 'axes.labelsize': '18',\n",
" 'axes.titlesize': '18', 'xtick.labelsize': '18',\n",
" 'ytick.labelsize': '18', 'lines.linewidth': 2,\n",
" 'axes.linewidth': 2, 'animation.html': 'html5',\n",
" 'figure.figsize': (8, 6)}\n",
"plt.rcParams.update(params)\n",
"plt.rcParams.update({'figure.max_open_warning': 0})"
]
},
{
"cell_type": "markdown",
"id": "1fdb9a16",
"metadata": {},
"source": [
"### Versions:"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "c36d248d",
"metadata": {},
"outputs": [],
"source": [
"print(\"jdaviz:\", jdaviz.__version__)\n",
"print(\"astropy:\", astropy.__version__)\n",
"print(\"specutils:\", specutils.__version__)"
]
},
{
"cell_type": "markdown",
"id": "16081eb5",
"metadata": {},
"source": [
"\n",
"## Fetch the example data and the mystery data\n",
"\n",
"Here we download a spectrum extracted from a NIRSpec IFU simulation, a mystery dataset you will explore later in an exercise, and a model spectrum we will use as template for the redshift measurement. The template is based on a combination of Simple Stellar Population models including emission lines as done in [Pacifici et al. (2012)](https://ui.adsabs.harvard.edu/abs/2012MNRAS.421.2002P/abstract)."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "55d54dd2",
"metadata": {},
"outputs": [],
"source": [
"fn_nirspec = download_file('https://stsci.box.com/shared/static/llbr4w8ucqpz6g81hi41v4tnszwi3bq4.fits', cache=True)\n",
"fn_mystery = download_file('https://stsci.box.com/shared/static/h51qakmow49j3isdlxeln4fcxrw7fudc.fits', cache=True)\n",
"fn_template = download_file('https://stsci.box.com/shared/static/3rkurzwl0l79j70ddemxafhpln7ljle7.dat', cache=True)"
]
},
{
"cell_type": "markdown",
"id": "68fa01fc",
"metadata": {},
"source": [
"\n",
"## \"By eye\" redshift measurement with Specviz\n",
"\n",
"Specviz will allow you to match line wavelengths to the emission lines you see in the spectrum. You will be able to do this using the [redshift slider](https://jdaviz.readthedocs.io/en/latest/specviz/redshift.html) at the top of the interface. But first, let us [open the spectrum in Specviz](https://jdaviz.readthedocs.io/en/latest/specviz/import_data.html)."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "2dbeb862",
"metadata": {},
"outputs": [],
"source": [
"# Call the app\n",
"viz = Specviz()\n",
"viz.app"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "51224ebc",
"metadata": {},
"outputs": [],
"source": [
"# Load spectrum\n",
"viz.load_spectrum(fn_nirspec, data_label=\"NIRSpec\")"
]
},
{
"cell_type": "markdown",
"id": "8b96dfda",
"metadata": {},
"source": [
"Now we need to:\n",
"- open the [\"line list\" plugin](https://jdaviz.readthedocs.io/en/latest/specviz/plugins.html?highlight=plugin#line-lists)\n",
"- choose pre-loaded lines or add custom lines (the lines will not show in the viewer because they are plotted at restframe)\n",
"- input a guess redshift into the redshift box\n",
"- move the slider to get a better match\n",
"- use the [zoom tool](https://jdaviz.readthedocs.io/en/latest/specviz/displaying.html#pan-zoom) to get an even better match"
]
},
{
"cell_type": "markdown",
"id": "41f7eb34",
"metadata": {},
"source": [
"### Exercise: measure the redshift on the mystery data\n",
"\n",
"You can open a new instance of Specviz and re-do all of the above.\n",
"\n",
"*Hint:* \n",
"*viz2 = SpecViz()* \n",
"*viz2.app*"
]
},
{
"cell_type": "markdown",
"id": "f9da086b",
"metadata": {},
"source": [
"The mystery dataset is a JWST simulation created with [MIRAGE](https://mirage-data-simulator.readthedocs.io/en/latest/) for the NIRISS wide field slitless mode and run through the [JWST Calibration Pipeline](https://jwst-docs.stsci.edu/jwst-data-reduction-pipeline). It is a combination of three filters (F115W, F150W, and F200W)."
]
},
{
"cell_type": "markdown",
"id": "79a5d1b3",
"metadata": {},
"source": [
"\n",
"## Redshift measurement with cross-correlation method\n",
"\n",
"It is very common in astronomy to measure a redshift using a cross-correlation algorithm. IRAF uses this methodology in its [xcsao](http://tdc-www.harvard.edu/iraf/rvsao/xcsao/xcsao.html) task. Here, we use the Specutils [template cross-correlation](https://specutils.readthedocs.io/en/stable/analysis.html#template-cross-correlation) function to derive the redshift of our source. There are a couple of things that we need to do before we run the correlation algorithm:\n",
"- get a template spectrum for the correlation\n",
"- subtract the continuum from both the template and the observed spectrum\n",
"- make sure the spectra have some overlap in wavelength"
]
},
{
"cell_type": "markdown",
"id": "53f4ab27",
"metadata": {},
"source": [
"### Get a template and prepare it for use"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "52a6979a",
"metadata": {},
"outputs": [],
"source": [
"# Define unit\n",
"spec_unit = u.erg / (u.s * u.cm**2 * u.AA)\n",
"\n",
"# Read spectrum with the ascii function\n",
"template = ascii.read(fn_template)\n",
"# Create Spectrum1D object\n",
"template = Spectrum1D(spectral_axis=template['col1']/1E4*u.um, \n",
" flux=template['col2']*spec_unit)"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "8348a12e",
"metadata": {},
"outputs": [],
"source": [
"# Cut to useful range - template and obs MUST overlap, so we go to 1.1um\n",
"use_tmp = (template.spectral_axis.value > 0.3) & (template.spectral_axis.value < 1.1)\n",
"template_cut = Spectrum1D(spectral_axis=template.spectral_axis[use_tmp], flux=template.flux[use_tmp])"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "ae483af5",
"metadata": {},
"outputs": [],
"source": [
"# Look at spectrum\n",
"plt.figure(figsize=[10, 6])\n",
"plt.plot(template_cut.spectral_axis, template_cut.flux)\n",
"plt.xlabel(\"wavelength ({:latex})\".format(template_cut.spectral_axis.unit))\n",
"plt.ylabel(\"flux ({:latex})\".format(template_cut.flux.unit))\n",
"plt.title(\"Template\")\n",
"plt.show()"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "a1488f36",
"metadata": {},
"outputs": [],
"source": [
"# Subtract continuum\n",
"mask_temp = ((template_cut.spectral_axis.value > 0.31) & (template_cut.spectral_axis.value < 0.37) |\n",
" (template_cut.spectral_axis.value > 0.40) & (template_cut.spectral_axis.value < 0.47) |\n",
" (template_cut.spectral_axis.value > 0.52) & (template_cut.spectral_axis.value < 0.62) |\n",
" (template_cut.spectral_axis.value > 0.70) & (template_cut.spectral_axis.value < 1.05))\n",
"\n",
"template_forcont = Spectrum1D(spectral_axis=template_cut.spectral_axis[mask_temp], flux=template_cut.flux[mask_temp])\n",
"\n",
"# Use fit_generic_continuum\n",
"fit_temp = fit_generic_continuum(template_forcont, model=Polynomial1D(5))\n",
"cont_temp = fit_temp(template_cut.spectral_axis)\n",
"template_sub = template_cut - cont_temp"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "4e8ad9cd",
"metadata": {},
"outputs": [],
"source": [
"# Print Spectrum1D object\n",
"print(template_sub)"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "5c36de57",
"metadata": {},
"outputs": [],
"source": [
"# Look at spectrum\n",
"plt.figure(figsize=[10, 6])\n",
"plt.plot(template_sub.spectral_axis, template_sub.flux)\n",
"plt.xlabel(\"wavelength ({:latex})\".format(template_sub.spectral_axis.unit))\n",
"plt.ylabel(\"flux ({:latex})\".format(template_sub.flux.unit))\n",
"plt.title(\"Continuum subtracted template\")\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"id": "ddcef91f",
"metadata": {},
"source": [
"### Subtract the continuum from the observed spectrum\n",
"\n",
"We can use a different approach and do it with [SpectralRegion](https://specutils.readthedocs.io/en/stable/spectral_regions.html) here. We also need to include an uncertianty to the observed spectrum, if it is not included."
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "dc598e8a",
"metadata": {},
"outputs": [],
"source": [
"# Read the spectrum directly as Spectrum1D\n",
"spec1d = Spectrum1D.read(fn_nirspec)\n",
"# Add the uncertainty\n",
"spec1d = Spectrum1D(spectral_axis=spec1d.spectral_axis, flux=spec1d.flux, \n",
" uncertainty=StdDevUncertainty(0.05*(spec1d.flux)))\n",
"\n",
"print(spec1d)\n",
"\n",
"# Define Spectral Region\n",
"region = SpectralRegion(1.3*u.um, 1.6*u.um)\n",
"# Extract region\n",
"spec1d_cont = extract_region(spec1d, region)\n",
"# Run fitting function\n",
"fit_obs = fit_generic_continuum(spec1d_cont, model=Linear1D(5))\n",
"# Apply to spectral axis\n",
"cont_obs = fit_obs(spec1d.spectral_axis)\n",
"\n",
"# Subtract continuum\n",
"spec1d_sub = spec1d - cont_obs"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "e8ae87d2",
"metadata": {},
"outputs": [],
"source": [
"# Look at spectrum\n",
"plt.figure(figsize=[10, 6])\n",
"plt.plot(spec1d_sub.spectral_axis, spec1d_sub.flux)\n",
"plt.xlabel(\"wavelength ({:latex})\".format(spec1d_sub.spectral_axis.unit))\n",
"plt.ylabel(\"flux ({:latex})\".format(spec1d_sub.flux.unit))\n",
"plt.title(\"Continuum subtracted observed spectrum\")\n",
"plt.show()"
]
},
{
"cell_type": "markdown",
"id": "4e4e6187",
"metadata": {},
"source": [
"### Exercise: fit the continuum using Specviz\n",
"\n",
"We can find and fit the continuum also in Specviz! You will have to:\n",
"- open a new instance of Specviz,\n",
"- load the template spectrum,\n",
"- [define the regions](https://jdaviz.readthedocs.io/en/latest/specviz/displaying.html#defining-spectral-regions) where there are no emission lines,\n",
"- and use the plugin [Model Fitting](https://jdaviz.readthedocs.io/en/latest/specviz/plugins.html#model-fitting) to measure the continuum."
]
},
{
"cell_type": "markdown",
"id": "8195cee8",
"metadata": {},
"source": [
"### Run the cross correlation function"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "586e9e52",
"metadata": {},
"outputs": [],
"source": [
"# Call the function\n",
"corr, lag = correlation.template_correlate(spec1d_sub, template_sub, lag_units=u.one)\n",
"\n",
"# Plot the correlation\n",
"plt.plot(lag, corr)\n",
"plt.xlabel(\"Redshift\")\n",
"plt.ylabel(\"Correlation\")\n",
"plt.show()"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "56659ed6",
"metadata": {},
"outputs": [],
"source": [
"# Redshift based on maximum\n",
"index_peak = np.argmax(corr)\n",
"z = lag[index_peak]\n",
"\n",
"print(\"Redshift from peak maximum: \", z)"
]
},
{
"cell_type": "markdown",
"id": "2e3e4101",
"metadata": {},
"source": [
"\n",
"## Template fitting using a library of spectra\n",
"\n",
"There are two other functions that are suitable for redshift fitting: [template_match](https://specutils.readthedocs.io/en/stable/api/specutils.analysis.template_match.html#specutils.analysis.template_match) and [template_redshift](https://specutils.readthedocs.io/en/stable/api/specutils.analysis.template_redshift.html#specutils.analysis.template_redshift). Here we use template_match which lives within [template_comparison](https://specutils.readthedocs.io/en/stable/analysis.html#template-comparison). We need to specify a library of spectra and the range of potential redshifts. The function computes a chi^2 for each model and returns the best metching model, either at fixed redshift or in a range of redshifts. The time this function takes depends heavily on the size of the library and the length of the redshift array. The library contains 100 model spectral from the same approach as our previous model template."
]
},
{
"cell_type": "markdown",
"id": "886dcb1b",
"metadata": {},
"source": [
"### Get the library of spectra"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "1a13c9f1",
"metadata": {},
"outputs": [],
"source": [
"# Download the library and unzip\n",
"boxlink = 'https://stsci.box.com/shared/static/1d68n0cg38it7ja3rgyvv8gxjezdqjgh.zip'\n",
"boxfile = './templates.zip'\n",
"urllib.request.urlretrieve(boxlink, boxfile)\n",
"\n",
"zf = zipfile.ZipFile(boxfile, 'r')\n",
"zf.extractall()\n",
"\n",
"templatedir = './templates/'"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "47933f29",
"metadata": {},
"outputs": [],
"source": [
"# Create the list of Spectrum1Ds for the templates\n",
"templatelist = []\n",
"for i in range(1, 100):\n",
" template_file = \"{0}{1:05d}.dat\".format(templatedir, i)\n",
" template = ascii.read(template_file)\n",
" temp1d = Spectrum1D(spectral_axis=template['col1']/1.E4*u.um, flux=template['col2']*spec_unit)\n",
" templatelist.append(temp1d)"
]
},
{
"cell_type": "markdown",
"id": "4eef5c68",
"metadata": {},
"source": [
"### Run the template_match function\n",
"\n",
"This takes a little while...it is time for your favorite drink!"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "6040c88c",
"metadata": {},
"outputs": [],
"source": [
"# You can either use a fixed redshift\n",
"zz = 1.501\n",
"\n",
"# Or you can define a range of redshifts (uncomment to use this option)\n",
"#min_redshift = 1.500\n",
"#max_redshift = 1.502\n",
"#delta_redshift = 0.0005\n",
"#zz = np.arange(min_redshift, max_redshift+delta_redshift, delta_redshift)\n",
"\n",
"# Define the resample method\n",
"resample_method = \"flux_conserving\"\n",
"\n",
"# Run the function\n",
"tm_results = template_comparison.template_match(observed_spectrum=spec1d, \n",
" spectral_templates=templatelist, \n",
" resample_method=resample_method, \n",
" redshift=zz) # Change the redshift variable for list or fixed"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "adff2994",
"metadata": {},
"outputs": [],
"source": [
"# Print results\n",
"print(tm_results)"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "a6a1e6c0",
"metadata": {},
"outputs": [],
"source": [
"# Load the best matching spectrum in Specviz\n",
"viz2 = Specviz()\n",
"viz2.app"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "49d13dca",
"metadata": {},
"outputs": [],
"source": [
"viz2.load_spectrum(spec1d, data_label='Observed spectrum')\n",
"viz2.load_spectrum(tm_results[0], data_label='Redshifted best template')"
]
},
{
"cell_type": "code",
"execution_count": null,
"id": "8a2ce17a",
"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": 5
}