MOS Spectroscopy of Extragalactic Field#

Use case: emission-line measurements and template matching on 1D spectra.
Data: CEERS NIRSpec observations
Tools: specutils, astropy, matplotlib, jdaviz.
Cross-intrument:
Documentation: This notebook is part of a STScI’s larger post-pipeline Data Analysis Tools Ecosystem.

Introduction#

In this notebook, we will inspect a set of spectra and perform a seris of spectroscopic analyses on an example spectrum, including continuum fitting and subtraction, line identification, centroiding and flux measurements, gaussian fitting, equivalent widths, and template fitting. We will do so using the interactive jdaviz package and the command line. We will use a JWST/NIRSpec spectroscopic dataset from the CEERS program.

Objective of the notebook#

The aim of the notebook is to showcase how to use the visualization tool Jdaviz or the combination Specutils+Matplotlib to measure the properties of the OII emission line in a NIRSpec spectrum.

Workflow#

  • visualize the spectroscopic dataset in Mosviz

  • select one galaxy (s02904) and visualize it in Specviz2d

  • perform the 1D extraction of the bright companion using the Spectral Extraction plugin in Specviz2d

  • attach the wavelength axis to the extracted 1D spectrum

  • select the OII emission line and measure

    • the redshift of the source

    • the properties of the emission line

  • fit a model with continuum + a Gaussian to the OII emission line

  • perform the same tasks using Specutils and Matplotlib instead of Jdaviz

  • find the best-fitting template of the observed spectrum

System requirements#

First, we create an environment with jdaviz which contains all the spectroscopic packages we need.

conda create -n jdaviz python
conda activate jdaviz
pip install jdaviz

Imports#

# general os
import zipfile
import urllib.request
from pathlib import Path

# general plotting
from matplotlib import pyplot as plt

# table/math handling
import numpy as np

# astropy
import astropy
import astropy.units as u
from astropy.io import fits, ascii
from astropy.nddata import StdDevUncertainty
from astropy.modeling import models
from astropy.visualization import quantity_support

# specutils
import specutils
from specutils import Spectrum1D, SpectralRegion
from specutils.fitting import fit_generic_continuum
from specutils.fitting import find_lines_threshold
from specutils.fitting import fit_lines
from specutils.manipulation import extract_region
from specutils.analysis import centroid
from specutils.analysis import line_flux
from specutils.analysis import equivalent_width
from specutils.analysis import template_comparison

# jdaviz
import jdaviz
from jdaviz import Mosviz, Specviz2d, Specviz  # noqa

# glue
from glue.core.roi import XRangeROI

np.seterr(all='ignore')  # hides irrelevant warnings about divide-by-zero, etc
quantity_support()  # auto-recognizes units on matplotlib plots

# Matplotlib parameters
params = {'legend.fontsize': '18', 
          'axes.labelsize': '18',
          'axes.titlesize': '18', 
          'xtick.labelsize': '18',
          'ytick.labelsize': '18', 
          'lines.linewidth': 2, 
          'axes.linewidth': 2, 
          'animation.html': 'html5'}
plt.rcParams.update(params)
plt.rcParams.update({'figure.max_open_warning': 0})

Check versions. Latest working environment is:#

Numpy: 2.1.0
Astropy: 6.1.2
Specutils: 1.16.0
Jdaviz: 3.10.3

print("Numpy: ", np.__version__)
print("Astropy: ", astropy.__version__)
print("Specutils: ", specutils.__version__)
print("Jdaviz: ", jdaviz.__version__)
Numpy:  2.1.3
Astropy:  6.1.6
Specutils:  1.19.0
Jdaviz:  4.0.0

Open data in Mosviz#

In Mosviz we can explore all the data products in the folder. This cell takes a minute or two to run. Since we are not including images, we can expand the 2D/1D spectra viewers to use the full width of the GUI. We can also keep the plugin tray open on metadata to check specifics of the files we are looking at.

Mosviz with its 2D spectrum, 1D spectrum, and table viewer.

Choose one galaxy#

We choose s02904 because we want to extract the very bright companion below the target.

file1d = pathtodata / 'jw01345-o064_s02904_nirspec_f100lp-g140m_x1d.fits'
file2d = pathtodata / 'jw01345-o064_s02904_nirspec_f100lp-g140m_s2d.fits'

Use Specviz2d to extract a better 1D spectrum#

specviz2d = Specviz2d()
specviz2d.load_data(file2d)
specviz2d.show()

Developer note
Is there a way to get out an uncertainty array from the extraction?
This is being worked on.

We open the Spectral Extraction plugin and select the appropriate trace (Polynomial, order 3, on pixel 2), background (Manual, on pixel 8, width 2, statistic average), and extraction (From Plugin, Horne). We then click Extract and inspect the extracted spectrum in the 1D viewer.

Specviz2d with spectral extraction plugin open showing the trace parameters. Specviz2d with spectral extraction plugin open showing the background parameters. Specviz2d with spectral extraction plugin open showing the extraction parameters.
# Get out extracted spectrum from Specviz2d
spectra = specviz2d.get_data('Spectrum 1D')
spectra
<Spectrum1D(flux=[0.0 ... 0.0] MJy (shape=(1325,), mean=0.00000 MJy); spectral_axis=<SpectralAxis 
   (observer to target:
      radial_velocity=0.0 km / s
      redshift=0.0)
  [0.96262388 0.96262388 0.96262388 ... 1.89519827 1.89583322 1.89583322] um> (length=1325))>
# Include some fake uncertainty for now
spec1d = Spectrum1D(spectral_axis=spectra.spectral_axis,
                    flux=spectra.flux,
                    uncertainty=StdDevUncertainty((np.zeros(len(spectra.flux)) + 1E-13) * spectra.unit))
spec1d
<Spectrum1D(flux=[0.0 ... 0.0] MJy (shape=(1325,), mean=0.00000 MJy); spectral_axis=<SpectralAxis 
   (observer to target:
      radial_velocity=0.0 km / s
      redshift=0.0)
  [0.96262388 0.96262388 0.96262388 ... 1.89519827 1.89583322 1.89583322] um> (length=1325); uncertainty=StdDevUncertainty)>
# And open in Specviz
specviz = Specviz()
specviz.load_data(spec1d, data_label='spec1d calibrated')
specviz.show()

There are still some artifacts in the data, but we can select a subset masking the artifacts and get out a spectrum without unwanted spikes. We can do so using the tool to select a subset with the “add” option (in the top bar) to select multiple regions as part of a single subset.

Specviz showing how to select a subset and choose the add option.
# Create a subset in the area of interest if it has not been created manually
try:
    region1 = specviz.get_data(data_label='spec1d calibrated', spectral_subset='Subset 1')
    print(region1)
    region1_exists = True
except Exception:
    print("There are no subsets selected.")
    region1_exists = False
    
# Spectral region for masking artifacts just around [OII]
if not region1_exists:
    sv = specviz.app.get_viewer('spectrum-viewer')
    sv.toolbar_active_subset.selected = []
    sv.apply_roi(XRangeROI(1., 1.27))  
