1. Load NIRISS pipeline output#

We get the data directly from AWS by searching for proposal 2079 and NIRISS instrument in WFSS mode. The filename we are looking for is jw02079-o001_s000000680_niriss_f150w-gr150c_x1d.fits. Uncomment the line with download_products if you want to save a local copy of the file.

# Simply call the `enable_cloud_dataset` method from `Observations`.
# The default provider is `AWS`, but we will write it in manually for this example:
Observations.enable_cloud_dataset(provider='AWS')

# Getting the cloud URIs
obs_table = Observations.query_criteria(proposal_id=['2079'],
                                        instrument_name=['NIRISS/WFSS'],
                                        dataproduct_type=['spectrum'])
products = Observations.get_product_list(obs_table)
filtered = Observations.filter_products(products, productFilename=['jw02079-o001_s000000680_niriss_f150w-gr150c_x1d.fits'])
s3_uris = Observations.get_cloud_uris(filtered)
print(s3_uris)

# Download the data locally if you want
# Observations.download_products(filtered, cloud_only=True)
# will save to ./mastDownload/JWST/jw02079-o001_s000000680_niriss_f150w-gr150c/jw02079-o001_s000000680_niriss_f150w-gr150c_x1d.fits

Observations.disable_cloud_dataset()

INFO: Using the S3 STScI public dataset [astroquery.mast.cloud]

['s3://stpubdata/jwst/public/jw02079/L3/s/o001/jw02079-o001_s000000680_niriss_f150w-gr150c_x1d.fits']

# Open file in readonly mode and get the data
with fits.open(s3_uris[0], 'readonly', fsspec_kwargs={"anon": True}) as hdu:
    hdu.info()
    # Create Spectrum1D object (will not be needed when Jdaviz can open file directly from S3)
    wave = hdu[1].data['WAVELENGTH'] * u.Unit(hdu[1].header['TUNIT1'])
    flux = hdu[1].data['FLUX'] * u.Unit(hdu[1].header['TUNIT2'])
    std = StdDevUncertainty(hdu[1].data['FLUX_ERROR'] * u.Unit(hdu[1].header['TUNIT3']))
    spec1d = Spectrum1D(spectral_axis=wave,
                        flux=flux,
                        uncertainty=std)

Filename: <class 's3fs.core.S3File'>
No.    Name      Ver    Type      Cards   Dimensions   Format
  0  PRIMARY       1 PrimaryHDU     355   ()      
  1  EXTRACT1D     1 BinTableHDU     80   114R x 18C   [D, D, D, D, D, D, D, D, D, D, D, J, D, D, D, D, D, D]   
  2  EXTRACT1D     2 BinTableHDU     80   115R x 18C   [D, D, D, D, D, D, D, D, D, D, D, J, D, D, D, D, D, D]   
  3  EXTRACT1D     3 BinTableHDU     80   114R x 18C   [D, D, D, D, D, D, D, D, D, D, D, J, D, D, D, D, D, D]   
  4  ASDF          1 BinTableHDU     11   1R x 1C   [29963B]   

2. Open Specviz and load the 1D spectra we are interested in#

We use the Spectrum1D object we created until Specviz can stream directly from AWS.

viz = Specviz()
viz.load_data(spec1d, data_label='NIRISS F150W')
viz.show()

# Some zoom to see it better
plg_plot = viz.plugins['Plot Options']
plg_plot.y_min = -1E-6
plg_plot.y_max = 2E-5

3. Select the emission lines using the GUI and in the notebook#

I select the region spanning the emission lines from roughly 1.58 to 1.63 microns.

Instructions: https://jdaviz.readthedocs.io/en/latest/specviz/displaying.html#defining-spectral-regions

collection = viz.data_labels
print(collection)
dataout = viz.get_data('NIRISS F150W')
dataout

['NIRISS F150W']

<Spectrum1D(flux=[nan ... nan] Jy (shape=(114,), mean=0.00002 Jy); spectral_axis=<SpectralAxis [1.76217806 1.75747955 1.75278091 ... 1.24063718 1.23593867 1.23124003] um> (length=114); uncertainty=StdDevUncertainty)>

try:
    plg_subsets = viz.plugins['Subset Tools']
    region = plg_subsets.get_regions()
    print(region['Subset 1'])
except KeyError:
    print("No region defined in the GUI")

No region defined in the GUI

cont_region = SpectralRegion(1.598*u.um, 1.621*u.um)
spec1d_el_code = extract_region(dataout, cont_region)
print(spec1d_el_code)

Spectrum1D (length=4)
Flux=[6.18657876e-06 1.17189457e-05 1.24067720e-05 7.14436535e-06] Jy,  mean=0.00001 Jy
Spectral Axis=[1.61652279 1.61182415 1.60712564 1.60242712] um,  mean=1.60947 um
Uncertainty=StdDevUncertainty ([1.10100516e-07 1.18055113e-07 1.18216998e-07 1.09990686e-07] Jy)

