WFSS Spectra Part 3: Emission Line Map

WFSS Spectra Part 3: Emission Line Map#

Use case: optimal extraction of grism spectra; redshift measurement; emission-line maps. Simplified version of JDox Science Use Case # 33.
Data: JWST simulated NIRISS images from MIRAGE, run through the JWST calibration pipeline; galaxy cluster.
Tools: specutils, astropy, pandas, emcee, lmfit, corner, h5py.
Cross-intrument: NIRSpec
Documentation: This notebook is part of a STScI’s larger post-pipeline Data Analysis Tools Ecosystem.

Introduction#

This notebook is 4 of 4 in a set focusing on NIRISS WFSS data: 1. 1D optimal extraction since the JWST pipeline only provides a box extraction. Optimal extraction improves S/N of spectra for faint sources. 2. Combine and normalize 1D spectra. 3. Cross correlate galaxy with template to get redshift. 4. Spatially resolved emission line map.

This notebook obtains an emission line map from the 2D-cutout of NIRISS WFSS spectra. This notebook will require the 2D-cutout of NIRISS WFSS spectrum (e.g., l3_nis_f150w_G150C_s00004_cal.fits; Pipeline spec3 product).

Note: For this example, the notebook uses only a 2D rectified spectra at one dither position, directly obtained from the current versionn of the pipeline (Build 7.5), but one could start with spectra stacked with all dither positions too, to improve S/Ns.

%matplotlib inline

import os

import numpy as np

from astropy.io import fits
import astropy.units as u
from astropy.modeling import models
from astropy.modeling.polynomial import Chebyshev1D

from specutils.spectra import Spectrum1D
from specutils.fitting import fit_lines
from specutils.fitting import continuum

from astropy import __version__ as astropy_version
print('astropy', astropy_version)

import matplotlib.pyplot as plt
import matplotlib as mpl

mpl.rcParams['savefig.dpi'] = 80
mpl.rcParams['figure.dpi'] = 80

0. Download notebook 01 products#

These can be also obtained by running the notebooks.

if not os.path.exists('./output'):
    import zipfile
    import urllib.request
    boxlink = 'https://data.science.stsci.edu/redirect/JWST/jwst-data_analysis_tools/NIRISS_lensing_cluster/output.zip'
    boxfile = './output.zip'
    urllib.request.urlretrieve(boxlink, boxfile)
    zf = zipfile.ZipFile(boxfile, 'r')
    zf.extractall()
else:
    print('Already exists')

# Download files, if not exists yet.
if not os.path.exists('./pipeline_products'):
    import zipfile
    import urllib.request
    boxlink = 'https://data.science.stsci.edu/redirect/JWST/jwst-data_analysis_tools/NIRISS_lensing_cluster/pipeline_products.zip'
    boxfile = './pipeline_products.zip'
    urllib.request.urlretrieve(boxlink, boxfile)
    zf = zipfile.ZipFile(boxfile, 'r')
    zf.extractall()

# Which target?
DIR_OUT = './output/'
filt = 'f200w'

grism = 'G150C'
# grism = 'G150R'

id = '00004'

1.Get Halpha emission line map;#

# Open two-dimensional file;
# can be downloaded in No.01a.
DIR_DATA = './pipeline_products/'

# Zero-indexed number for dither --- the test data has two dither positions.
ndither = 0

# File
file_2d = f'{DIR_DATA}l3_nis_{filt}_{grism}_s{id}_cal.fits'
hdu_2d = fits.open(file_2d)

# Align grism direction
#   - x-direction = Dispersion (wavelength) direction.
#   - y-direction = Cross-dispersion.
# in this notebook.
    
if grism == 'G150C':
    # If spectrum is horizontal;
    data_2d = hdu_2d[ndither*7+1].data
    dq_2d = hdu_2d[ndither*7+2].data
    err_2d = hdu_2d[ndither*7+3].data
    wave_2d = hdu_2d[ndither*7+4].data