There are no subsets selected.
# Get spectrum out with mask
spec1d_region = specviz.get_spectral_regions()
spec1d_masked = extract_region(spec1d, spec1d_region['Subset 1'], return_single_spectrum=True)
# Load in specviz
specviz.load_data(spec1d_masked, data_label='spec1d masked')
# Write the extracted spectrum to a fits file
file_extracted = Path('./extracted_spectrum.fits')
spec1d_masked.write(file_extracted, overwrite=True)
# Check that it has everything
hdu = fits.open(file_extracted)
hdu.info()
Filename: extracted_spectrum.fits
No.    Name      Ver    Type      Cards   Dimensions   Format
  0  PRIMARY       1 PrimaryHDU       4   ()      
  1  DATA          1 BinTableHDU     18   424R x 3C   [D, D, D]   
hdu[1].data
FITS_rec([(1.00020414, 2.20226191e-13, 1.e-13),
          (1.0008411 , 7.11024936e-14, 1.e-13),
          (1.00149188, 3.35875760e-13, 1.e-13),
          (1.00212885, 2.73494495e-13, 1.e-13),
          (1.00276581, 2.33289362e-13, 1.e-13),
          (1.00340278, 2.08752581e-13, 1.e-13),
          (1.00403974, 3.60731324e-13, 1.e-13),
          (1.00467671, 2.91254772e-13, 1.e-13),
          (1.00531368, 3.54084183e-13, 1.e-13),
          (1.00595064, 1.58528613e-13, 1.e-13),
          (1.00658761, 1.00099923e-13, 1.e-13),
          (1.00722458, 1.93899402e-13, 1.e-13),
          (1.00786154, 2.71433068e-13, 1.e-13),
          (1.00849851, 5.44897721e-13, 1.e-13),
          (1.00913548, 2.77380637e-13, 1.e-13),
          (1.00977245, 5.20373636e-13, 1.e-13),
          (1.01040941, 3.18918981e-13, 1.e-13),
          (1.01104638, 2.26125943e-13, 1.e-13),
          (1.01168335, 3.19381799e-13, 1.e-13),
          (1.01232032, 3.14886887e-13, 1.e-13),
          (1.01295729, 2.44143994e-13, 1.e-13),
          (1.01359426, 4.88193151e-13, 1.e-13),
          (1.01423123, 2.53360983e-13, 1.e-13),
          (1.0148682 , 1.27069879e-13, 1.e-13),
          (1.01550516, 4.64169909e-13, 1.e-13),
          (1.01614213, 2.08450736e-13, 1.e-13),
          (1.0167791 , 2.73294743e-13, 1.e-13),
          (1.01741607, 3.41261670e-13, 1.e-13),
          (1.01806689, 1.28366805e-13, 1.e-13),
          (1.01870387, 2.49072251e-13, 1.e-13),
          (1.01934084, 1.26547081e-13, 1.e-13),
          (1.01997781, 1.97832808e-13, 1.e-13),
          (1.02061478, 1.62252138e-13, 1.e-13),
          (1.02125175, 1.70550677e-13, 1.e-13),
          (1.02188873, 4.56591770e-13, 1.e-13),
          (1.0225257 , 3.02788739e-13, 1.e-13),
          (1.02316267, 2.49439997e-13, 1.e-13),
          (1.02379965, 3.97524604e-13, 1.e-13),
          (1.02443662, 2.53102488e-13, 1.e-13),
          (1.02507359, 3.03111230e-13, 1.e-13),
          (1.02571057, 6.63408469e-14, 1.e-13),
          (1.02634754, 4.50355456e-14, 1.e-13),
          (1.02698451, 1.39524997e-13, 1.e-13),
          (1.02762149, 4.23035735e-13, 1.e-13),
          (1.02825846, 1.30840441e-13, 1.e-13),
          (1.02889544, 2.38858164e-13, 1.e-13),
          (1.02953241, 1.90624463e-13, 1.e-13),
          (1.03016938, 1.56144866e-13, 1.e-13),
          (1.03080636, 3.47522079e-13, 1.e-13),
          (1.03144333, 3.05209873e-13, 1.e-13),
          (1.03208031, 2.66962446e-13, 1.e-13),
          (1.03271728, 4.45226992e-13, 1.e-13),
          (1.03335426, 3.47891271e-13, 1.e-13),
          (1.03399123, 2.72119563e-13, 1.e-13),
          (1.03462821, 2.85247721e-13, 1.e-13),
          (1.03526518, 1.95890246e-13, 1.e-13),
          (1.03590216, 3.64113112e-13, 1.e-13),
          (1.03653913, 5.22469764e-13, 1.e-13),
          (1.03717611, 3.54889614e-13, 1.e-13),
          (1.03781308, 3.96473035e-13, 1.e-13),
          (1.03845006, 2.99263516e-13, 1.e-13),
          (1.03908704, 7.15882448e-14, 1.e-13),
          (1.03972401, 1.87265535e-13, 1.e-13),
          (1.04036099, 1.80539887e-13, 1.e-13),
          (1.04099796, 5.28958335e-13, 1.e-13),
          (1.04163494, 1.51612494e-13, 1.e-13),
          (1.04227191, 2.32234510e-13, 1.e-13),
          (1.04290889, 4.02132384e-13, 1.e-13),
          (1.04354587, 3.60423596e-13, 1.e-13),
          (1.04418284, 4.07250692e-13, 1.e-13),
          (1.04481982, 2.05864649e-13, 1.e-13),
          (1.0454568 , 8.58324407e-14, 1.e-13),
          (1.04609377, 3.64159481e-13, 1.e-13),
          (1.04673075, 1.32323051e-13, 1.e-13),
          (1.04736772, 2.38120456e-13, 1.e-13),
          (1.0480047 , 1.67463119e-13, 1.e-13),
          (1.04864168, 2.68398002e-13, 1.e-13),
          (1.04927865, 2.70872084e-13, 1.e-13),
          (1.04991563, 4.49884622e-13, 1.e-13),
          (1.05055261, 2.90386076e-13, 1.e-13),
          (1.05118958, 1.84278016e-13, 1.e-13),
          (1.05182656, 1.97502237e-13, 1.e-13),
          (1.05246354, 1.79450524e-13, 1.e-13),
          (1.05310051, 3.16555189e-13, 1.e-13),
          (1.05373749, 2.83794762e-13, 1.e-13),
          (1.05437446, 3.58840671e-13, 1.e-13),
          (1.05501144, 1.01462592e-13, 1.e-13),
          (1.05564842, 2.98397874e-13, 1.e-13),
          (1.05628539, 3.26516193e-13, 1.e-13),
          (1.05692237, 1.01659088e-13, 1.e-13),
          (1.05755935, 3.04582978e-13, 1.e-13),
          (1.05819632, 3.33097866e-13, 1.e-13),
          (1.0588333 , 1.78090829e-13, 1.e-13),
          (1.05947028, 2.45500165e-13, 1.e-13),
          (1.06010725, 4.97942802e-13, 1.e-13),
          (1.06074423, 4.88389341e-13, 1.e-13),
          (1.06139514, 4.00679570e-13, 1.e-13),
          (1.06203211, 3.34417938e-13, 1.e-13),
          (1.06266909, 1.70919045e-13, 1.e-13),
          (1.06330607, 1.65230117e-13, 1.e-13),
          (1.06394305, 4.36515868e-13, 1.e-13),
          (1.06458002, 2.81471403e-13, 1.e-13),
          (1.065217  , 3.37830819e-13, 1.e-13),
          (1.06585398, 4.27911614e-13, 1.e-13),
          (1.06649096, 2.76036147e-13, 1.e-13),
          (1.06712793, 3.11917183e-13, 1.e-13),
          (1.06776491, 3.31327270e-13, 1.e-13),
          (1.06840189, 1.81340533e-13, 1.e-13),
          (1.06903886, 4.90888931e-13, 1.e-13),
          (1.06967584, 2.67252941e-13, 1.e-13),
          (1.07031282, 2.32996195e-13, 1.e-13),
          (1.07094979, 1.66874857e-13, 1.e-13),
          (1.07158677, 1.53642859e-13, 1.e-13),
          (1.07222375, 3.56233580e-13, 1.e-13),
          (1.07286072, 4.79251782e-13, 1.e-13),
          (1.0734977 , 5.22163418e-13, 1.e-13),
          (1.07413468, 4.41704329e-13, 1.e-13),
          (1.07477165, 3.26611081e-13, 1.e-13),
          (1.07540863, 2.21805842e-13, 1.e-13),
          (1.0760456 , 2.31779499e-13, 1.e-13),
          (1.07668258, 3.57814161e-13, 1.e-13),
          (1.07731955, 4.05631637e-13, 1.e-13),
          (1.07795653, 4.09391098e-15, 1.e-13),
          (1.07860747, 7.48259747e-14, 1.e-13),
          (1.07924445, 3.82162280e-13, 1.e-13),
          (1.07988142, 2.79952022e-13, 1.e-13),
          (1.0805184 , 3.00224870e-13, 1.e-13),
          (1.08115538, 2.80047300e-13, 1.e-13),
          (1.08179235, 3.50226498e-13, 1.e-13),
          (1.08242933, 2.81212596e-13, 1.e-13),
          (1.0830663 , 1.69329891e-13, 1.e-13),
          (1.08370328, 2.98331475e-13, 1.e-13),
          (1.08434026, 2.02522585e-13, 1.e-13),
          (1.08497723, 2.38028432e-13, 1.e-13),
          (1.08561421, 2.94623541e-13, 1.e-13),
          (1.08625118, 1.98378133e-13, 1.e-13),
          (1.08688816, 2.79912249e-13, 1.e-13),
          (1.08752513, 4.76740118e-13, 1.e-13),
          (1.08816211, 3.27518460e-13, 1.e-13),
          (1.08879908, 3.77980941e-13, 1.e-13),
          (1.08943605, 1.56084517e-13, 1.e-13),
          (1.09007303, 4.87805261e-13, 1.e-13),
          (1.09071   , 4.12587606e-13, 1.e-13),
          (1.09134698, 4.41010049e-13, 1.e-13),
          (1.09198395, 3.51925040e-13, 1.e-13),
          (1.09262092, 4.43774334e-13, 1.e-13),
          (1.0932579 , 5.36644863e-13, 1.e-13),
          (1.09389487, 3.53419892e-13, 1.e-13),
          (1.09453184, 3.02046708e-13, 1.e-13),
          (1.09516882, 4.60493155e-13, 1.e-13),
          (1.09580579, 4.85328159e-13, 1.e-13),
          (1.09644276, 4.03150577e-13, 1.e-13),
          (1.09707973, 4.66453043e-13, 1.e-13),
          (1.0977167 , 4.20896813e-13, 1.e-13),
          (1.09835368, 4.34047978e-13, 1.e-13),
          (1.09899065, 2.69152230e-13, 1.e-13),
          (1.09962762, 5.57799040e-13, 1.e-13),
          (1.10026459, 3.75521510e-13, 1.e-13),
          (1.10090156, 1.60157684e-13, 1.e-13),
          (1.10153853, 3.90056351e-13, 1.e-13),
          (1.1021755 , 4.49673205e-13, 1.e-13),
          (1.10281247, 4.17112815e-13, 1.e-13),
          (1.10344944, 5.57487304e-13, 1.e-13),
          (1.10408641, 4.72369144e-13, 1.e-13),
          (1.10472338, 3.83054709e-13, 1.e-13),
          (1.10536035, 4.03284496e-13, 1.e-13),
          (1.10599732, 3.49443909e-13, 1.e-13),
          (1.10663429, 4.14028603e-13, 1.e-13),
          (1.10727126, 2.98529630e-13, 1.e-13),
          (1.10790822, 2.23500643e-13, 1.e-13),
          (1.10854519, 2.72224139e-13, 1.e-13),
          (1.10918216, 2.75031227e-13, 1.e-13),
          (1.10981913, 4.49926698e-13, 1.e-13),
          (1.1104561 , 5.17033445e-13, 1.e-13),
          (1.11109306, 3.77205709e-13, 1.e-13),
          (1.11173003, 5.14777195e-13, 1.e-13),
          (1.112367  , 3.45728148e-13, 1.e-13),
          (1.11300396, 4.78515955e-13, 1.e-13),
          (1.11364093, 4.41022423e-13, 1.e-13),
          (1.11427789, 3.15633798e-13, 1.e-13),
          (1.11491486, 4.60366515e-13, 1.e-13),
          (1.11555182, 4.59537498e-13, 1.e-13),
          (1.11618879, 3.57125795e-13, 1.e-13),
          (1.11682575, 4.71355585e-13, 1.e-13),
          (1.11746272, 5.39477170e-13, 1.e-13),
          (1.11809968, 4.71055690e-13, 1.e-13),
          (1.11873664, 2.58177902e-13, 1.e-13),
          (1.11937361, 1.64687286e-13, 1.e-13),
          (1.12001057, 4.36668862e-13, 1.e-13),
          (1.12064753, 3.92515244e-13, 1.e-13),
          (1.1212845 , 2.27476382e-13, 1.e-13),
          (1.12192146, 3.18172037e-13, 1.e-13),
          (1.12255842, 4.06247783e-13, 1.e-13),
          (1.12319538, 4.13594863e-13, 1.e-13),
          (1.12383234, 3.45983392e-13, 1.e-13),
          (1.1244693 , 4.04505634e-13, 1.e-13),
          (1.12510626, 3.43614649e-13, 1.e-13),
          (1.12574322, 6.68088553e-13, 1.e-13),
          (1.12639424, 7.26797843e-13, 1.e-13),
          (1.1270312 , 1.02647395e-12, 1.e-13),
          (1.12766816, 1.35566790e-12, 1.e-13),
          (1.12830512, 1.11219154e-12, 1.e-13),
          (1.12894208, 7.09701231e-13, 1.e-13),
          (1.12957904, 5.88944366e-13, 1.e-13),
          (1.130216  , 3.44075404e-13, 1.e-13),
          (1.13085296, 4.79117942e-13, 1.e-13),
          (1.13148991, 3.68676437e-13, 1.e-13),
          (1.13212687, 4.17866740e-13, 1.e-13),
          (1.13276383, 3.21985728e-13, 1.e-13),
          (1.13340079, 3.89816289e-13, 1.e-13),
          (1.13403774, 3.65957736e-13, 1.e-13),
          (1.1346747 , 3.33283158e-13, 1.e-13),
          (1.13531166, 3.55083953e-13, 1.e-13),
          (1.13594861, 5.91551296e-13, 1.e-13),
          (1.13658557, 4.47070259e-13, 1.e-13),
          (1.13722252, 5.73471459e-13, 1.e-13),
          (1.13785948, 5.68602253e-13, 1.e-13),
          (1.13849643, 5.00111671e-13, 1.e-13),
          (1.13913339, 4.35416450e-13, 1.e-13),
          (1.13977034, 2.84418065e-13, 1.e-13),
          (1.14040729, 4.52535860e-13, 1.e-13),
          (1.14104425, 5.45833725e-13, 1.e-13),
          (1.1416812 , 4.96246994e-13, 1.e-13),
          (1.14231815, 3.36758488e-13, 1.e-13),
          (1.1429551 , 5.51705227e-13, 1.e-13),
          (1.14359205, 7.45589405e-13, 1.e-13),
          (1.144229  , 6.98371932e-13, 1.e-13),
          (1.14486595, 6.56632122e-13, 1.e-13),
          (1.14551699, 7.82987859e-13, 1.e-13),
          (1.14615394, 6.46901991e-13, 1.e-13),
          (1.14679089, 5.42259047e-13, 1.e-13),
          (1.14742784, 4.51682342e-13, 1.e-13),
          (1.14806479, 4.81228676e-13, 1.e-13),
          (1.14870174, 5.44249697e-13, 1.e-13),
          (1.14933869, 3.04091916e-13, 1.e-13),
          (1.14997564, 4.05739575e-13, 1.e-13),
          (1.15061259, 4.56585692e-13, 1.e-13),
          (1.15124953, 4.88541726e-13, 1.e-13),
          (1.15188648, 5.04249058e-13, 1.e-13),
          (1.15252343, 5.77207197e-13, 1.e-13),
          (1.15316037, 5.67643641e-13, 1.e-13),
          (1.15379732, 7.93678610e-13, 1.e-13),
          (1.15443426, 7.54212908e-13, 1.e-13),
          (1.15507121, 8.45107557e-13, 1.e-13),
          (1.15570815, 8.26016970e-13, 1.e-13),
          (1.1563451 , 8.62169408e-13, 1.e-13),
          (1.15698204, 6.18060482e-13, 1.e-13),
          (1.15761898, 6.24817299e-13, 1.e-13),
          (1.15825592, 5.32213718e-13, 1.e-13),
          (1.15889287, 6.99702682e-13, 1.e-13),
          (1.15952981, 6.52652994e-13, 1.e-13),
          (1.16016675, 5.38442732e-13, 1.e-13),
          (1.16080369, 4.31362150e-13, 1.e-13),
          (1.16144063, 3.40311832e-13, 1.e-13),
          (1.16207757, 4.29430448e-13, 1.e-13),
          (1.16271451, 6.20946174e-13, 1.e-13),
          (1.16335144, 6.21408236e-13, 1.e-13),
          (1.16398838, 6.93829160e-13, 1.e-13),
          (1.16462532, 7.22670083e-13, 1.e-13),
          (1.16526226, 6.57540181e-13, 1.e-13),
          (1.16589919, 6.15419791e-13, 1.e-13),
          (1.16653613, 6.91254856e-13, 1.e-13),
          (1.16717306, 6.55687150e-13, 1.e-13),
          (1.16781   , 9.88001222e-13, 1.e-13),
          (1.16844693, 8.89964926e-13, 1.e-13),
          (1.16908387, 9.15438730e-13, 1.e-13),
          (1.1697208 , 9.79759846e-13, 1.e-13),
          (1.17035773, 9.94234731e-13, 1.e-13),
          (1.17099467, 1.01260039e-12, 1.e-13),
          (1.1716316 , 9.01097130e-13, 1.e-13),
          (1.17226853, 1.06732534e-12, 1.e-13),
          (1.17290546, 8.26004028e-13, 1.e-13),
          (1.17354239, 7.37443931e-13, 1.e-13),
          (1.17417932, 6.34469860e-13, 1.e-13),
          (1.17481625, 6.11196615e-13, 1.e-13),
          (1.17545318, 6.07083641e-13, 1.e-13),
          (1.1760901 , 4.84819533e-13, 1.e-13),
          (1.17672703, 2.52542002e-13, 1.e-13),
          (1.17736396, 3.62636719e-13, 1.e-13),
          (1.17800088, 3.58695834e-13, 1.e-13),
          (1.17863781, 5.30236188e-13, 1.e-13),
          (1.17927474, 5.89884451e-13, 1.e-13),
          (1.17991166, 5.28850538e-13, 1.e-13),
          (1.18054858, 4.87772924e-13, 1.e-13),
          (1.18118551, 6.05768150e-13, 1.e-13),
          (1.18182243, 8.04558113e-13, 1.e-13),
          (1.18245935, 1.15401964e-12, 1.e-13),
          (1.18309628, 8.73469938e-13, 1.e-13),
          (1.1837332 , 8.90579680e-13, 1.e-13),
          (1.18437012, 8.77589475e-13, 1.e-13),
          (1.18500704, 8.13590424e-13, 1.e-13),
          (1.18564396, 8.72015191e-13, 1.e-13),
          (1.18628088, 8.68059519e-13, 1.e-13),
          (1.18691779, 1.07781649e-12, 1.e-13),
          (1.18755471, 8.81699998e-13, 1.e-13),
          (1.18819163, 5.80767481e-13, 1.e-13),
          (1.18882855, 5.71175568e-13, 1.e-13),
          (1.18946546, 5.93353996e-13, 1.e-13),
          (1.19010238, 6.01578988e-13, 1.e-13),
          (1.19073929, 4.78701849e-13, 1.e-13),
          (1.19137621, 7.66148406e-13, 1.e-13),
          (1.19201312, 7.66950918e-13, 1.e-13),
          (1.19265003, 6.60555840e-13, 1.e-13),
          (1.19328694, 7.04460744e-13, 1.e-13),
          (1.19392386, 6.15648398e-13, 1.e-13),
          (1.19456077, 8.09748504e-13, 1.e-13),
          (1.19519768, 6.85802060e-13, 1.e-13),
          (1.19583459, 7.09985663e-13, 1.e-13),
          (1.1964715 , 7.36797096e-13, 1.e-13),
          (1.19710841, 6.76080645e-13, 1.e-13),
          (1.19774531, 6.74077883e-13, 1.e-13),
          (1.19838222, 6.16646818e-13, 1.e-13),
          (1.19903332, 5.78564532e-13, 1.e-13),
          (1.19967022, 3.85148788e-13, 1.e-13),
          (1.20030713, 3.58329469e-13, 1.e-13),
          (1.20094404, 6.88216371e-13, 1.e-13),
          (1.20158094, 2.01044678e-13, 1.e-13),
          (1.20221785, 2.91833933e-13, 1.e-13),
          (1.20285475, 5.08608957e-13, 1.e-13),
          (1.20349165, 6.75650298e-13, 1.e-13),
          (1.20412856, 7.19989477e-13, 1.e-13),
          (1.20476546, 8.48415740e-13, 1.e-13),
          (1.20540236, 6.47505770e-13, 1.e-13),
          (1.20603926, 6.06798682e-13, 1.e-13),
          (1.20667616, 5.65372771e-13, 1.e-13),
          (1.20731306, 6.66720700e-13, 1.e-13),
          (1.20794996, 6.64022075e-13, 1.e-13),
          (1.20858686, 7.96373765e-13, 1.e-13),
          (1.20922375, 7.83601084e-13, 1.e-13),
          (1.20986065, 6.28620782e-13, 1.e-13),
          (1.21049755, 6.65820264e-13, 1.e-13),
          (1.21113444, 8.10991477e-13, 1.e-13),
          (1.21177134, 1.12044211e-12, 1.e-13),
          (1.21240823, 1.23266757e-12, 1.e-13),
          (1.21304512, 1.30867330e-12, 1.e-13),
          (1.21368202, 1.11527540e-12, 1.e-13),
          (1.21431891, 7.02505830e-13, 1.e-13),
          (1.2149558 , 7.58094836e-13, 1.e-13),
          (1.21559269, 8.70156171e-13, 1.e-13),
          (1.21622958, 6.99432686e-13, 1.e-13),
          (1.21686647, 1.02199437e-12, 1.e-13),
          (1.21750336, 8.14122151e-13, 1.e-13),
          (1.21814024, 9.30280485e-13, 1.e-13),
          (1.21877713, 9.36998475e-13, 1.e-13),
          (1.21941402, 8.79586976e-13, 1.e-13),
          (1.22006513, 1.09553431e-12, 1.e-13),
          (1.22070202, 1.08331052e-12, 1.e-13),
          (1.2213389 , 1.15951308e-12, 1.e-13),
          (1.22197579, 1.01063784e-12, 1.e-13),
          (1.22261267, 8.00462840e-13, 1.e-13),
          (1.22324955, 9.45203676e-13, 1.e-13),
          (1.22388643, 8.56253876e-13, 1.e-13),
          (1.22452332, 1.03363840e-12, 1.e-13),
          (1.2251602 , 8.09953287e-13, 1.e-13),
          (1.22579708, 8.30163275e-13, 1.e-13),
          (1.22643396, 8.10837625e-13, 1.e-13),
          (1.22707083, 8.64276105e-13, 1.e-13),
          (1.22770771, 8.32194438e-13, 1.e-13),
          (1.22834459, 9.38323343e-13, 1.e-13),
          (1.22898146, 8.48516739e-13, 1.e-13),
          (1.22961834, 9.24011939e-13, 1.e-13),
          (1.23025521, 7.40206422e-13, 1.e-13),
          (1.23089209, 8.71513488e-13, 1.e-13),
          (1.23152896, 1.02335713e-12, 1.e-13),
          (1.23216583, 1.08035828e-12, 1.e-13),
          (1.23280271, 1.02178728e-12, 1.e-13),
          (1.23343958, 8.96742101e-13, 1.e-13),
          (1.23407645, 7.25233319e-13, 1.e-13),
          (1.23471332, 9.40153242e-13, 1.e-13),
          (1.23535018, 8.89781131e-13, 1.e-13),
          (1.23598705, 7.44040690e-13, 1.e-13),
          (1.23662392, 7.71931994e-13, 1.e-13),
          (1.23726078, 5.98854781e-13, 1.e-13),
          (1.23789765, 9.15759433e-13, 1.e-13),
          (1.23853451, 7.35240220e-13, 1.e-13),
          (1.23917138, 5.97180351e-13, 1.e-13),
          (1.23980824, 5.53652167e-13, 1.e-13),
          (1.2404451 , 5.19108648e-13, 1.e-13),
          (1.24108196, 6.35684304e-13, 1.e-13),
          (1.24171882, 8.64535789e-13, 1.e-13),
          (1.24235568, 5.54112249e-13, 1.e-13),
          (1.24299254, 6.85167108e-13, 1.e-13),
          (1.2436294 , 9.77341773e-13, 1.e-13),
          (1.24426626, 7.49452161e-13, 1.e-13),
          (1.24490311, 7.16886097e-13, 1.e-13),
          (1.24553997, 1.00827328e-12, 1.e-13),
          (1.24617683, 5.48718152e-13, 1.e-13),
          (1.24681368, 7.68007073e-13, 1.e-13),
          (1.24745053, 6.22569841e-13, 1.e-13),
          (1.24808738, 6.85560648e-13, 1.e-13),
          (1.24872424, 1.13240405e-12, 1.e-13),
          (1.24936109, 9.34135506e-13, 1.e-13),
          (1.24999794, 6.59118950e-13, 1.e-13),
          (1.25063478, 6.61413653e-13, 1.e-13),
          (1.25127163, 8.38442593e-13, 1.e-13),
          (1.25190848, 8.00919548e-13, 1.e-13),
          (1.25254533, 8.79331486e-13, 1.e-13),
          (1.25318217, 8.47567667e-13, 1.e-13),
          (1.25381902, 9.14539723e-13, 1.e-13),
          (1.25445586, 9.67988537e-13, 1.e-13),
          (1.2550927 , 8.84173404e-13, 1.e-13),
          (1.25572955, 8.71136028e-13, 1.e-13),
          (1.25636639, 9.17617424e-13, 1.e-13),
          (1.25700323, 7.46035154e-13, 1.e-13),
          (1.25764007, 1.10956686e-12, 1.e-13),
          (1.2582769 , 1.06695331e-12, 1.e-13),
          (1.25891374, 8.98392197e-13, 1.e-13),
          (1.25955058, 7.14786070e-13, 1.e-13),
          (1.26018742, 9.03643173e-13, 1.e-13),
          (1.26082425, 9.52279988e-13, 1.e-13),
          (1.26146108, 8.62454019e-13, 1.e-13),
          (1.26209792, 9.21385679e-13, 1.e-13),
          (1.26273475, 1.11029095e-12, 1.e-13),
          (1.26337158, 8.66299543e-13, 1.e-13),
          (1.26400841, 9.89356609e-13, 1.e-13),
          (1.26464524, 9.28098908e-13, 1.e-13),
          (1.26528207, 9.11408630e-13, 1.e-13),
          (1.2659189 , 8.56994858e-13, 1.e-13),
          (1.26655573, 7.50266908e-13, 1.e-13),
          (1.26719255, 7.39715298e-13, 1.e-13),
          (1.26782938, 7.64548293e-13, 1.e-13),
          (1.2684662 , 1.06161603e-12, 1.e-13),
          (1.26910303, 9.15607191e-13, 1.e-13),
          (1.26973985, 9.27040003e-13, 1.e-13)],
         dtype=(numpy.record, [('wavelength', '>f8'), ('flux', '>f8'), ('uncertainty', '>f8')]))

