WFSS Spectra Part 0: Optimal Extraction#

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 1 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 start with post-pipeline products of NIRISS WFSS, 2D rectified spectra, from spec level3.

Optimal extraction requires source morphology along the cross-dispersion direction, which will be retrieved from direct images taken along with WFSS observations. Morphology along dispersion direction is also essential to infer the spectral resolution, which will be obtained using template fitting to get redshift and stellar population in notebook #3 of this set.

Note: We here assume reduction of the 2D rectified spectrum has been performed at a decent level, i.e. there is no contaminating flux from other sources on the target 2D spectrum, and that background is already subtracted.

%matplotlib inline
import os
import numpy as np
from scipy.ndimage import rotate
from scipy.optimize import curve_fit

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

from specutils import Spectrum1D
from astropy.nddata import StdDevUncertainty

import specutils
print('specutils', specutils.__version__)

import astropy
print('astropy', astropy.__version__)
specutils 1.12.0
astropy 6.0.0

The version should be#

  • specutils 1.11.0

  • astropy 5.3.4

import matplotlib.pyplot as plt
import matplotlib as mpl

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

0.Download and load data:#

These include pipeline processed data for NIRISS, as well as photometric catalog from image step3.

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()
DIR_DATA = './pipeline_products/'

# Output directory;
DIR_OUT = './output/'
if not os.path.exists(DIR_OUT):
    os.mkdir(DIR_OUT)

# Filter for detection and science image;
filt_det = 'f200w'

# Image array from direct image. This is for optimal extraction and masking.
# This image should already be sky-subtracted; otherwise, you will encounter a wrong result with optimal extraction.

infile = f'{DIR_DATA}l3_nis_{filt_det}_i2d_skysub.fits'
hdu = fits.open(infile)

# This is just for error array;
infile = f'{DIR_DATA}l3_nis_{filt_det}_i2d.fits'
hdu_err = fits.open(infile)

data = hdu[0].data
imwcs = wcs.WCS(hdu[0].header, hdu)

err = hdu_err[2].data
weight = 1/np.square(err)

# Segmentation map;
# This can be prepared by running Photutils, if the pipeline does not generate one.
segfile = f'{DIR_DATA}l3_nis_{filt_det}_i2d_seg.fits'
seghdu = fits.open(segfile)
segdata = seghdu[0].data
/tmp/ipykernel_2017/874489203.py:25: RuntimeWarning: divide by zero encountered in divide
  weight = 1/np.square(err)
