<jupyter_start><jupyter_text>MOS Spectroscopy of Typical Extragalactic Field**Use case:** emission-line measurements and template matching on 1D spectra.**Data:** [LEGA-C](https://www2.mpia-hd.mpg.de/home/legac/) spectra and galaxy template spectra; optical rest-frame.**Tools:**  specutils, astropy, matplotlib, pandas.**Cross-intrument:** NIRSpec **Documentation:** This notebook is part of a STScI's larger [post-pipeline Data Analysis Tools Ecosystem](https://jwst-docs.stsci.edu/jwst-post-pipeline-data-analysis). IntroductionThis notebook will perform a seris of spectroscopic analyses on multiple spectra, including smoothing, continuum fitting and subtraction, line identification, centroiding and flux measurements, gaussian fitting, equivalent widths, and template fitting.**Note:** This notebook is intended to ultimately be compatible with the final data products (1D and 2D spectra) from the JWST pipeline. These data products are not available yet, so the notebook uses LEGA-C data (van der Wel et al. 2016, Straatmann et al. 2018) for now.LEGA-C is a galaxy survey of about 3000 galaxies at z~0.8 and M* > 10^10 M_sun in the COSMOS field. The spectra sample the rest-frame optical between ~3000A and 5000A at high resolution and very high signal-to-noise ratio. More information about the survey can be found here: http://www.mpia.de/home/legac/ First, set the environment with astroconda and specutils.`conda create -n astroutils stsci``source activate astroutils`  Imports<jupyter_code>#general os
import os
import zipfile
import urllib.request

#general plotting
from matplotlib import pyplot as plt

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})

#table/math handling
import pandas as pd
import numpy as np
np.seterr(all='ignore')  # hides irrelevant warnings about divide-by-zero, etc

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

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

quantity_support();  # auto-recognizes units on matplotlib plots<jupyter_output><empty_output><jupyter_text>Check versions. Should be:Pandas:  1.0.1Numpy:  1.18.1Astropy:  4.0Specutils:  0.7<jupyter_code>print("Pandas: ",pd.__version__)
print("Numpy: ",np.__version__)
print("Astropy: ",astropy.__version__)
print("Specutils: ",specutils.__version__)<jupyter_output><empty_output><jupyter_text>Set path to data and download from box link<jupyter_code>boxlink = 'https://data.science.stsci.edu/redirect/JWST/jwst-data_analysis_tools/mos_spectroscopy/mos_spectroscopy.zip'
boxfile = './mos_spectroscopy.zip'
urllib.request.urlretrieve(boxlink, boxfile)

zf = zipfile.ZipFile(boxfile, 'r')
zf.extractall()

observedfiles = './mos_spectroscopy/observed/'<jupyter_output><empty_output><jupyter_text>Choose one galaxy<jupyter_code>file1d = observedfiles + 'legac_M1_v3.7_spec1d_130902.fits'
file1dwht = observedfiles + 'legac_M1_v3.7_wht1d_130902.fits'
file2d = observedfiles + 'legac_M1_v3.7_spec2d_130902.fits'<jupyter_output><empty_output><jupyter_text>Inspect its 2D spectrum<jupyter_code>hdu2d = fits.open(file2d)
hdu2d.info()
hdu2d[0].header<jupyter_output><empty_output><jupyter_text>Developer noteI would appreciate the interactive tools here to zoom and pan through the 2D spectrum. Hoovering to know the precise wavelength of a feature would also be very useful. With that, the interactive tool could show automatically the calibration in wavelength reading it from the header.<jupyter_code>plt.figure(figsize=[20,20])
plt.imshow(hdu2d[0].data)
plt.xlim(2000,3000) #spec is very big, plot just a bit<jupyter_output><empty_output><jupyter_text>Now work with 1D spectrum Calibrate (in wavelength), inspect, and write in Spectrum1D object<jupyter_code>hdu1d = fits.open(file1d)
hdu1dwht = fits.open(file1dwht)
hdu1d.info()
hdu1d[0].header<jupyter_output><empty_output><jupyter_text>If i want to do it in Pandas.<jupyter_code>flux = hdu1d[0].data.byteswap().newbyteorder()
wht = hdu1dwht[0].data.byteswap().newbyteorder()
unc = 1./ np.sqrt(wht)
wave = np.arange(flux.shape[0])*hdu1d[0].header['CD1_1'] + hdu1d[0].header['CRVAL1']

##for certain functions, I need to cut the spectrum where the weight is 0.
d = {'wavelength':wave, 'flux':flux, 'weight':wht, 'uncertainty':unc}
dataspec = pd.DataFrame(data=d)