else:
    data_2d = rotate(hdu_2d[ndither*7+1].data, 90)
    dq_2d = rotate(hdu_2d[ndither*7+2].data, 90)
    err_2d = rotate(hdu_2d[ndither*7+3].data, 90)
    wave_2d = rotate(hdu_2d[ndither*7+4].data, 90)

plt.imshow(data_2d[:, :])

Subtract continuum at each y-position:#

There are two options, depending on the brightness of the source;

  • 1.Subtract continuum estimated at each position of cross-dispersion column. This works for a bright source.

  • 2.Subtract a single continuum (“master continuum”), assuming the shape is same across cross-dispersion direction as the one from a 1d extraction. This is required for a faint source.

Since we clearly see continuum across the source, here we demonstrate the option 1.

# At which y position do we want to see spectrum?
# as an example;
yy = 10

spec_unit = u.dimensionless_unscaled

mask_line = ((wave_2d[yy, :] > 1.75) & (wave_2d[yy, :] < 1.97)) | ((wave_2d[yy, :] > 2.08) & (wave_2d[yy, :] < 2.23))

obs = Spectrum1D(spectral_axis=wave_2d[yy, :][mask_line]*u.um, flux=data_2d[yy, :][mask_line]*spec_unit)
cont = continuum.fit_generic_continuum(obs, model=Chebyshev1D(7))

plt.plot(wave_2d[yy, :], data_2d[yy, :], color='r', label=f'Observed at y={yy}')
plt.plot(wave_2d[yy, :][mask_line]*u.um, cont(wave_2d[yy, :][mask_line]*u.um), color='b', label='Fitted continuum')

plt.legend(loc=0)
plt.xlim(1.6, 2.3)

cont(wave_2d[yy, :][mask_line]*u.um).value

# Repeat this along y-axis;
flux_cont2d = data_2d[:, :] * 0
for yy in range(len(data_2d[:, 0])):

    mask_line = ((wave_2d[yy, :] > 1.75) & (wave_2d[yy, :] < 1.97)) | ((wave_2d[yy, :] > 2.08) & (wave_2d[yy, :] < 2.23))

    obs = Spectrum1D(spectral_axis=wave_2d[yy, :][mask_line]*u.um, flux=data_2d[yy, :][mask_line]*spec_unit)
    cont = continuum.fit_generic_continuum(obs, model=Chebyshev1D(7))

    flux_cont2d[yy, :] = cont(wave_2d[yy, :]*u.um).value

plt.imshow(flux_cont2d[:, :])
plt.title('Fitted continuum')

plt.imshow(data_2d[:, :] - flux_cont2d[:, :])
plt.title('Continuum subtracted spectrum')

Extract Halpha emission from a continuum-subtracted spectrum;#

# Cutout Ha line map;
rsq = data_2d.shape[0]
cut_ha = np.zeros((rsq, rsq), 'float32')

zin = 2.1 # From Notebook No.3, cross-correlation;
lamcen = 0.6564 * (1. + zin)

for yy in range(len(data_2d[:, 0])):
    # This has to be done at each y pixel, as wavelength array can be tilted.
    index_lamcen = np.argmin(np.abs(lamcen - wave_2d[yy, :]))
    cut_ha[yy, :] = (data_2d - flux_cont2d)[yy, int(index_lamcen-rsq/2.):int(index_lamcen+rsq/2.)]

plt.imshow(cut_ha)
plt.title('H$\\alpha$ map')

2.Get Hbeta and OIII maps#

This is more challenging, as these lines locate close to each other. Ideally, iteration process will be preferred, but here we use Specutils’ double gaussian component fitting, in a similar way for Ha.

filt = 'f150w'

file_2d = f'{DIR_DATA}l3_nis_{filt}_{grism}_s{id}_cal.fits'
hdu_2d = fits.open(file_2d)

# Align grism direction
#   - x-direction = Dispersion (wavelength) direction.
#   - y-direction = Cross-dispersion.
# in this notebook.
    
if grism == 'G150C':
    # If spectrum is horizontal;
    data_2d = hdu_2d[ndither*7+1].data
    dq_2d = hdu_2d[ndither*7+2].data
    err_2d = hdu_2d[ndither*7+3].data
    wave_2d = hdu_2d[ndither*7+4].data