# Load catalog from image level3;
# to obtain source position in pixel coordinate.
catfile = f'{DIR_DATA}l3_nis_{filt_det}_cat.ecsv'
fd = ascii.read(catfile)
fd
QTable length=7
idxcentroidycentroidsky_centroidaper_bkg_fluxaper_bkg_flux_erraper30_fluxaper30_flux_erraper50_fluxaper50_flux_erraper70_fluxaper70_flux_erraper_total_fluxaper_total_flux_erraper30_abmagaper30_abmag_erraper50_abmagaper50_abmag_erraper70_abmagaper70_abmag_erraper_total_abmagaper_total_abmag_erraper30_vegamagaper30_vegamag_erraper50_vegamagaper50_vegamag_erraper70_vegamagaper70_vegamag_erraper_total_vegamagaper_total_vegamag_errCI_30_50CI_50_70CI_30_70is_starsharpnessroundnessnn_distnn_abmagisophotal_fluxisophotal_flux_errisophotal_abmagisophotal_abmag_errisophotal_vegamagisophotal_vegamag_errisophotal_areasemimajor_sigmasemiminor_sigmaellipticityorientationsky_orientationsky_bbox_llsky_bbox_ulsky_bbox_lrsky_bbox_ur
pixpixdeg,degJyJyJyJyJyJyJyJyJyJypixJyJypix2pixpixdegdegdeg,degdeg,degdeg,degdeg,deg
int64float64float64SkyCoordfloat64float64float64float64float64float64float64float64float64float64float64float64float64float64float64float64float64float64float64float64float64float64float64float64float64float64float64float64float64float64float64float32float64float64float64float32float64float64float64float64float64float64float64float64float64float64SkyCoordSkyCoordSkyCoordSkyCoord
11901.4791616.66173.569653944921018,-30.4050701867237763.483045e-08nan8.130363e-08nan1.921404e-07nan3.765590e-07nannannan26.624725nan25.690953nan24.960418nannannan24.915445nan23.981673nan23.251138nannannan0.93380.73051.6643nan0.6430530.050833687.14285622.2450993.217875e-06nan22.631077nan20.921797nan83.02.3000092.2493400.022030-68.27240522.3448333.569758007135362,-30.4051652899355463.5697605503538763,-30.4049647348617683.56952547091111,-30.4051630962729043.5695280146070885,-30.404962541203844
21215.5473657.44113.5841636319830568,-30.404462674434824.591611e-08nan1.550019e-07nan3.638838e-07nan6.725310e-07nannannan25.924158nan24.997593nan24.330719nannannan24.214878nan23.288313nan22.621439nannannan0.92660.66691.5934nan0.6664810.076069687.14285622.6310774.591563e-06nan22.245099nan20.535819nan95.02.9689641.7782620.40105089.747234180.3644713.5842476045730627,-30.4045900237988253.5842510319996697,-30.4043165392158523.5840573485032308,-30.4045882500167653.58406077646256,-30.404314765438812
3381.9250984.51653.601859711428753,-30.3986624377804484.268341e-08nan1.404209e-07nan3.172360e-07nan6.277053e-07nannannan26.031420nan25.146544nan24.405610nannannan24.322140nan23.437264nan22.696330nannannan0.88490.74091.6258nan0.6721850.331872673.40698616.3206584.647145e-06nan22.232034nan20.522754nan101.02.8800191.9896740.309146-89.2374631.3797753.6019516675347565,-30.398791223651383.601955045183127,-30.398517738611013.6017402841066257,-30.3987892812738743.60174366234683,-30.398515796238758
41100.56771408.44083.5867654957418873,-30.3907928519617564.182255e-08nan1.977441e-07nan4.634087e-07nan8.557904e-07nannannan25.659741nan24.735090nan24.069081nannannan23.950461nan23.025810nan22.359801nannannan0.92470.66601.5907nan0.6758910.239616759.75060622.2450994.479109e-06nan22.272021nan20.562741nan78.02.8643311.5241970.46787089.968483180.5857203.586849890460872,-30.3909201996267733.586853310121526,-30.3906467149754253.586659660970119,-30.3909184296174933.586663081163341,-30.39064494497111
5662.42881685.17023.596088874688556,-30.3858332726452253.108490e-08nan6.850841e-07nan1.445057e-06nan2.486505e-06nannannan24.310640nan23.500287nan22.911026nannannan22.601360nan21.791007nan21.201746nannannan0.81040.58931.3996nan0.584634-0.067064334.07677220.7272655.441945e-06nan22.060615nan20.351335nan71.01.8741131.8202290.02875210.431093101.0483303.5961919899951726,-30.385919383231273.5961942522168715,-30.3857370599751853.5959594989853887,-30.385917236390353.5959617616410102,-30.38573491313818
6319.49901885.30353.60338188710553,-30.3822511019631463.326749e-08nan6.893335e-07nan1.514847e-06nan2.523116e-06nannannan24.303927nan23.449078nan22.895157nannannan22.594647nan21.739798nan21.185877nannannan0.85480.55391.4088nan0.5817810.158721111.50410922.0283285.666923e-06nan22.016632nan20.307352nan75.01.8944241.8663220.01483414.267744104.8849813.603486459620283,-30.382339651593853.603488932995949,-30.3821390958841353.6032539769794103,-30.382337517665473.6032564508323546,-30.382136961959965
7266.91821983.63173.604515260247416,-30.3804685451213073.310798e-08nan7.155692e-07nan1.476385e-06nan2.520043e-06nannannan24.263371nan23.477001nan22.896480nannannan22.554091nan21.767721nan21.187200nannannan0.78640.58051.3669nan0.5937630.061278111.50410922.0166325.606203e-06nan22.028328nan20.319048nan74.01.8938071.8479210.02423089.641888180.2591263.6046074836711655,-30.3805629644681233.6046099547002646,-30.3803624087416543.6043961396388418,-30.380561026353113.6043986111017956,-30.38036047063047

Make a broadband flux catalog.#

  • For a convenient reason, here we compile catalogs into a flux catalog, which will be used in the following notebook (01b).

  • To run this cell, you will need a photometric catalog of sources, that list sources position and flux for each filter. For now, I use this catalog prepared in another notebook. (“sources_extend_01.cat”)

  • This catalog can also be used for generic phot-z/SED fitting softwares, like EAZY and gsf (see notebook No.04).

For now, we use an input catalog.#