Workflow via API calls#

I can do some analysis on the spectrum using the plugins in the GUI. For reproducibility, I can do the same thing from the API, changing the parameters in the plugins programmatically.

Select regions#

# Select a region in the spectrum
sv = specviz.app.get_viewer('spectrum-viewer')
# Region with just line for line analysis
sv.toolbar_active_subset.selected = []
sv.apply_roi(XRangeROI(1.124, 1.131))  
# Region with some continuum for gaussian fit
sv.toolbar_active_subset.selected = []
sv.apply_roi(XRangeROI(1.05, 1.25))  

Line analysis#

Specviz showing the region on the [OII] line and the line analysis plugin.
# Open line analysis plugin
plugin_la = specviz.plugins['Line Analysis']
plugin_la.open_in_tray()
# List what's in the data menu
specviz.data_labels
['spec1d calibrated', 'spec1d masked']
# Input the appropriate spectrum and region
plugin_la.dataset = 'spec1d masked'
plugin_la.spectral_subset = 'Subset 2'
# Input the values for the continuum
plugin_la.continuum = 'Surrounding'
plugin_la.continuum_width = 7
# Return line analysis results
plugin_la.get_results()
[{'function': 'Line Flux', 'result': ''},
 {'function': 'Equivalent Width', 'result': ''},
 {'function': 'Gaussian Sigma Width', 'result': ''},
 {'function': 'Gaussian FWHM', 'result': ''},
 {'function': 'Centroid', 'result': ''}]

