Redshift and Template Fitting

This notebook covers basic examples on how a user can measure the redshift of a source using the visualization tool Jdaviz or programmatically with Specutils.

Use case: measure the redshift of a galaxy from its spectrum using 2 different methods.
Data: JWST/NIRSpec spectrum from program 2736.
Tools: jdaviz, specutils.
Cross-instrument: NIRISS, NIRCam.

Content

• Resources and documentation

• Installation

• Imports

• Fetch the example data

• “By eye” redshift measurement with Specviz

• Redshift measurement with cross-correlation method

  • Get a template and prepare it for use

  • Subtract the continuum from the observed spectrum

  • Clean up the spectrum

  • Run the cross correlation function

Author: Camilla Pacifici (cpacifici@stsci.edu)
Updated: November 18, 2024

Resources and documentation

This notebook uses functionality from Specutils and Jdaviz. Developers at the Space Telescope Science Institute are available to answer questions and resolve problems through the JWST Help Desk. If you wish to send feedback or report problems, you can also submit an issue directly on Github, both for Specutils and for Jdaviz.

Installation

This notebook was extracted from the JWebbinar material. To run this notebook, you will need to create an environment that includes the jdaviz package with the following instructions.
conda create -n jdaviz python=3.11
conda activate jdaviz
from the latest release
pip install jdaviz
or from git
pip install git+https://github.com/spacetelescope/jdaviz.git

Imports

# general os
import tempfile

# numpy
import numpy as np

# astroquery
from astroquery.mast import Observations

# specviz
import jdaviz # this is needed to get the version number later
from jdaviz import Specviz

# astropy
import astropy # again for the version number
import astropy.units as u
from astropy.io import ascii, fits
from astropy.utils.data import download_file
from astropy.modeling.models import Linear1D
from astropy.nddata import StdDevUncertainty

# specutils
import specutils # again for the version number
from specutils import Spectrum1D, SpectralRegion
from specutils.fitting import fit_generic_continuum
from specutils.analysis import correlation
from specutils.manipulation import extract_region

# glue
from glue.core.roi import XRangeROI

# matplotlib
from matplotlib import pyplot as plt

# display
from IPython.display import display, HTML

# customization of matplotlib style
plt.rcParams["figure.figsize"] = (10, 5)
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',
          'figure.figsize': (8, 6)}
plt.rcParams.update(params)
plt.rcParams.update({'figure.max_open_warning': 0})
# This ensures that our notebook is using the full width of the browser
display(HTML("<style>.container { width:100% !important; }</style>"))

Versions:

print("jdaviz:", jdaviz.__version__)
print("astropy:", astropy.__version__)
print("specutils:", specutils.__version__)

jdaviz: 4.0.0
astropy: 6.1.6
specutils: 1.19.0

Fetch the example data

Here we download a spectrum from the Early Release Observation data program 2736 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).

# Select a specific directory on your machine or a temporary directory
data_dir = tempfile.gettempdir()
# Get the file from MAST
fn = "jw02736-o007_s000009239_nirspec_f170lp-g235m_x1d.fits"
result = Observations.download_file(f"mast:JWST/product/{fn}", local_path=f'{data_dir}/{fn}')

fn_template = download_file('https://stsci.box.com/shared/static/3rkurzwl0l79j70ddemxafhpln7ljle7.dat', cache=True)

Downloading URL https://mast.stsci.edu/api/v0.1/Download/file?uri=mast:JWST/product/jw02736-o007_s000009239_nirspec_f170lp-g235m_x1d.fits to /tmp/jw02736-o007_s000009239_nirspec_f170lp-g235m_x1d.fits ...

 [Done]

Jdaviz will read by default the surface brightness extension (in MJy/sr), but I prefer the flux extension (in Jy). I create a Spectrum1D object reading the file myself.

hdu = fits.open(f'{data_dir}/{fn}')
wave = hdu[1].data['WAVELENGTH'] * u.Unit(hdu[1].header['TUNIT1'])
flux = hdu[1].data['FLUX'] * u.Unit(hdu[1].header['TUNIT2'])
# std = hdu[1].data['FLUX_ERROR'] * u.Unit(hdu[1].header['TUNIT3']) # These are all nan. Define an artificial uncertainty

fluxobs = Spectrum1D(spectral_axis=wave,
                     flux=flux,
                     uncertainty=StdDevUncertainty(0.1*flux))