filts = ['f115w', 'f150w', 'f200w', 'f090w', 'f435w', 'f606w', 'f814w', 'f105w', 'f125w', 'f140w', 'f160w']
eazy_filts = [309, 310, 311, 308, 1, 4, 6, 202, 203, 204, 205]
magzp = 25.0 # magnitude zeropoint in the catalog.

# Read catalog;
fd_input = ascii.read(f'{DIR_DATA}sources_extend_01.cat')
ra_input = fd_input['x_or_RA']
dec_input = fd_input['y_or_Dec']

# Header;
fw = open(f'{DIR_OUT}l3_nis_flux.cat', 'w')
fw.write('# id')
for ff in range(len(filts)):
    fw.write(f' F{eazy_filts[ff]} E{eazy_filts[ff]}')
fw.write('\n')

# Contents;
for ii in range(len(fd['id'])):
    
    rtmp = np.sqrt((fd['sky_centroid'].ra.value[ii] - ra_input[:])**2 + (fd['sky_centroid'].dec.value[ii] - dec_input[:])**2)
    iix = np.argmin(rtmp)
    
    for ff in range(len(filts)):
        if ff == 0:
            fw.write(f"{fd['id'][ii]}")

        mag = fd_input[f'niriss_{filts[ff]}_magnitude'][iix]
        flux_nu = 10**((mag-magzp)/(-2.5))

        # Currently, the catalog does not provide proper error;
        # Assuming a 5% error for flux.
        
        flux_err_nu = flux_nu * 0.05

        fw.write(f' {flux_nu:.5e} {flux_err_nu:.5e}') 

    fw.write('\n')
fw.close()

1.Load 2D spectrum;#

# Which filter, grating, and object?
filt = 'f200w'

# grism = 'G150R'
grism = 'G150C'

id = '00004'

# Zero-indexed number for dither --- the test data here has two dither positions, so 0 or 1.
ndither = 0

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)

# Get position angle of observation;
hd_2d = hdu_2d[1].header
PA_V3 = hd_2d['PA_V3']
PA_V3
0.0
plt.imshow(data_2d, vmin=0, vmax=1000)
<matplotlib.image.AxesImage at 0x7fa114afba90>
../../_images/cb5a831ca2c6ac1c09fd75756b8f652fd03500ef6d9c69b9918648f322fba6a5.png
# Get light profile of the source;

# Again, y is for cross-dispersion, and x is for dispersion directions.
y2d, x2d = data_2d.shape[:]

# Cut out segmentation map;
iix = np.where(fd['id'] == int(id))[0][0]

# Target position from image 3 catalog;
ycen = fd['ycentroid'][iix].value
xcen = fd['xcentroid'][iix].value

# Cutout size = y direction of 2D spectrum;
rsq = y2d

sci_cut = data[int(ycen-rsq/2.+0.5):int(ycen+rsq/2.+0.5), int(xcen-rsq/2.+0.5):int(xcen+rsq/2.+0.5)]
seg_cut = segdata[int(ycen-rsq/2.+0.5):int(ycen+rsq/2.+0.5), int(xcen-rsq/2.+0.5):int(xcen+rsq/2.+0.5)]

# Rotate image for PA of Grism observation;
if grism == 'G150C':
    sci_rot = rotate(sci_cut, PA_V3)
else:
    sci_rot = rotate(sci_cut, PA_V3+90)

WFSS grism is dispersed in a direction of x-axis in the plot below.#

plt.imshow(sci_rot, vmin=0, vmax=1.0)
plt.title('Direct image')
plt.xlabel('Wavelength direction >>>', color='r', fontsize=18)
plt.ylabel('Cross-dispersion direction >>>', color='r', fontsize=18)
Text(0, 0.5, 'Cross-dispersion direction >>>')
../../_images/82ed68640d2491ef8924b0e63d12c874e2b3f27e33cf1a1065c998fbf1301fbe.png

2.Get light profile at different x position — This will be used for optimal extraction.#

for ii in range(sci_rot.shape[1]):
    flux_tmp = sci_rot[:, ii]
    xx_tmp = np.arange(0, len(sci_rot[:, ii]), 1)
    plt.plot(xx_tmp, flux_tmp, label='x={}'.format(ii))