##define subset where wht>0
dataspec_sub = dataspec[dataspec['weight'] > 0.].reset_index(drop=True)
dataspec_sub.head()<jupyter_output><empty_output><jupyter_text>If I want to do it in astropy Tables.<jupyter_code>flux = hdu1d[0].data
wht = hdu1dwht[0].data
unc = 1./ np.sqrt(wht)
wave = np.arange(flux.shape[0])*hdu1d[0].header['CD1_1'] + hdu1d[0].header['CRVAL1']

spec_unit = u.Unit('10^-19 erg s^-1 cm^-2 angstrom^-1')
dataspec = QTable([wave*u.angstrom, flux*spec_unit, wht, unc*spec_unit], 
                   names=('wavelength','flux','weight','uncertainty'))
dataspec_sub = dataspec[dataspec['weight']>0.]
dataspec_sub
plt.figure(figsize=[10,6])
plt.fill_between(dataspec_sub['wavelength'],
                 dataspec_sub['flux']+dataspec_sub['uncertainty'],
                 dataspec_sub['flux']-dataspec_sub['uncertainty'],
                 color='gray',label="uncertainty")
plt.plot(dataspec_sub['wavelength'],dataspec_sub['flux'], label="data")
plt.xlabel("wavelength ({:latex})".format(dataspec_sub['wavelength'].unit))
plt.ylabel("flux ({:latex})".format(dataspec_sub['flux'].unit))
plt.title("Observed spectrum")
plt.legend()
plt.show()<jupyter_output><empty_output><jupyter_text>Go with specutils<jupyter_code>#write Spectrum1D object
spec1d = Spectrum1D(spectral_axis=dataspec_sub['wavelength'], 
                    flux=dataspec_sub['flux'], 
                    uncertainty=StdDevUncertainty(dataspec_sub['uncertainty']))<jupyter_output><empty_output><jupyter_text>Developer noteFor supported datasets (like final JWST data products), this will be as simple as:```spec1d = Spectrum1d.read('datafile.fits')```  Developer noteImplemented but not yet released: `snr_threshold`, which will allow cutting the spectrum using that function.  Smooth to better inspect the features  Developer noteThe uncertainty is not carried over, but would be useful. Same comment on the interactive tool as before: it would be very useful to be able to zoom, pan, hoover, etc. on the spectrum.<jupyter_code>spec1d_gsmooth = gaussian_smooth(spec1d, stddev=5)
plt.figure(figsize=[10,6])
plt.plot(spec1d_gsmooth.spectral_axis,spec1d_gsmooth.flux)
plt.xlabel("wavelength ({:latex})".format(spec1d_gsmooth.spectral_axis.unit))
plt.ylabel("flux ({:latex})".format(spec1d_gsmooth.flux.unit))
plt.title("Smoothed observed spectrum")
plt.show()<jupyter_output><empty_output><jupyter_text>Back to the non-smoothed spectrum to find linesDocumentation says I need a continuum subtracted spectrum. So fit continuum first<jupyter_code>cont_spec1d = fit_generic_continuum(spec1d)
cont_fit = cont_spec1d(spec1d.spectral_axis)
plt.figure(figsize=[10,6])
plt.plot(spec1d.spectral_axis, spec1d.flux, label="data")
plt.plot(spec1d.spectral_axis, cont_fit, label="modeled continuum")
plt.xlabel("wavelength ({:latex})".format(spec1d.spectral_axis.unit))
plt.ylabel("flux ({:latex})".format(spec1d.flux.unit))
plt.legend()
plt.title("Observed spectrum and fitted continuum")
plt.show()

plt.figure(figsize=[10,6])
plt.plot(spec1d.spectral_axis, spec1d.uncertainty.array, label="data")
plt.xlabel("wavelength ({:latex})".format(spec1d.spectral_axis.unit))
plt.ylabel("uncertainty ({:latex})".format(spec1d.uncertainty.unit))
plt.legend()
plt.title("Uncertianty of observed spectrum")
plt.show()<jupyter_output><empty_output><jupyter_text>Creating the continuum-subtracted spectrumSpecutils will figure out what to do with the uncertainty!<jupyter_code>spec1d_sub = spec1d - cont_fit
spec1d_sub
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()<jupyter_output><empty_output><jupyter_text>Now look for the lines<jupyter_code>lines = find_lines_derivative(spec1d_sub, flux_threshold=50)
lines
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/absorption lines')
for line in lines:
    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.flux.unit))