else:
    data_2d = rotate(hdu_2d[ndither*7+1].data, 90)
    dq_2d = rotate(hdu_2d[ndither*7+2].data, 90)
    err_2d = rotate(hdu_2d[ndither*7+3].data, 90)
    wave_2d = rotate(hdu_2d[ndither*7+4].data, 90)
    
# !! Note that the extracted spectra has flipped wavelength direction !!
plt.imshow(data_2d[:, ::-1])
plt.title(f'{filt}')

In the plot above, you can see Oiii doublet, and Hbeta.#

2a. Get continuum;#

yy = 10 # as an example;

spec_unit = u.dimensionless_unscaled

mask_line = ((wave_2d[yy, :] > 1.35) & (wave_2d[yy, :] < 1.48)) | ((wave_2d[yy, :] > 1.6) & (wave_2d[yy, :] < 1.65))

obs = Spectrum1D(spectral_axis=wave_2d[yy, :][mask_line]*u.um, flux=data_2d[yy, :][mask_line]*spec_unit)
cont = continuum.fit_generic_continuum(obs, model=Chebyshev1D(7))

plt.plot(wave_2d[yy, :], data_2d[yy, :], color='r', label=f'Observed at y={yy}')

plt.plot(wave_2d[yy, :][mask_line]*u.um, cont(wave_2d[yy, :][mask_line]*u.um), color='b', label='Fitted continuum')

plt.legend(loc=0)
plt.xlim(1.2, 1.8)

# Repeat this along y-axis;
flux_cont2d_150 = data_2d[:, :] * 0
for yy in range(len(data_2d[:, 0])):

    mask_line = ((wave_2d[yy, :] > 1.35) & (wave_2d[yy, :] < 1.48)) | ((wave_2d[yy, :] > 1.6) & (wave_2d[yy, :] < 1.65))

    obs = Spectrum1D(spectral_axis=wave_2d[yy, :][mask_line]*u.um, flux=data_2d[yy, :][mask_line]*spec_unit)
    cont = continuum.fit_generic_continuum(obs, model=Chebyshev1D(7))

    flux_cont2d_150[yy, :] = cont(wave_2d[yy, :]*u.um).value

plt.imshow(flux_cont2d_150[:, ::-1])

line_2d = data_2d[:, :] - flux_cont2d_150[:, :]

plt.imshow(line_2d[:, ::-1])

2b. Fit emissionn, with multi-conponent gaussian;#

yy = 10

# Fit the spectrum
con = (1.4 < wave_2d[yy, :]) & (wave_2d[yy, :] < 1.65)

spectrum_cut = Spectrum1D(flux=line_2d[yy, :][con]*spec_unit,
                          spectral_axis=wave_2d[yy, :][con]*u.um)

# !!! Some tweaks may be needed for initial value, to successfuully run the fit;

# For Hb
g1_init = models.Gaussian1D(amplitude=40*spec_unit, mean=1.5*u.um, stddev=0.005*u.um)

# For O3 blue
g2_init = models.Gaussian1D(amplitude=50.*spec_unit, mean=1.535*u.um, stddev=0.002*u.um)

# For O3 red
g3_init = models.Gaussian1D(amplitude=45.*spec_unit, mean=1.55*u.um, stddev=0.001*u.um)

g123_fit = fit_lines(spectrum_cut, g1_init+g2_init+g3_init, window=[0.01*u.um, 0.001*u.um, 0.001*u.um])
y_fit = g123_fit(wave_2d[yy, :]*u.um)

print(g123_fit)

# Plot separately?
plt.plot(wave_2d[yy, :], line_2d[yy, :], marker='.', ls='', color='r', label=f'Observed at y={yy}')
plt.plot(wave_2d[yy, :], y_fit, color='b', label='Fit', zorder=-2, alpha=0.4, lw=5)

y_fit1 = g123_fit[0](wave_2d[yy, :]*u.um)
plt.plot(wave_2d[yy, :], y_fit1, color='g', label='1')