plt.legend(loc=0, fontsize=8)
<matplotlib.legend.Legend at 0x7fa114a4c3d0>
../../_images/74c89de1f3b441d0aa3f643983b87b4c0f124f12d996b29044cdfdf69fd607dc.png
# Sum along x (disperse) direction
flux_y = np.zeros(len(sci_rot[:, 0]), 'float')
for ii in range(sci_rot.shape[0]):
    flux_y[ii] = np.sum(sci_rot[ii, :])

# Sky subtraction, if needed.
# sky = np.mean([flux_y[0], flux_y[-1]])

# Normalize;
flux_y[:] /= flux_y.sum()

plt.plot(xx_tmp, flux_y)
plt.xlabel('y-position')
plt.ylabel('Source Flux')
Text(0, 0.5, 'Source Flux')
../../_images/1101f87f64286f25bed25a5e380301bfdce4483684ef2f5b49becffc949e1429.png

3.One-dimensional extraction;#

Show pipeline 1D extraction as an example;#

# Normal extraction;
flux_disp1 = np.zeros(x2d, 'float')
err_disp1 = np.zeros(x2d, 'float')
wave_disp1 = np.zeros(x2d, 'float')
    
for ii in range(x2d): # Wavelength direction.
    mask_tmp = (dq_2d[:, ii] == 0) & (err_2d[:, ii] > 0)

    # Sum within a box;
    flux_disp1[ii] = np.sum(data_2d[:, ii][mask_tmp]) 
    err_disp1[ii] = np.sqrt(np.sum(err_2d[:, ii][mask_tmp]**2)) 
    wave_disp1[ii] = wave_2d[0, ii]

plt.errorbar(wave_disp1, flux_disp1, yerr=err_disp1)
plt.xlim(1.7, 2.3)
(1.7, 2.3)
../../_images/eef9e5ca249f8af795677bb564c1d157ff513f9d5937380da00dbdf3460ad452.png

Optimal extraction;#

# Following Horne(1986, PASP, 98, 609);
flux_disp = np.zeros(x2d, 'float')
err_disp = np.zeros(x2d, 'float')
wave_disp = np.zeros(x2d, 'float')

# Sigma clipping.
sig = 5.0

for ii in range(x2d): # ii : wavelength element.
    # Mask; 
    # 1. DQ array
    # 2. error value
    # 3. CR detection
    mask_tmp = (dq_2d[:, ii] == 0) & (err_2d[:, ii] > 0) & ((data_2d[:, ii] - flux_y[:] * flux_disp1[ii])**2 < sig**2 * err_2d[:, ii]**2)
    ivar = 1. / err_2d[:, ii]**2

    num = flux_y[:] * data_2d[:, ii] * ivar
    den = flux_y[:]**2 * ivar
    flux_disp[ii] = num[mask_tmp].sum(axis=0) / den[mask_tmp].sum(axis=0)
    err_disp[ii] = np.sqrt(1./den[mask_tmp].sum(axis=0))
    wave_disp[ii] = wave_2d[0, ii]
    
plt.errorbar(wave_disp, flux_disp, yerr=err_disp)
plt.xlim(1.7, 2.3)
(1.7, 2.3)
../../_images/6c1909b190e100c5dcff86abe4a9814e44b2161dae91536eff6818fddaa8596c.png
# Compare;
plt.errorbar(wave_disp, flux_disp, yerr=err_disp, color='r', label='Optimal')
plt.errorbar(wave_disp1, flux_disp1, yerr=err_disp1, color='b', alpha=0.5, label='Box')
plt.ylim(-10, 20000)
plt.legend(loc=0)
plt.xlabel('Wavelength')
Text(0.5, 0, 'Wavelength')
../../_images/d9e33745d8b5413e945dd7ea1d118a1b1b73754b279e8d3472782226ed4cf73e.png

4.Write 1d spectrum out to a file;#

file_1d = f'{DIR_OUT}l3_nis_{filt}_{grism}_s{id}_1d_opt.fits'

# Now make it into a Spectrum1D instance.
obs = Spectrum1D(spectral_axis=wave_disp*u.um,
                 flux=flux_disp*u.MJy,
                 uncertainty=StdDevUncertainty(err_disp), unit='MJy')
obs.write(file_1d, format='tabular-fits', overwrite=True)

5.Light profile along x-axis = Resolution of dispersed spectrum;#

As WFSS does not have a slit, any dispersed spectrum is affected by source morphology. The estimate on the effective spectral resolution will be needed in the following notebook. And we here try to estimate it beforehand;