Line lists#

Specviz showing the line list plugin and the redshift identified with the [OII] doublet.
# Open line list plugin
plugin_ll = specviz.plugins['Line Lists']
plugin_ll.open_in_tray()

Developer note
The line list plugin cannot yet be accessed by the notebook. I can do it in the GUI though.

Open line list plugin. Select the SDSS IV line list. Load the Oxygen II lines and Hb. Go back to line analysis plugin and associate the Oxygen II line with the line we just analyzed.

Model fitting#

Specviz showing a model fit to the continuum and the [OII] line.
# Open model fitting plugin
plugin_mf = specviz.plugins['Model Fitting']
plugin_mf.open_in_tray()
# Input the appropriate datasets
plugin_mf.dataset = 'spec1d masked'
plugin_mf.spectral_subset = 'Subset 3'
# Input the model components
plugin_mf.create_model_component(model_component='Polynomial1D',
                                 poly_order=2,
                                 model_component_label='P2')
plugin_mf.create_model_component(model_component='Gaussian1D',
                                 model_component_label='G')
plugin_mf.get_model_component('G')
{'model_type': 'Gaussian1D',
 'parameters': {'amplitude': {'value': np.float64(9.28889869794063e-13),
   'unit': 'MJy',
   'std': nan,
   'fixed': False},
  'mean': {'value': np.float64(1.1695679004508002),
   'unit': 'um',
   'std': nan,
   'fixed': False},
  'stddev': {'value': np.float64(0.04850103358252829),
   'unit': 'um',
   'std': nan,
   'fixed': False}}}