y_fit1 = g123_fit[1](wave_2d[yy, :]*u.um)
plt.plot(wave_2d[yy, :], y_fit1, color='orange', label='2')

y_fit1 = g123_fit[2](wave_2d[yy, :]*u.um)
plt.plot(wave_2d[yy, :], y_fit1, color='purple', label='3')

plt.xlim(1.4, 1.6)

plt.title('Single fit peak')
plt.grid(True)
plt.legend(loc=2)

print(g123_fit[0])
print(g123_fit[1])
print(g123_fit[2])

Fit to the central array looks good. Repeat this along y-axis and get emission line maps, as a same way as for Halpha#

# Cutout Hb, Oiii line maps;
rsq = data_2d.shape[0]
cut_hb = np.zeros((rsq, rsq), 'float32')
cut_o3b = np.zeros((rsq, rsq), 'float32')
cut_o3r = np.zeros((rsq, rsq), 'float32')

zin = 2.1 # Redshift estimate from Notebook No.2, cross-correlation;
lamcen_hb = 0.4862680 * (1. + zin)
lamcen_o3b = 0.4960295 * (1. + zin)
lamcen_o3r = 0.5008240 * (1. + zin)

for yy in range(len(data_2d[:, 0])):
    # Fit the spectrum
    con = (1.4 < wave_2d[yy, :]) & (wave_2d[yy, :] < 1.65)
    spectrum_cut = Spectrum1D(flux=line_2d[yy, :][con]*spec_unit, 
                              spectral_axis=wave_2d[yy, :][con]*u.um)

    # !!! Some tweaks may be needed for initial value, to successfuully run the fit;

    # For Hb
    g1_init = models.Gaussian1D(amplitude=40*spec_unit, mean=1.5*u.um, stddev=0.005*u.um)
    # For O3 blue
    g2_init = models.Gaussian1D(amplitude=50.*spec_unit, mean=1.535*u.um, stddev=0.002*u.um)
    # For O3 red
    g3_init = models.Gaussian1D(amplitude=45.*spec_unit, mean=1.55*u.um, stddev=0.001*u.um)
    g123_fit = fit_lines(spectrum_cut, g1_init+g2_init+g3_init, window=[0.01*u.um, 0.001*u.um, 0.001*u.um])
    y_fit = g123_fit(wave_2d[yy, :]*u.um)

    # This has to be done at each y pixel, as wavelength array can be tilted.
    index_lamcen_hb = np.argmin(np.abs(lamcen_hb - wave_2d[yy, :]))
    cut_hb[yy, :] = g123_fit[0](wave_2d[yy, :]*u.um)[int(index_lamcen_hb-rsq/2.):int(index_lamcen_hb+rsq/2.)]

    index_lamcen_o3b = np.argmin(np.abs(lamcen_o3b - wave_2d[yy, :]))
    cut_o3b[yy, :] = g123_fit[1](wave_2d[yy, :]*u.um)[int(index_lamcen_o3b-rsq/2.):int(index_lamcen_o3b+rsq/2.)]

    index_lamcen_o3r = np.argmin(np.abs(lamcen_o3r - wave_2d[yy, :]))
    cut_o3r[yy, :] = g123_fit[2](wave_2d[yy, :]*u.um)[int(index_lamcen_o3r-rsq/2.):int(index_lamcen_o3r+rsq/2.)]

plt.imshow(cut_hb)
plt.title('H$\\beta$ map')

plt.imshow(cut_o3b)
plt.title('Oiii 4960')

plt.imshow(cut_o3r)
plt.title('Oiii 5008')

Summary;#

As seen above, regions except for the very center do not look right. This is due to failure of multi-component fit, especially for Oiii doublet. To improve the fit, one can either;

  • fit and inspect the fit repeatedly at each y-axis until it converges,

  • use MCMC for more intensive fitting, which also enables to fix the ratio of two Oiii lines by setting up a prior.

  • or reduce the number of components, especially for Oiii doublet at the edge of the source position, where the lines are blended.

Contents