for ii in range(sci_rot.shape[0]):
    flux_tmp = sci_rot[ii, :]
    xx_tmp = np.arange(0, len(sci_rot[ii, :]), 1)
    plt.plot(xx_tmp, flux_tmp, label=f'y={ii}')
    
plt.legend(loc=1, fontsize=8)
plt.xlabel('Wavelength direction')
plt.title('Source light profile along dispersed direction', fontsize=14)
Text(0.5, 1.0, 'Source light profile along dispersed direction')
../../_images/6b2750ff01e2d5e5d50f7c3b6a55e19a51740c92885b1d447c836e1930602061.png

*Unless you are interested in spatially resolved spectra, you can stack and get light profile as a good approximation.#

# Sum along cross-disperse direction
flux_x = np.zeros(len(sci_rot[0, :]), 'float')
for ii in range(sci_rot.shape[0]):
    flux_x[ii] = np.sum(sci_rot[:, ii])

# Normalize;
flux_x[:] /= flux_x.sum()

plt.plot(xx_tmp, flux_x, label='Convolution kernel')
plt.legend(loc=2)
<matplotlib.legend.Legend at 0x7fa10f2b3910>
../../_images/3119a76fdb7a641a32407e180f96e9a9ed995f1aae391ce5bce6af785bcf6711.png

Fit with a moffat function;#

# Fitting function with Moffat
# Moffat fnc.

def moffat(xx, A, x0, gamma, alp):
    yy = A * (1. + (xx-x0)**2/gamma**2)**(-alp)
    return yy


def fit_mof(xx, lsf):
    # xx = lsf * 0
    # for ii in range(len(lsf)):
    #    xx[ii] = ii - len(lsf)/2.
    popt, pcov = curve_fit(moffat, xx, lsf)
    return popt


def LSF_mof(xsf, lsf, f_plot=True):
    '''
    Input:
    =======
    xsf : x axis for the profile.
    lsf : light profile.    
    '''
    
    # for ii in range(len(sci[0, :])):
    #    lsf[ii] = np.mean(sci_rot[int(height/2.-5):int(height/2.+5), ii])
    #    xsf[ii] = ii - len(lsf)/2.

    try:
        A, xm, gamma, alpha = fit_mof(xsf, lsf)
    except RuntimeError:
        print('Fitting failed.')
        A, xm, gamma, alpha = -1, -1, -1, -1
        pass

    if A > 0:
        lsf_mod = moffat(xsf, A, 0, gamma, alpha)
        
    if f_plot:
        yy = moffat(xsf, A, xm, gamma, alpha)
        plt.plot(xsf, yy, 'r.', ls='-', label='Data')
        plt.plot(xsf, lsf_mod, 'b+', ls='-', label=f'Model: gamma={gamma:2f}\nalpha={alpha:2f}')
        plt.legend()
        plt.show()
    
    return A, xm, gamma, alpha
# LSF, line spread function
iix_peak = np.argmax(flux_x)
xx_tmp_shift = xx_tmp - xx_tmp[iix_peak]
A, xm, gamma, alpha = LSF_mof(xx_tmp_shift, flux_x)
../../_images/d3c444b70dfc1c56bf9c99db324bf4a6416b58f970f022da5b45a9b321378fa7.png
# Write it down;
# Tha parameters are in unit of pixel.
fm = open(f'{DIR_OUT}l3_nis_{filt}_{grism}_s{id}_moffat.txt', 'w')
fm.write('# A x0 gamma alp\n')
fm.write('# Moffat function\n')
fm.write(f'{A:.3f} {xm:.3f} {gamma:.3f} {alpha:.3f}\n')

fm.close()

Repeat for other filters, other objects.#

The following big colum executes the same processes above for other filters and dither position.#

grism = 'G150C'
id = '00004'
DIR_OUT = './output/'
if not os.path.exists(DIR_OUT):
    os.mkdir(DIR_OUT)

filts = ['f115w', 'f150w', 'f200w']
ndithers = np.arange(0, 2, 1)

sig = 5.0