plugin_mf.set_model_component('G', 'stddev', 0.001)
plugin_mf.set_model_component('G', 'mean', 1.128)
{'name': 'mean', 'value': 1.128, 'unit': 'um', 'fixed': False}
# Model equation gets populated automatically
plugin_mf.equation = 'P2+G'
# After we run this, we go to the GUI and check that the fit makes sense
plugin_mf.calculate_fit()
(<CompoundModel(c0_0=-0. MJy, c1_0=0. MJy / um, c2_0=-0. MJy / um2, amplitude_1=0. MJy, mean_1=1.12768736 um, stddev_1=0.00080152 um)>,
 <Spectrum1D(flux=[1.804397812792741e-14 ... 9.378888006579073e-13] MJy (shape=(424,), mean=0.00000 MJy); spectral_axis=<SpectralAxis [1.00020414 1.0008411  1.00149188 ... 1.2684662  1.26910303 1.26973985] um> (length=424))>)
specviz.get_model_parameters()
{'spec1d masked model': {'c0_0': <Quantity -6.08177385e-12 MJy>,
  'c1_0': <Quantity 8.21429761e-12 MJy / um>,
  'c2_0': <Quantity -2.11529292e-12 MJy / um2>,
  'amplitude_1': <Quantity 8.29374162e-13 MJy>,
  'mean_1': <Quantity 1.12768736 um>,
  'stddev_1': <Quantity 0.00080152 um>}}

