WFSS Spectra Part 1: Combine and Normalize 1D Spectra#

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 2 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 will first combine optimally extracted 1D spectra (from #1) at different dither positions. The combined spectrum will then be normalized by flux estimated from direct images (i.e. broadband photometry).

This notebook will start with an optimally extracted 1D spectrum in the previous notebook, saved in a file with a name similar to “l3_nis_f115w_G150C_s00002_ndither0_1d_opt.fits”. The various one-dimensional spectra from each dither will be combined and saved into a single file, named like “l3_nis_f115w_G150C_s00004_combine_1d_opt.fits”.

Note: The procedure is not intended to combine spectra from the two different orientations (C&R), as each has a different spectral resolution depending on source morphology.

Note: Normalization of grism spectrum to each broadband filter photometry will improve zeropoint calibration, which is critical in continuum fitting over multiple filters.

%matplotlib inline
import os
import numpy as np
from scipy.integrate import simps
from urllib.request import urlretrieve
import tarfile

from astropy.io import fits
import astropy.units as u
from astropy.io import ascii

import astropy
print('astropy', astropy.__version__)
astropy 5.3.4
import matplotlib.pyplot as plt
import matplotlib as mpl

# These gymnastics are needed to make the sizes of the figures
# be the same in both the inline and notebook versions
%config InlineBackend.print_figure_kwargs = {'bbox_inches': None}

mpl.rcParams['savefig.dpi'] = 80
mpl.rcParams['figure.dpi'] = 80
import specutils
from specutils.spectra.spectrum1d import Spectrum1D
from astropy.nddata import StdDevUncertainty

print("Specutils: ", specutils.__version__)
Specutils:  1.12.0

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')
# Which data set?;
DIR_OUT = './output/'
filt = 'f200w'
grism = 'G150C'
id = '00004'
ndither = 2 # Number of dithering positions.

1.Combine spectra from different dithers;#

dithers = np.arange(0, ndither, 1)
for dither in dithers:
    file_1d = f'{DIR_OUT}l3_nis_{filt}_{grism}_s{id}_ndither{dither}_1d_opt.fits'
    fd = fits.open(file_1d)[1].data
    if dither == 0:
        wave = fd['wavelength']
        flux = np.zeros((ndither, len(wave)), 'float')
        flux_err = np.zeros((ndither, len(wave)), 'float')
        flux[dither, :] = fd['flux']
        flux_err[dither, :] = fd['uncertainty']
    else:
        wave_tmp = fd['wavelength']
        flux_tmp = fd['flux']
        flux_err_tmp = fd['uncertainty']
        
        flux[dither, :] = np.interp(wave, wave_tmp, flux_tmp)
        flux_err[dither, :] = np.sqrt(np.interp(wave, wave_tmp, flux_err_tmp[:]**2))
plt.plot(wave, flux[0, :], color='r', label='dither1')
plt.plot(wave, flux[1, :], color='b', label='dither2')
plt.legend(loc=2)
<matplotlib.legend.Legend at 0x7fb3bf84efd0>
../../_images/eae46b2c44a6ccc9c6c2ec669ea09db5e4f5a39e1eed6425241362bf32596300.png

Make weighted-average of the two spectra from different dither positions;#

flux_combine = np.zeros(len(wave), 'float')
flux_err_combine = np.zeros(len(wave), 'float')

for ii in range(len(wave)):
    flux_combine[ii] = np.sum(flux[:, ii] * 1 / flux_err[:, ii]**2) / np.sum(1 / flux_err[:, ii]**2) 
    flux_err_combine[ii] = np.sqrt(np.sum(flux_err[:, ii]**2)) / len(flux_err[:, ii])

plt.plot(wave, flux[0, :], color='gray', label='dither1')
plt.plot(wave, flux[1, :], color='gray', label='dither2')
plt.plot(wave, flux_combine[:], color='r', label='Combined')
plt.legend(loc=2)
<matplotlib.legend.Legend at 0x7fb3bf6f4410>
../../_images/cf2f42eec739ca18bcd914490fbc6bbc08e7305264a05254dd066f720cd58872.png

2.Continuum normalization;#

Since NIRISS spectra may be affected background subtraction & contamination, which could lead to mismatch between filters, here we aim to normaliza spectra from filter to its broadband magnitude.

# Open broadband flux catalog from Notebook 01a;

# Flux are already in Fnu, with magzp = 25.0;
file = DIR_OUT + 'l3_nis_flux.cat'

fd_cat = ascii.read(file)
id_cat = fd_cat['id']
magzp = 25.0

fd_cat
Table length=7
idF309E309F310E310F311E311F308E308F1E1F4E4F6E6F202E202F203E203F204E204F205E205
int64float64float64float64float64float64float64float64float64float64float64float64float64float64float64float64float64float64float64float64float64float64float64
14.176380.20881913.93160.69657815.84890.7924472.398830.1199420.291340.0145671.457470.07287362.224330.1112173.496230.1748126.115050.30575210.56820.52840914.92790.746397
24.176380.20881913.93160.69657815.84890.7924472.398830.1199420.291340.0145671.457470.07287362.224330.1112173.496230.1748126.115050.30575210.56820.52840914.92790.746397
34.176380.20881913.93160.69657815.84890.7924472.398830.1199420.291340.0145671.457470.07287362.224330.1112173.496230.1748126.115050.30575210.56820.52840914.92790.746397
411.29280.56463816.20320.81015915.84890.7924477.684220.3842115.016490.2508256.486340.3243176.86120.343068.863390.4431714.02170.70108416.2480.812416.21810.810905
57.231030.3615516.292160.3146086.309570.3154795.60790.2803950.1752270.008761330.2066330.01033172.711440.1355727.085980.3542996.867520.3433766.468450.3234226.367960.318398
67.231030.3615516.292160.3146086.309570.3154795.60790.2803950.1752270.008761330.2066330.01033172.711440.1355727.085980.3542996.867520.3433766.468450.3234226.367960.318398
77.231030.3615516.292160.3146086.309570.3154795.60790.2803950.1752270.008761330.2066330.01033172.711440.1355727.085980.3542996.867520.3433766.468450.3234226.367960.318398
# Retrieve filter curves for NIRISS images;
url = 'https://jwst-docs.stsci.edu/files/97978948/97978953/1/1596073225227/NIRISS_Filters_2017May.tar.gz'

filename = 'tmp.tar.gz'
urlretrieve(url, filename)
my_tar = tarfile.open(filename)
my_tar.extractall('./')
# Load fiter response:
DIR_FIL = './NIRISS_Filters/DATA/'

# The number corresponds to F200W in EAZY filter response;
eazy_filt = 311

# Read transmission data;
filt_data = ascii.read(f'{DIR_FIL}/NIRISS_{filt.upper()}.txt')

print(filt_data)
wave_filt = filt_data['Wavelength']
flux_filt = filt_data['FilterTrans']

plt.plot(wave_filt, flux_filt, ls='-', label=f'NIRISS {filt.upper()}')
plt.xlabel('wavelength')
plt.ylabel('response')
plt.legend(loc=1)
Wavelength FilterTrans     PCE    
---------- ----------- -----------
       0.5         0.0         0.0
     0.501         0.0         0.0
     0.502         0.0         0.0
     0.503         0.0         0.0
     0.504         0.0         0.0
     0.505         0.0         0.0
     0.506         0.0         0.0
     0.507         0.0         0.0
     0.508         0.0         0.0
     0.509         0.0         0.0
       ...         ...         ...
      5.49 1.14149e-06 5.89206e-08
     5.491 1.15901e-06 5.86208e-08
     5.492 1.17562e-06 5.83295e-08
     5.493 1.19224e-06 5.80189e-08
     5.494 1.20426e-06 5.75614e-08
     5.495 1.21578e-06 5.70597e-08
     5.496  1.2271e-06 5.66351e-08
     5.497 1.23762e-06 5.61508e-08
     5.498 1.24814e-06 5.56559e-08
     5.499 1.25629e-06 5.51494e-08
       5.5 1.26205e-06 5.46377e-08
Length = 5001 rows
<matplotlib.legend.Legend at 0x7fb3bf64ead0>
../../_images/507eb8ce3bb682144ef23605047803a813c52442206e58ea760b9e940c37f35d.png
# Define a small function for flux convolution with filters;
def filconv(lfil, ffil, l0, f0, DIR='', c=3e18):
    '''
    lfil : Wave array for filter response curve.
    ffil : Flux array for filter response curve.

    l0 : Wave array for spectrum, in f_nu, not f_lam
    f0 : Flux array for spectrum, in f_nu, not f_lam

    '''

    fhalf = np.max(ffil)/2.0
    con = (ffil > fhalf)
    lfwhml = np.min(lfil[con])
    lfwhmr = np.max(lfil[con])
    lcen = (lfwhmr + lfwhml)/2.

    lamS, spec = l0, f0  # Two columns with wavelength and flux density
    lamF, filt = lfil, ffil  # Two columns with wavelength and response in the range [0,1]
    filt_int = np.interp(lamS, lamF, filt)  # Interpolate Filter to common(spectra) wavelength axis

    if len(lamS) > 0:
        I1 = simps(spec / lamS**2 * c * filt_int * lamS, lamS)  # Denominator for Fnu
        I2 = simps(filt_int/lamS, lamS) # Numerator
        fnu = I1/I2/c # Average flux density
    else:
        I1 = 0
        I2 = 0
        fnu = 0

    return lcen, fnu
# Convolve observed spectrum with the filter response.
# The input flux need to be in Fnu with magzp=25;

# Cut edge of observed flux;
con = (wave > 1.75) & (wave < 2.25)
lcen, fnu = filconv(wave_filt, flux_filt, wave[con], flux_combine[con], DIR='./')
print('Central wavelength and total flux are ;', lcen, fnu)
Central wavelength and total flux are ; 1.9885 3130.4558232154595
# Get normalization factor;
iix = np.where(id_cat[:] == int(id))
Cnorm = fd_cat[f'F{eazy_filt}'][iix] / fnu 
Cnorm
<Column name='F311' dtype='float64' length=1>
0.005062809026872241
# Wirte the normalized spectrum:
# Make it into a Spectrum1D instance.
file_1d = f'{DIR_OUT}l3_nis_{filt}_{grism}_s{id}_ndither{dither}_1d_opt.fits'

obs = Spectrum1D(spectral_axis=wave*u.um,
                 flux=flux_combine * Cnorm * u.dimensionless_unscaled,
                 uncertainty=StdDevUncertainty(flux_err_combine * Cnorm), unit='')
obs.write(file_1d, format='tabular-fits', overwrite=True)

Repeat 1&2 for other filters too;#

DIR_FIL = './NIRISS_Filters/DATA/'
filts = ['f115w', 'f150w', 'f200w']
# The number corresponds to NIRISS filters in EAZY filter response;
eazy_filts = [309, 310, 311]

# Masks for problematic flux at edge;
mask_lw = [1.05, 1.35, 1.75]
mask_uw = [1.25, 1.65, 2.23]

for ff, filt in enumerate(filts):
    for dither in dithers:
        file_1d = f'{DIR_OUT}l3_nis_{filt}_{grism}_s{id}_ndither{dither}_1d_opt.fits'
        fd = fits.open(file_1d)[1].data
        if dither == 0:
            wave = fd['wavelength']
            flux = np.zeros((ndither, len(wave)), 'float')
            flux_err = np.zeros((ndither, len(wave)), 'float')

            flux[dither, :] = fd['flux']
            flux_err[dither, :] = fd['uncertainty']
        else:
            wave_tmp = fd['wavelength']
            flux_tmp = fd['flux']
            flux_err_tmp = fd['uncertainty']

            flux[dither, :] = np.interp(wave, wave_tmp, flux_tmp)
            flux_err[dither, :] = np.sqrt(np.interp(wave, wave_tmp, flux_err_tmp[:]**2))

    # Combine;
    flux_combine = np.zeros(len(wave), 'float')
    flux_err_combine = np.zeros(len(wave), 'float')

    for ii in range(len(wave)):
        flux_combine[ii] = np.sum(flux[:, ii] * 1 / flux_err[:, ii]**2) / np.sum(1 / flux_err[:, ii]**2) 
        flux_err_combine[ii] = np.sqrt(np.sum(flux_err[:, ii]**2)) / len(flux_err[:, ii])

    # Normalize;
    filt_data = ascii.read(f'{DIR_FIL}/NIRISS_{filt.upper()}.txt')
    wave_filt = filt_data['Wavelength']
    flux_filt = filt_data['FilterTrans']

    con = (wave > mask_lw[ff]) & (wave < mask_uw[ff])
    lcen, fnu = filconv(wave_filt, flux_filt, wave[con], flux_combine[con])
    iix = np.where(id_cat[:] == int(id))
    Cnorm = fd_cat[f'F{eazy_filts[ff]}'][iix] / fnu

    # Write:
    file_1d = f'{DIR_OUT}l3_nis_{filt}_{grism}_s{id}_combine_1d_opt.fits' 
    
    obs = Spectrum1D(spectral_axis=wave*u.um,
                     flux=flux_combine * Cnorm * u.dimensionless_unscaled,
                     uncertainty=StdDevUncertainty(flux_err_combine * Cnorm), unit='')
    obs.write(file_1d, format='tabular-fits', overwrite=True)

Results;#

The spectra are align with each other, as they are matched to BB photometry.

filts = ['f115w', 'f150w', 'f200w']
cols = ['lightblue', 'orange', 'purple']

for ff, filt in enumerate(filts):
    file_1d = f'{DIR_OUT}l3_nis_{filt}_{grism}_s{id}_combine_1d_opt.fits' 
    fd = fits.open(file_1d)[1].data
    wave = fd['wavelength']
    flux = fd['flux']
    flux_err = fd['uncertainty']

    con = (wave > mask_lw[ff]) & (wave < mask_uw[ff])
    plt.plot(wave[con], flux[con], ls='-', label=f'{grism} {filts[ff]}', color=cols[ff])

    filt_data = ascii.read(f'{DIR_FIL}/NIRISS_{filt.upper()}.txt')
    wave_filt = filt_data['Wavelength']
    flux_filt = filt_data['FilterTrans']

    con = (wave > mask_lw[ff]) & (wave < mask_uw[ff])
    lcen, fnu = filconv(wave_filt, flux_filt, wave[con], flux[con])
    iix = np.where(id_cat[:] == int(id))
    plt.scatter(lcen, fd_cat[f'F{eazy_filts[ff]}'][iix], marker='s', s=50, edgecolor='k', color=cols[ff], label=f'BB {filts[ff]}')


plt.xlabel('wavelength')
plt.ylabel('flux')
plt.ylim(0, 60)
plt.legend(loc=0)
<matplotlib.legend.Legend at 0x7fb3bd53ec50>
../../_images/a788ea71569c2a3ec99805b404fcd963f24d563bea5c4204af22bc1575d17df7.png

Repeat for another target;#

id = '00003'

filts = ['f115w', 'f150w', 'f200w']
eazy_filts = [309, 310, 311]
DIR_FIL = './NIRISS_Filters/DATA/'

# Masks for problematic flux at edge;
mask_lw = [1.05, 1.35, 1.75]
mask_uw = [1.25, 1.65, 2.23]

for ff, filt in enumerate(filts):
        
    for dither in dithers:
        file_1d = f'{DIR_OUT}l3_nis_{filt}_{grism}_s{id}_ndither{dither}_1d_opt.fits'
        fd = fits.open(file_1d)[1].data
        if dither == 0:
            wave = fd['wavelength']
            flux = np.zeros((ndither, len(wave)), 'float')
            flux_err = np.zeros((ndither, len(wave)), 'float')

            flux[dither, :] = fd['flux']
            flux_err[dither, :] = fd['uncertainty']
        else:
            wave_tmp = fd['wavelength']
            flux_tmp = fd['flux']
            flux_err_tmp = fd['uncertainty']

            flux[dither, :] = np.interp(wave, wave_tmp, flux_tmp)
            flux_err[dither, :] = np.sqrt(np.interp(wave, wave_tmp, flux_err_tmp[:]**2))

    # Combine;
    flux_combine = np.zeros(len(wave), 'float')
    flux_err_combine = np.zeros(len(wave), 'float')

    for ii in range(len(wave)):
        flux_combine[ii] = np.sum(flux[:, ii] * 1 / flux_err[:, ii]**2) / np.sum(1 / flux_err[:, ii]**2) 
        flux_err_combine[ii] = np.sqrt(np.sum(flux_err[:, ii]**2)) / len(flux_err[:, ii])

    # Normalize;
    filt_data = ascii.read(f'{DIR_FIL}/NIRISS_{filt.upper()}.txt')
    wave_filt = filt_data['Wavelength']
    flux_filt = filt_data['FilterTrans']

    con = (wave > mask_lw[ff]) & (wave < mask_uw[ff])
    lcen, fnu = filconv(wave_filt, flux_filt, wave[con], flux_combine[con])
    iix = np.where(id_cat[:] == int(id))
    Cnorm = fd_cat[f'F{eazy_filts[ff]}'][iix] / fnu
        
    # Write:
    file_1d = f'{DIR_OUT}l3_nis_{filt}_{grism}_s{id}_combine_1d_opt.fits'
    obs = Spectrum1D(spectral_axis=wave*u.um,
                     flux=flux_combine * Cnorm * u.dimensionless_unscaled,
                     uncertainty=StdDevUncertainty(flux_err_combine * Cnorm), unit='')
    obs.write(file_1d, format='tabular-fits', overwrite=True)
filts = ['f115w', 'f150w', 'f200w']
cols = ['lightblue', 'orange', 'purple']

for ff, filt in enumerate(filts):
    file_1d = f'{DIR_OUT}l3_nis_{filt}_{grism}_s{id}_combine_1d_opt.fits' 
    fd = fits.open(file_1d)[1].data
    wave = fd['wavelength']
    flux = fd['flux']
    flux_err = fd['uncertainty']

    con = (wave > mask_lw[ff]) & (wave < mask_uw[ff])
    plt.plot(wave[con], flux[con], ls='-', label=f'{grism} {filts[ff]}', color=cols[ff])

    filt_data = ascii.read(f'{DIR_FIL}/NIRISS_{filt.upper()}.txt')
    wave_filt = filt_data['Wavelength']
    flux_filt = filt_data['FilterTrans']

    con = (wave > mask_lw[ff]) & (wave < mask_uw[ff])
    lcen, fnu = filconv(wave_filt, flux_filt, wave[con], flux[con])
    iix = np.where(id_cat[:] == int(id))
    plt.scatter(lcen, fd_cat[f'F{eazy_filts[ff]}'][iix], marker='s', s=50, edgecolor='k', color=cols[ff], label=f'BB {filts[ff]}')

plt.xlabel('wavelength')
plt.ylabel('flux')
plt.ylim(0, 60)
plt.legend(loc=0)
<matplotlib.legend.Legend at 0x7fb3bd2742d0>
../../_images/74e7b8fe3319ccd61403cd6d6b0bb7a0b2e6a4ec8e1f7726636c77aeaf63566c.png