for filt in filts:
    print(filt)
    
    # 2d spectrum;
    file_2d = f'{DIR_DATA}l3_nis_{filt}_{grism}_s{id}_cal.fits'
    hdu_2d = fits.open(file_2d)

    for ndither in ndithers:
        print(ndither)

        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)

        y2d, x2d = data_2d.shape[:]

        plt.close()
        plt.imshow(data_2d, vmin=0, vmax=300)
        plt.show()

        # Re-extract 2d image;
        # if ndither == 0:
        rsq = y2d
        sci_cut = data[int(ycen-rsq/2.+0.5):int(ycen+rsq/2.+0.5), int(xcen-rsq/2.+0.5):int(xcen+rsq/2.+0.5)]
        seg_cut = segdata[int(ycen-rsq/2.+0.5):int(ycen+rsq/2.+0.5), int(xcen-rsq/2.+0.5):int(xcen+rsq/2.+0.5)]

        # Not sure if the offset in extractioin box is bug ;
        if grism == 'G150C':
            sci_rot = rotate(sci_cut, PA_V3+0)
        else:
            sci_rot = rotate(sci_cut, PA_V3+0+90)

        # This is for spectral resolution;
        # Get light profile along the x-axis
        for ii in range(sci_rot.shape[0]):
            flux_tmp = sci_rot[ii, :]
            xx_tmp = np.arange(0, len(sci_rot[ii, :]), 1)

        # Sum along cross-disperse direction
        flux_x = np.zeros(len(sci_rot[0, :]), 'float')
        for ii in range(sci_rot.shape[0]):
            flux_x[ii] = np.sum(sci_rot[ii, :])

        # Normalize;
        flux_x[:] /= flux_x.sum()

        # LSF
        iix_peak = np.argmax(flux_x)
        xx_tmp_shift = xx_tmp - xx_tmp[iix_peak]
        A, xm, gamma, alpha = LSF_mof(xx_tmp_shift, flux_x)

        if ndither == 0:
            # Write it down;
            fm = open(f'{DIR_OUT}l3_nis_{filt}_{grism}_s{id}_moffat.txt', 'w')
            fm.write('# A x0 gamma alp\n')
            fm.write('# Moffat function\n')
            fm.write(f'{A:.3f} {xm:.3f} {gamma:.3f} {alpha:.3f}\n')
            fm.close()

        # This is for Optimal extraction;
        # Sum along x (disperse) direction
        flux_y = np.zeros(len(sci_rot[:, 0]), 'float')
        for ii in range(sci_rot.shape[0]):
            flux_y[ii] = np.sum(sci_rot[ii, :])
            
        # Normalize;
        flux_y[:] /= flux_y.sum()

        # Following Horne;
        flux_disp = np.zeros(x2d, 'float')
        err_disp = np.zeros(x2d, 'float')
        wave_disp = np.zeros(x2d, 'float')

        for ii in range(x2d):
            # Mask; 
            # 1. DQ array
            # 2. error value
            # 3. CR detection
            mask_tmp = (dq_2d[:, ii] == 0) & (err_2d[:, ii] > 0)
            ivar = 1. / err_2d[:, ii]**2

            num = flux_y[:] * data_2d[:, ii] * ivar 
            den = flux_y[:]**2 * ivar
            flux_disp[ii] = num[mask_tmp].sum(axis=0)/den[mask_tmp].sum(axis=0)
            err_disp[ii] = np.sqrt(1./den[mask_tmp].sum(axis=0))
            wave_disp[ii] = wave_2d[0, ii]

        plt.close()
        con_plot = (wave_disp > 0)
        plt.errorbar(wave_disp[con_plot], flux_disp[con_plot], yerr=err_disp[con_plot])
        plt.ylim(-0, 3000)
        plt.show()

        # Wirte:
        # Now make it into a Spectrum1D instance.
        file_1d = f'{DIR_OUT}l3_nis_{filt}_{grism}_s{id}_ndither{ndither}_1d_opt.fits'

        if wave_disp[1] - wave_disp[0] < 0:
            obs = Spectrum1D(spectral_axis=wave_disp[::-1]*u.um,
                             flux=flux_disp[::-1]*u.MJy,
                             uncertainty=StdDevUncertainty(err_disp[::-1]), unit='MJy')
        else:
            obs = Spectrum1D(spectral_axis=wave_disp*u.um,
                             flux=flux_disp*u.MJy,
                             uncertainty=StdDevUncertainty(err_disp), unit='MJy')
            
        obs.write(file_1d, format='tabular-fits', overwrite=True)