plt.legend()
plt.title("Continuum-subtracted spectrum and marked lines using find_lines_derivative")
plt.show()<jupyter_output><empty_output><jupyter_text>This works for cases where you understand the thresholds well, but doesn't automate as well with noisy spectra.  Works better with find_lines_threshold<jupyter_code>lines = find_lines_threshold(spec1d_sub, noise_factor=6)
lines<jupyter_output><empty_output><jupyter_text>Plot lines on the spectrum.<jupyter_code>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/absorption lines')
for line in lines:
    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()<jupyter_output><empty_output><jupyter_text>Developer noteIt would be useful to have a tool to cycle through the lines, show a zoom of the spectrum, and inspect how good the line identification is. For now I do it by hand on a single line.<jupyter_code>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:
    plt.axvline(line['line_center'].value, alpha=0.5, color='red')
plt.xlim(6500,6600)
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()<jupyter_output><empty_output><jupyter_text>Measure line centroids and fluxesThese too need spectra continuum subtracted.<jupyter_code>#example with just one line
centroid(spec1d_sub, SpectralRegion(6540*u.AA, 6580*u.AA))
sline = centroid(spec1d_sub, SpectralRegion(6540*u.AA, 6580*u.AA))

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(6500,6600)
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()
line_flux(spec1d_sub, SpectralRegion(6540*u.AA, 6570*u.AA))<jupyter_output><empty_output><jupyter_text>Fit the line with a Gaussian Developer noteFitting lines with Gaussians is so common that it might make sense to have a one-line to iterate this over the list of lines, with sensible initialization of parameters and fit ranges for each line so that the fits generally work.Also, the interactive tool here could allow one to select the continuum region and position of the line by clicking on the spectrum.<jupyter_code>g_init = models.Gaussian1D(amplitude= 3 * 1e-19 * u.erg / u.s / u.cm**2 / u.AA, mean=6554*u.AA, stddev=2.*u.AA)
g_fit = fit_lines(spec1d_sub, g_init)
spec1d_fit = g_fit(spec1d_sub.wavelength)
g_fit
#calculate the velocity dispertion from the stddev
vel = ((3.71890256/6554.415) * const.c.to('km/s').value)*u.km/u.s
print(vel)
plt.figure(figsize=[10,6])
plt.plot(spec1d_sub.wavelength,spec1d_sub.flux,label='data')
plt.plot(spec1d_sub.wavelength,spec1d_fit,color='darkorange',label='Gaussian fit')
plt.xlim(6500,6600)
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()<jupyter_output><empty_output><jupyter_text>Measure the equivalent width of the linesThis needs the spectrum continuum normalized.<jupyter_code>spec1d_norm = spec1d / 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(6500,6600)
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(6500,6600)
plt.legend()
plt.title("Uncertainty of continuum-normalized spectrum, zoomed on [OII]")
plt.show()
equivalent_width(spec1d_norm, regions=SpectralRegion(6540*u.AA, 6570*u.AA))<jupyter_output><empty_output><jupyter_text>Find the best-fitting templateIt 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 noteMaybe 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.<jupyter_code>templatedir = './mos_spectroscopy/templates/'
zz = 0.758

templatelist = []
for i in range (1, 101):
    template_file = "{0}{1:05d}.dat".format(templatedir,i)
    template = ascii.read(template_file)
    temp1d = Spectrum1D(spectral_axis=template['col1']*u.AA,flux=template['col2']*u.erg/u.s/u.AA)
    templatelist.append(temp1d)
    
    
tm_results = template_comparison.template_match(observed_spectrum=spec1d, spectral_templates=templatelist, resample_method="flux_conserving", redshift=zz)
tm_results[0]
plt.figure(figsize=[10, 6])
plt.plot(spec1d.wavelength, spec1d.flux, label="data")
plt.plot(tm_results[0].wavelength, tm_results[0].flux,color='r',alpha=0.5,label='model')
plt.xlim(6000, 9000)
plt.xlabel("wavelength ({:latex})".format(spec1d_norm.spectral_axis.unit))
plt.ylabel("flux ({:latex})".format(spec1d_sub.flux.unit))
plt.legend()
plt.title("Observed spectrum and best-fitting model template")
plt.show()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1298170731707317, 'fraction_numerical': 0.027134146341463414, 'mean_word_length': 7.171898355754858, 'pattern_counts': {'":': 0, '<': 85, '<?xml version=': 0, '>': 89, 'https://': 3, 'lorem ipsum': 0, 'www.': 2, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 8, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '21', 'dir': 'jupyter-structured-clean-dedup', 'id': '601733', 'license': 'Apache-2.0', 'n_tokens_mistral': 6119, 'n_tokens_neox': 5613, 'n_words': 1238, 'path': '_notebooks/mos-spectroscopy/MOSspec_sv06_revised.ipynb', 'repo_name': 'TheRealZoidberg/fastpagesTest', 'size': '16400'}