fluxobs

<Spectrum1D(flux=[nan ... nan] Jy (shape=(1300,), mean=0.00000 Jy); spectral_axis=<SpectralAxis [1.64721847 1.64828265 1.64934671 ... 3.1772368  3.17829227 3.17934775] um> (length=1300); uncertainty=StdDevUncertainty)>

“By eye” redshift measurement with Specviz

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 in the Line List plugin. But first, let us open the spectrum in Specviz.

# Call the app
viz = Specviz()
viz.show()

# Load spectrum
# viz.load_data(f'{data_dir}/{fn}', data_label="NIRSpec") # To load the surface brightness directly from file
viz.load_data(fluxobs, data_label='NIRSpec')

Now we need to:

• Open the “line list” plugin

• Choose pre-loaded lines or add custom lines (the lines will not show in the viewer because they are plotted at restframe)

  • Hint: select Ha and the surrounding NII lines

• Input a guess redshift

• Move the slider to get a better match

• Use the zoom tool to get an even better match

Redshift measurement with cross-correlation method

It is very common in astronomy to measure a redshift using a cross-correlation algorithm. IRAF uses this methodology in its xcsao task. Here, we use the Specutils 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:

• Get a template spectrum for the correlation

• Subtract the continuum from both the template and the observed spectrum

• Make sure the spectra have some overlap in wavelength

Get a template and prepare it for use

The template is used for cross correlation, so it can be renormalized for convenience. The units have to match the units of the observed spectrum. We can do the continuum subtraction in erg/(s cm^2 A) since the continuum is close to linear and then run it by Jdaviz to get the appropriate conversion.

# Define unit
spec_unit = u.erg / u.s / u.cm**2 / u.AA

# Read spectrum with the ascii function
template = ascii.read(fn_template)
# Create Spectrum1D object
template = Spectrum1D(spectral_axis=template['col1']/1E4*u.um, 
                      flux=(template['col2']/1E24)*spec_unit) # Normalize to something close to the observed spectrum

# Cut to useful range - template and obs MUST overlap, so we go to 2.4um.
use_tmp = (template.spectral_axis.value > 0.45) & (template.spectral_axis.value < 2.4)
template_cut = Spectrum1D(spectral_axis=template.spectral_axis[use_tmp], flux=template.flux[use_tmp])

# Look at spectrum
plt.figure(figsize=[10, 6])
plt.plot(template_cut.spectral_axis, template_cut.flux)
plt.xlabel("wavelength ({:latex})".format(template_cut.spectral_axis.unit))
plt.ylabel("flux ({:latex})".format(template_cut.flux.unit))
plt.title("Template")
plt.show()

This diagram shows the template spectrum (flux vs wavelength) out to 2.4 micron to allow for some wavelength overlap with the observed spectrum.

# Subtract continuum
mask_temp = ((template_cut.spectral_axis.value > 0.70) & (template_cut.spectral_axis.value < 2.40))

template_forcont = Spectrum1D(spectral_axis=template_cut.spectral_axis[mask_temp], flux=template_cut.flux[mask_temp])

# Use fit_generic_continuum
fit_temp = fit_generic_continuum(template_forcont, model=Linear1D())
cont_temp = fit_temp(template_cut.spectral_axis)
template_sub = template_cut - cont_temp


plt.figure(figsize=[10, 6])
plt.plot(template_cut.spectral_axis, template_cut.flux)
plt.plot(template_cut.spectral_axis, cont_temp)
plt.xlabel("wavelength ({:latex})".format(template_sub.spectral_axis.unit))
plt.ylabel("flux ({:latex})".format(template_sub.flux.unit))
plt.title("Plot template and continuum")
plt.show()

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.

The diagram shows the template spectrum and the best-fitting continuum.

# Print Spectrum1D object
print(template_sub)

Spectrum1D (length=6069)
Flux=[ 3.48301592e-16  3.48211648e-16  3.48231704e-16 ... -1.14961276e-15
      -1.15740032e-15 -1.16516288e-15] erg / (Angstrom s cm2),  mean=0.00000 erg / (Angstrom s cm2)
Spectral Axis=[0.45008 0.45017 0.45026 ... 2.375   2.385   2.395  ] um,  mean=0.75618 um