f115w
0
../../_images/997b557d9a29626e7c614eac55800573d0a04d9bc5fe3e7c26fb83bd0a75e6b0.png ../../_images/a723234642e585224846c1c837653bc4506022317fba9ee0234d65b1950a4422.png ../../_images/7675af5bfbfd9a75b4146f6744e234d8256f041ff05945ac21fdab5ed497cd09.png
1
../../_images/0ddbf85e886117a7faf9d0b55b493c160f8c03990a69efe218d8225520b9c9e0.png ../../_images/a723234642e585224846c1c837653bc4506022317fba9ee0234d65b1950a4422.png ../../_images/b9182a1959dd0f794872383815f58cf86a73e2927d9751f43aa8368536327f38.png
f150w
0
../../_images/46649887ac522caed63289d26914d5c72f1ce79850d85f8cef09caa330f5759a.png ../../_images/a723234642e585224846c1c837653bc4506022317fba9ee0234d65b1950a4422.png ../../_images/14d14ef0cb0608e94595ef10e1db2e15ce18c1cd39fe585e08a0af99f4606387.png
1
../../_images/ce16a1b0aa374159e358b86ad9fcf279fb1083b11bfb549dc194178aa8f8cc19.png ../../_images/a723234642e585224846c1c837653bc4506022317fba9ee0234d65b1950a4422.png ../../_images/3273b579ccc136a4eba796dad619fbf5d589d18ff798ea1b72aa92676eea5c3b.png
f200w
0
../../_images/543d61d5992615e80aa43a5e051fd14ab98fdbb33757853af12823e0649e45b7.png ../../_images/8171a4498fb03e5af348610d41136a75186235150a768f568ac13246c31254b6.png ../../_images/d6b60e86af3ebe93897066a349855817bb94abdd3981d9216780990c485ab2d3.png
1
../../_images/8f46d3b13ba1eb8c255d9c09793522f94acb67f745d4908fbd4834ab77b59fa8.png ../../_images/a723234642e585224846c1c837653bc4506022317fba9ee0234d65b1950a4422.png ../../_images/af95714152a04bab8d7547977f7e6ee42e5625ba56dadc20b8b4a1fded82bce0.png

Another object;#

Absorption line galaxy

grism = 'G150C'
id = '00003'
DIR_OUT = './output/'

filts = ['f115w', 'f150w', 'f200w']
ndithers = np.arange(0, 2, 1) # There are four dithers in the data set;

sig = 5.0

for filt in filts:
    print(filt)
    # 2d spectrum;
    file_2d = f'{DIR_DATA}l3_nis_{filt}_{grism}_s{id}_cal.fits'
    hdu_2d = fits.open(file_2d)

    for ndither in ndithers:
        print(ndither)
        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)

        y2d, x2d = data_2d.shape[:]

        plt.close()
        plt.imshow(data_2d, vmin=0, vmax=20)
        plt.show()

        # Re-extract 2d image;
        # if ndither == 0:
        rsq = y2d
        sci_cut = data[int(ycen-rsq/2.+0.5):int(ycen+rsq/2.+0.5), int(xcen-rsq/2.+0.5):int(xcen+rsq/2.+0.5)]
        seg_cut = segdata[int(ycen-rsq/2.+0.5):int(ycen+rsq/2.+0.5), int(xcen-rsq/2.+0.5):int(xcen+rsq/2.+0.5)]

        # Not sure if the offset in extraction box is bug ;
        if grism == 'G150C':
            sci_rot = rotate(sci_cut, PA_V3+0)
        else:
            sci_rot = rotate(sci_cut, PA_V3+0+90)

        # This is for spectral resolution;
        # Get light profile along the x-axis
        for ii in range(sci_rot.shape[0]):
            flux_tmp = sci_rot[ii, :]
            xx_tmp = np.arange(0, len(sci_rot[ii, :]), 1)

        # Sum along cross-disperse direction
        flux_x = np.zeros(len(sci_rot[0, :]), 'float')
        for ii in range(sci_rot.shape[0]):
            flux_x[ii] = np.sum(sci_rot[ii, :])

        # Normalize;
        flux_x[:] /= flux_x.sum()

        # LSF
        iix_peak = np.argmax(flux_x)
        xx_tmp_shift = xx_tmp - xx_tmp[iix_peak]
        A, xm, gamma, alpha = LSF_mof(xx_tmp_shift, flux_x)

        if ndither == 0:
            # Write it down;
            fm = open(f'{DIR_OUT}l3_nis_{filt}_{grism}_s{id}_moffat.txt', 'w')
            fm.write('# A x0 gamma alp\n')
            fm.write('# Moffat function\n')
            fm.write(f'{A:.3f} {xm:.3f} {gamma:.3f} {alpha:.3f}\n')
            fm.close()

        # This is for Optimal extraction;
        # Sum along x (disperse) direction
        flux_y = np.zeros(len(sci_rot[:, 0]), 'float')
        for ii in range(sci_rot.shape[0]):
            flux_y[ii] = np.sum(sci_rot[ii, :])
    
        # Normalize;
        flux_y[:] /= flux_y.sum()

        # Following Horne;
        flux_disp = np.zeros(x2d, 'float')
        err_disp = np.zeros(x2d, 'float')
        wave_disp = np.zeros(x2d, 'float')

        for ii in range(x2d):
            # Mask; 
            # 1. DQ array
            # 2. error value
            # 3. CR detection
            mask_tmp = (dq_2d[:, ii] == 0) & (err_2d[:, ii] > 0) 
            ivar = 1. / err_2d[:, ii]**2

            num = flux_y[:] * data_2d[:, ii] * ivar 
            den = flux_y[:]**2 * ivar
            flux_disp[ii] = num[mask_tmp].sum(axis=0)/den[mask_tmp].sum(axis=0)
            err_disp[ii] = np.sqrt(1./den[mask_tmp].sum(axis=0))
            wave_disp[ii] = wave_2d[0, ii]

        plt.close()
        con_plot = (wave_disp > 0)
        plt.errorbar(wave_disp[con_plot], flux_disp[con_plot], yerr=err_disp[con_plot])
        plt.ylim(-20, 100)
        plt.show()

        # Write:
        # Now, make it into a Spectrum1D instance.
        file_1d = f'{DIR_OUT}l3_nis_{filt}_{grism}_s{id}_ndither{ndither}_1d_opt.fits'

        if wave_disp[1] - wave_disp[0] < 0:
            obs = Spectrum1D(spectral_axis=wave_disp[::-1]*u.um,
                             flux=flux_disp[::-1]*u.MJy,
                             uncertainty=StdDevUncertainty(err_disp[::-1]), unit='MJy')
        else:
            obs = Spectrum1D(spectral_axis=wave_disp*u.um,
                             flux=flux_disp*u.MJy,
                             uncertainty=StdDevUncertainty(err_disp), unit='MJy')

        obs.write(file_1d, format='tabular-fits', overwrite=True)    