Same workflow with specutils (old workflow)#

The same workflow can be achieved using directly the package specutils (which is used under the hood in jdaviz) and using matplotlib for a static visualization.

Fit and substract the continuum#

cont_spec1d = fit_generic_continuum(spec1d_masked)
cont_fit = cont_spec1d(spec1d_masked.spectral_axis)
WARNING: Model is linear in parameters; consider using linear fitting methods. [astropy.modeling.fitting]
WARNING:astroquery:Model is linear in parameters; consider using linear fitting methods.
plt.figure(figsize=[10, 6])
plt.plot(spec1d_masked.spectral_axis, spec1d_masked.flux, label="data")
plt.plot(spec1d_masked.spectral_axis, cont_fit, label="modeled continuum")
plt.xlabel("wavelength ({:latex})".format(spec1d_masked.spectral_axis.unit))
plt.ylabel("flux ({:latex})".format(spec1d_masked.flux.unit))
plt.legend()
plt.title("Observed spectrum and fitted continuum")
plt.show()

plt.figure(figsize=[10, 6])
plt.plot(spec1d_masked.spectral_axis, spec1d_masked.uncertainty.array, label="data")
plt.xlabel("wavelength ({:latex})".format(spec1d_masked.spectral_axis.unit))
plt.ylabel("uncertainty ({:latex})".format(spec1d_masked.uncertainty.unit))
plt.legend()
plt.title("Uncertianty of observed spectrum")
plt.show()
../../../_images/3f1bced85b4436ccf7fd58934d8f5948620e6dd6a3ca7895978da097509fa533.png ../../../_images/9cbf125e50c1864dcbca2b18ef92a1430e7c97db7c10c84e0307a8a05fcaaab9.png

Creating the continuum-subtracted spectrum#

Specutils will figure out what to do with the uncertainty!