try:
    spec1d_el_viz = extract_region(dataout, region['Subset 1'])
    print(spec1d_el_viz)
except KeyError:
    print("Region was not defined in the GUI")
    # Define spec1d_el_viz as spec1d_el_code
    spec1d_el_viz = spec1d_el_code

Region was not defined in the GUI

Plot the spectrum and the subset with matplotlib#

plt.plot(dataout.spectral_axis, dataout.flux, label='data')
plt.plot(spec1d_el_viz.spectral_axis, spec1d_el_viz.flux, label='subset defined in tool')
plt.plot(spec1d_el_code.spectral_axis, spec1d_el_code.flux, label='subset defined in code')
plt.legend()
plt.xlabel("wavelength ({:latex})".format(dataout.spectral_axis.unit))
plt.ylabel("flux ({:latex})".format(dataout.flux.unit))
plt.ylim(-1E-6, 2E-5)
plt.title("NIRISS F150W")
plt.show()

4. Use the redshift slider in Specviz to find the redshift#

I can use an available line lists or define my own lines (I know I need Hb4862.68 and the [OIII]4960.29,5008.24 doublet) and play with the redshift slider to match the lines in the line list with the lines in the spectrum. The line list plugin can be found clicking the plugin icon on the upper right of the viewer. To input just the three lines, I can use the “Custom” menu.

Here is the documentation where line lists are explained: https://jdaviz.readthedocs.io/en/latest/specviz/plugins.html#line-lists

I can also define the lines of interest programmatically, as shown in the following cell.

lt = QTable()
lt['linename'] = ['Hb', '[OIII]1', '[OIII]2']
lt['rest'] = [4862.68, 4960.29, 5008.24]*u.AA

viz.load_line_list(lt)

The lines are not showing now because their rest value is outside the range plotted here. I can move the lines using the redshift slider in the line list plugin. It is best to first set the redshift to 2 in the box with the number and then move the slider to bring the lines on top of the observed emission lines.

spec1d_redshift = viz.get_spectra(apply_slider_redshift=True)["NIRISS F150W"]
print(spec1d_redshift)

print()

if spec1d_redshift.redshift != 0.0:
    print("NIRISS F150W redshift=", spec1d_redshift.redshift)
else:
    print("Redshift was not defined in GUI. Defining it here.")
    spec1d_redshift.set_redshift_to(2.2138)
    print("NIRISS F150W redshift=", spec1d_redshift.redshift)

Spectrum1D (length=114)
Flux=[nan nan nan ... nan nan nan] Jy,  mean=0.00002 Jy
Spectral Axis=[1.76217806 1.75747955 1.75278091 ... 1.24063718 1.23593867
               1.23124003] um,  mean=1.49671 um
Uncertainty=StdDevUncertainty ([0. 0. 0. ... 0. 0. 0.] Jy)

Redshift was not defined in GUI. Defining it here.
NIRISS F150W redshift= 2.2138000000000004

5. Model the continuum of the spectrum#

I can use the GUI to select the region where I see the continuum. Challenge: select a discontinuous subset that covers two intervals (1.35-1.55um and 1.63-1.65um). Hint: select “Add” at the top near the Subset dropdown.

I can then use the Model Fitting plugin under the plugin icon to fit a linear model to the selected region. Instructions can be found here: https://jdaviz.readthedocs.io/en/latest/specviz/plugins.html#model-fitting. The individual steps to complete this task are:

• Select Subset 1 under Data

• Select Linear1D or Polynomial (whatever you think is best) under Model

• Click Add Component

• Enter a name for the model under Output Data Label (I choose “continuum”)

• Click Fit Model

I can extract the model and its parameters from the datasets in use.

datalabels = viz.data_labels
print(datalabels)
if 'continuum' in datalabels:
    spectrum = viz.get_data('NIRISS F150W')
    continuum = viz.get_data('continuum')
    print(continuum)
else:
    print("Continuum has not been created. Setting it to 0")
    spectrum = viz.get_data('NIRISS F150W')
    continuum = Spectrum1D(spectral_axis=spectrum.spectral_axis, flux=0.*spectrum.flux)

['NIRISS F150W']
Continuum has not been created. Setting it to 0

I can do a continuum subtraction and plot the result with matplotlib. If the continuum has not been defined in the GUI, this operation returns the original spectrum unchanged.

spectrum_sub = spectrum - continuum

plt.plot(spectrum_sub.spectral_axis, spectrum_sub.flux)
plt.hlines(0, 1.3, 1.7, color='black')
plt.xlabel("wavelength ({:latex})".format(spectrum_sub.spectral_axis.unit))
plt.ylabel("flux ({:latex})".format(spectrum_sub.flux.unit))
plt.title("NIRISS F150W")
plt.ylim(-1E-6, 2E-5)
plt.show()