# Look at spectrum
plt.figure(figsize=[10, 6])
plt.plot(template_sub.spectral_axis, template_sub.flux)
plt.xlabel("wavelength ({:latex})".format(template_sub.spectral_axis.unit))
plt.ylabel("flux ({:latex})".format(template_sub.flux.unit))
plt.title("Continuum subtracted template")
plt.show()

This diagram shows the template spectrum (flux vs wavelength) after subtracting the fitted continuum.

Subtract the continuum from the observed spectrum

We can use a different approach and do it with SpectralRegion here. We also need to include an uncertianty to the observed spectrum, if it is not included.

# Define Spectral Region
region = SpectralRegion(2.0*u.um, 3.0*u.um)
# Extract region
spec1d_cont = extract_region(fluxobs, region)
# Run fitting function
fit_obs = fit_generic_continuum(spec1d_cont, model=Linear1D(5))
# Apply to spectral axis
cont_obs = fit_obs(fluxobs.spectral_axis)

# Subtract continuum
spec1d_sub = fluxobs - cont_obs

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.

# Look at spectrum
plt.figure(figsize=[10, 6])
plt.plot(spec1d_sub.spectral_axis, spec1d_sub.flux)
plt.xlabel("wavelength ({:latex})".format(spec1d_sub.spectral_axis.unit))
plt.ylabel("flux ({:latex})".format(spec1d_sub.flux.unit))
plt.title("Continuum subtracted observed spectrum")
plt.show()

This diagram shows the observed spectrum (flux vs wavelength) after subtracting the fitted continuum.

Clean up the spectrum

It is best to remove artifacts that can look like emission/absorption lines and remove big gaps. The selection of the clean part of the spectrum can be done using the GUI. If this is not performed manually, the following cell takes care of it programmatically.

viz2 = Specviz()
viz2.load_data(spec1d_sub, data_label='spectrum continuum subtracted')
viz2.show()

# Create a subset in the area of interest if it has not been created manually
try:
    region1 = viz2.get_spectra(data_label='spectrum continuum subtracted', subset_to_apply='Subset 1')
    print(region1)
    region1_exists = True
except Exception:
    print("There are no subsets selected.")
    region1_exists = False
    
# Spectral region for masking artifacts
if region1_exists is False:
    sv = viz2.app.get_viewer('spectrum-viewer')
    sv.toolbar_active_subset.selected = []
    sv.apply_roi(XRangeROI(2.24, 3.13))  

There are no subsets selected.

# Get spectrum out with mask
spec1d_region = viz2.get_spectral_regions()
spec1d_masked = extract_region(spec1d_sub, spec1d_region['Subset 1'], return_single_spectrum=True)

We visualize the observed and template continuum-subtracted spectra in a new instance of Specviz. Hit the Home button to see the entire wavelength range. The template spectrum will change unit to match the observed spectrum.

viz3 = Specviz()
viz3.load_data(spec1d_masked, data_label='observation')
viz3.load_data(template_sub, data_label='template')
viz3.show()

template_newunit = viz3.get_data('template', use_display_units=True)
template_newunit

<Spectrum1D(flux=[3.48301592239088e-16 ... -1.1651628832986837e-15] erg / (Angstrom s cm2) (shape=(6069,), mean=0.00000 erg / (Angstrom s cm2)); spectral_axis=<SpectralAxis [0.45008 0.45017 0.45026 ... 2.375   2.385   2.395  ] um> (length=6069))>

Run the cross correlation function

# Call the function
corr, lag = correlation.template_correlate(spec1d_masked, template_newunit, lag_units=u.one)

# Plot the correlation
plt.plot(lag, corr)
plt.xlabel("Redshift")
plt.ylabel("Correlation")
plt.show()

This diagram shows the correlation value vs the redshift. The spike (around redshift of 2.5) indicates the value where the observed and template spectra correlate best.

# Redshift based on maximum
index_peak = np.argmax(corr)
z = lag[index_peak]

print("Redshift from peak maximum: ", z)

Redshift from peak maximum:  2.4638739503966023

# Redshifted template_sub
template_sub_z = Spectrum1D(spectral_axis=template_sub.spectral_axis * (1. + z),
                            flux=template_sub.flux)

# Visualize the redshifted template and the observed spectrum
viz4 = Specviz()
viz4.load_data(spec1d_masked, data_label='Observed spectrum')
viz4.load_data(template_sub_z, data_label='Redshifted best template')
viz4.show()