spec1d_sub = spec1d_masked - cont_fit
spec1d_sub
<Spectrum1D(flux=[2.2022619120573399e-13 ... 9.270400033289505e-13] MJy (shape=(424,), mean=0.00000 MJy); spectral_axis=<SpectralAxis 
   (observer to target:
      radial_velocity=0.0 km / s
      redshift=0.0)
  [1.00020414 1.0008411  1.00149188 ... 1.2684662  1.26910303 1.26973985] um> (length=424); uncertainty=StdDevUncertainty)>
plt.figure(figsize=[10, 6])
plt.plot(spec1d_sub.spectral_axis, spec1d_sub.flux, label="data")
plt.xlabel("wavelength ({:latex})".format(spec1d_sub.spectral_axis.unit))
plt.ylabel("flux ({:latex})".format(spec1d_sub.flux.unit))
plt.legend()
plt.title("Continuum-subracted spectrum")
plt.show()

plt.figure(figsize=[10, 6])
plt.plot(spec1d_sub.spectral_axis, spec1d_sub.uncertainty.array, label="data")
plt.xlabel("wavelength ({:latex})".format(spec1d_sub.spectral_axis.unit))
plt.ylabel("uncertainty ({:latex})".format(spec1d_sub.uncertainty.unit))
plt.legend()
plt.title("Uncertainty of continuum-subracted spectrum")
plt.show()
../../../_images/004ec5a4a26889cf8fa5315857704c6993b5d16683d86ab00b86c70990c33741.png ../../../_images/83fd8ac1fe37923beea6a4aaff916563344f1caf22e07986f9a32208f07748a0.png

Find emission and absorption lines#

lines = find_lines_threshold(spec1d_sub, noise_factor=3)
lines
WARNING: Spectrum is not below the threshold signal-to-noise 0.01. This may indicate you have not continuum subtracted this spectrum (or that you have but it has high SNR features).

If you want to suppress this warning either type 'specutils.conf.do_continuum_function_check = False' or see http://docs.astropy.org/en/stable/config/#adding-new-configuration-items for other ways to configure the warning. [specutils.analysis.flux]
WARNING:astroquery:Spectrum is not below the threshold signal-to-noise 0.01. This may indicate you have not continuum subtracted this spectrum (or that you have but it has high SNR features).

If you want to suppress this warning either type 'specutils.conf.do_continuum_function_check = False' or see http://docs.astropy.org/en/stable/config/#adding-new-configuration-items for other ways to configure the warning.
QTable length=44
line_centerline_typeline_center_index
um
float64str8int64
1.001491881561837emission2
1.0040397430718908emission6
1.0053136752958713emission8
1.0084985099349837emission13
1.00977244535525emission15
1.0116833500823965emission18
1.0135942566493585emission21
1.015505164970485emission24
1.0174160749602992emission27
.........
1.0932578967661883emission146
1.099627618702622emission156
1.103449441777272emission162
1.1174627171235694emission184
1.1200105709685617emission188
1.1276681596792193emission200
1.1722685279834808emission270
1.1824593532966197emission286
1.2130451230160508emission334

Plot the emission lines on the spectrum.

plt.figure(figsize=[10, 6])
plt.plot(spec1d_sub.spectral_axis, spec1d_sub.flux, label="data")
plt.axvline(lines['line_center'][0].value, color="red", alpha=0.5, label='emission lines')
for line in lines:
    if line['line_type'] == 'emission':
        plt.axvline(line['line_center'].value, color='red', alpha=0.5)
plt.xlabel("wavelength ({:latex})".format(spec1d_sub.spectral_axis.unit))
plt.ylabel("flux ({:latex})".format(spec1d_sub.flux.unit))
plt.legend()
plt.title("Continuum-subtracted spectrum and marked lines using find_lines_threshold")
plt.show()
../../../_images/455dc558b8e9a9a6a9c394170b52295c9c0f938e5413c29df8e059fe038cf5df.png

Work by hand on a single line.

# Define limits for plotting
x_min = 1.1
x_max = 1.16

# Define limits for line region
line_min = 1.124*u.um
line_max = 1.131*u.um
plt.figure(figsize=[10, 6])
plt.plot(spec1d_sub.spectral_axis, spec1d_sub.flux, label="data")
plt.scatter(spec1d_sub.spectral_axis, spec1d_sub.flux, label=None)
plt.axvline(lines['line_center'][0].value, color="red", alpha=0.5, label='[OII]')
for line in lines:
    if line['line_type'] == 'emission':
        plt.axvline(line['line_center'].value, alpha=0.5, color='red')
plt.xlim(x_min, x_max)
plt.xlabel("wavelength ({:latex})".format(spec1d_sub.spectral_axis.unit))
plt.ylabel("flux ({:latex})".format(spec1d_sub.flux.unit))
plt.legend()
plt.title("Continuum-subtracted spectrum zoomed on [OII]")
plt.show()
../../../_images/0e319f3715c65fa353fbf5a555d15eb41ff208b8b6e09772d0efacd2150fbc7f.png

Measure line centroids and fluxes#

# Example with just one line
centroid(spec1d_sub, SpectralRegion(line_min, line_max))
\[1.1276811 \; \mathrm{\mu m}\]
sline = centroid(spec1d_sub, SpectralRegion(line_min, line_max))

plt.figure(figsize=[10, 6])
plt.plot(spec1d_sub.spectral_axis, spec1d_sub.flux, label="data")
plt.scatter(spec1d_sub.spectral_axis, spec1d_sub.flux, label=None)
plt.axvline(sline.value, color='red', label="[OII]")
plt.axhline(0, color='black', label='flux = 0')
plt.xlim(x_min, x_max)
plt.xlabel("wavelength ({:latex})".format(spec1d_sub.spectral_axis.unit))
plt.ylabel("flux ({:latex})".format(spec1d_sub.flux.unit))
plt.legend()
plt.title("Continuum-subtracted spectrum zoomed on [OII]")
plt.show()
../../../_images/841655d309ad9ed973edc2ebf0fd3cbf0ea0a3779c51278df2ac9d233e59aba9.png
line_flux(spec1d_sub, SpectralRegion(line_min, line_max))  
\[4.9520884 \times 10^{-15} \; \mathrm{MJy\,\mu m}\]

Fit the line with a Gaussian#

g_init = models.Gaussian1D(mean=1.1278909*u.um, stddev=0.001*u.um)
g_fit = fit_lines(spec1d_sub, g_init)
spec1d_fit = g_fit(spec1d_sub.spectral_axis)
g_fit
<Gaussian1D(amplitude=0. MJy, mean=1.1278909 um, stddev=0.001 um)>
plt.figure(figsize=[10, 6])
plt.plot(spec1d_sub.spectral_axis, spec1d_sub.flux, label='data')
plt.plot(spec1d_sub.spectral_axis, spec1d_fit, color='darkorange', label='Gaussian fit')
plt.xlim(x_min, x_max)
plt.xlabel("wavelength ({:latex})".format(spec1d_sub.spectral_axis.unit))
plt.ylabel("flux ({:latex})".format(spec1d_sub.flux.unit))
plt.legend()
plt.title('Gaussian fit to the [OII] line')
plt.show()
../../../_images/e53f21fbd6404bc6d9be3da1dd5d74865bf3846a8f11b2079faca4e3bfa7e30e.png