f115w
0
../../_images/f5e475476d82e7f5c05c65490cce234914101cbab2316d02f44889690d67146c.png ../../_images/2ce48fbe087381443d2d5106a44a99dcb62da7992bceab40084fa8d41f5a4d17.png ../../_images/d1904518e0931d38fbc625f6b4e5abdf9548def6f7321526652d434ba28b3c1b.png
1
../../_images/4386ec415bb0c461d97baad3006d6ef0cba4f725a3b76f46037c3c870953a643.png ../../_images/2ce48fbe087381443d2d5106a44a99dcb62da7992bceab40084fa8d41f5a4d17.png ../../_images/02dc990d3ea5a494ff88a8acf0432e63802aeced99492c646c8ada2f59d5e7ac.png
f150w
0
../../_images/c779ec22750e478febd3e5f9209d6626cea5003c0813bf4bcb094f24aafd88cd.png ../../_images/2ce48fbe087381443d2d5106a44a99dcb62da7992bceab40084fa8d41f5a4d17.png ../../_images/fb9b11a655a79843f1c9e068d1cf7f2124763f03372c0007b1abc9ac6e3df5bd.png
1
../../_images/655a48a7d11eb325370dc5f37fd3945114ca75e648d16ef2489b43312bff9220.png ../../_images/f339e5ccc733442c55f02c02f543f838dd35837b4dd0178b1d673fcac1241ed4.png ../../_images/0ad1239dabe110165b0fb3e7bd4b29a31d6270bfffd5655f172e480266fad721.png
f200w
0
../../_images/8901a714efb20ccc7bc3670c711b3aef301a316f036fe6712c6269bb7b1b6a6a.png ../../_images/2ce48fbe087381443d2d5106a44a99dcb62da7992bceab40084fa8d41f5a4d17.png ../../_images/a28e32ffd5ecb5ab5931151030fef3c11bdd33fa785d96f4fc30c4588a58d548.png
1
../../_images/03197a43f0dcceb67ec8f6611f4a5afa4cca38be010f7628f4023791cb96d8ac.png ../../_images/2ce48fbe087381443d2d5106a44a99dcb62da7992bceab40084fa8d41f5a4d17.png ../../_images/958ac7ca03809408b26193a8288942263f9a099df3558a58e1490f0c15658709.png