Measure the equivalent width of the lines#

This needs the spectrum continuum normalized.

spec1d_norm = spec1d_masked / cont_fit
plt.figure(figsize=[10, 6])
plt.plot(spec1d_norm.spectral_axis, spec1d_norm.flux, label='data')
plt.axhline(1, color='black', label='flux = 1')
plt.xlabel("wavelength ({:latex})".format(spec1d_norm.spectral_axis.unit))
plt.ylabel("flux (normalized)")
plt.xlim(x_min, x_max)
plt.legend()
plt.title("Continuum-normalized spectrum, zoomed on [OII]")
plt.show()

plt.figure(figsize=[10, 6])
plt.plot(spec1d_norm.spectral_axis, spec1d_norm.uncertainty.array, label='data')
plt.xlabel("wavelength ({:latex})".format(spec1d_norm.spectral_axis.unit))
plt.ylabel("uncertainty (normalized)")
plt.xlim(x_min, x_max)
plt.legend()
plt.title("Uncertainty of continuum-normalized spectrum, zoomed on [OII]")
plt.show()
../../../_images/0b5867e533be6c269d5a2cc96013144e825b2f016796059594d42542a092363f.png ../../../_images/e6e8c10681a1616acb90d9de8d6f2912df24ebb946412e755e10b2b53baccf71.png
equivalent_width(spec1d_norm, regions=SpectralRegion(line_min, line_max))
\[-\infty \; \mathrm{\mu m}\]

Find the best-fitting template#

It needs a list of templates and the redshift of the observed galaxy. For the templates, I am using a set of model SEDs generated with Bruzual & Charlot stellar population models, emission lines, and dust attenuation as described in Pacifici et al. (2012).

Developer note
Maybe there is a way to speed this up (maybe using astropy model_sets)? This fit is run with 100 models, but ideally, if we want to extract physical parameters from this, we would need at least 10,000 models. A dictionary structure with meaningful keys (which can be, e.g., tuples of the relevant physical parameters) could be better than a list? It could make later analysis much clearer than having to map from the list indices back to the relevant parameters.

templatedir = './mos_spectroscopy/templates/'
# Redshift taken from the Specviz analysis
zz = 2.0256

f_lamb_units = u.erg / u.s / (u.cm**2) / u.AA

templatelist = []
# Run on 30 out of 100 for speed
for i in range(1, 30):
    template_file = "{0}{1:05d}.dat".format(templatedir, i)
    template = ascii.read(template_file)
    temp1d = Spectrum1D(spectral_axis=(template['col1']/1E4)*u.um, flux=template['col2']*f_lamb_units)
    templatelist.append(temp1d)
# Change the units of the observed spectrum to match the template
spec1d_masked_flamb = spec1d_masked.new_flux_unit(f_lamb_units)
# The new_flux_unit function does not change the uncertainty and specviz complains that there is a mismatch
# so we re-add the uncertainty like we did a few cells above
spec1d_masked_flamb_unc = Spectrum1D(spectral_axis=spec1d_masked_flamb.spectral_axis,
                                     flux=spec1d_masked_flamb.flux,
                                     uncertainty=StdDevUncertainty((np.zeros(len(spec1d_masked_flamb.flux)) + 1E-20) * spec1d_masked_flamb.unit))
WARNING: AstropyDeprecationWarning: The new_flux_unit function is deprecated and may be removed in a future version.
        Use with_flux_unit instead. [warnings]
WARNING:astroquery:AstropyDeprecationWarning: The new_flux_unit function is deprecated and may be removed in a future version.
        Use with_flux_unit instead.
# Take a look at the observed spectrum and one of the templates at the correct redshift
mean_obs = np.mean(spec1d_masked_flamb_unc.flux)
mean_temp = np.mean(templatelist[0].flux)
temp_for_plot = Spectrum1D(spectral_axis=templatelist[0].spectral_axis * (1.+zz),
                           flux=templatelist[0].flux*mean_obs/mean_temp)

plt.figure(figsize=[10, 6])
plt.plot(spec1d_masked_flamb_unc.spectral_axis, spec1d_masked_flamb_unc.flux, label='data')
plt.plot(temp_for_plot.spectral_axis, temp_for_plot.flux, label='model', alpha=0.6)
plt.xlabel("wavelength ({:latex})".format(spec1d_masked_flamb_unc.spectral_axis.unit))
plt.ylabel("flux (normalized)")
plt.xlim(1.1, 1.7)
plt.ylim(0, 2e-18)
plt.legend()
plt.title("Observed spectrum compared to one template at correct redshift")
plt.show()
../../../_images/24c445176d2a81d5f220c68ad10b18bbd813ac3c8bce94db46ffd234219e8a14.png
tm_results = template_comparison.template_match(observed_spectrum=spec1d_masked_flamb_unc, 
                                                spectral_templates=templatelist, 
                                                resample_method="flux_conserving", 
                                                redshift=zz)
tm_results[0]
<Spectrum1D(flux=[0.0 ... 6.044037867431478e-22] erg / (Angstrom s cm2) (shape=(11283,), mean=0.00000 erg / (Angstrom s cm2)); spectral_axis=<SpectralAxis 
   (observer to target:
      radial_velocity=240744.80805001882 km / s
      redshift=2.0256)
  [ 0.03046779  0.0307401   0.0310124  ... 29.34832    29.469344
 29.590368  ] um> (length=11283))>
plt.figure(figsize=[10, 6])
plt.plot(spec1d_masked_flamb_unc.spectral_axis, spec1d_masked_flamb_unc.flux, label="data")
plt.plot(tm_results[0].spectral_axis, tm_results[0].flux, color='r', alpha=0.5, label='model')
plt.xlim(1.0, 1.7)
plt.ylim(0, 5e-19)
plt.xlabel("wavelength ({:latex})".format(spec1d_masked_flamb_unc.spectral_axis.unit))
plt.ylabel("flux ({:latex})".format(spec1d_masked_flamb_unc.flux.unit))
plt.legend()
plt.title("Observed spectrum and best-fitting model template")
plt.show()
../../../_images/64cdd1139f40b9beec7e6a1b6c38dc3dacbb0381206bfd3083073f4c4969df2c.png

New instance of Specviz with spectrum and template#

Passing the spectra with different (but compatible) units. Specviz adopts the first and converts the second spectrum appropriately.

specviz_2 = Specviz()
specviz_2.load_data(spec1d_masked, data_label='observed') # This is in MJy
specviz_2.load_data(tm_results[0], data_label='model') # This is in erg/(s cm^2 A)
specviz_2.show()
Space Telescope Logo

Notebook created by Camilla Pacifici (cpacifici@stsci.edu)
Updated on August 28, 2024