{ "cells": [ { "cell_type": "markdown", "id": "strong-airfare", "metadata": {}, "source": [ "\n", "# NIRSpec IFU pipeline processing" ] }, { "cell_type": "markdown", "id": "fluid-importance", "metadata": {}, "source": [ "**Author**: James Muzerolle\n", "\n", "Plotting function originally developed by Bryan Hilbert\n", "\n", "**Latest Update**: June 10, 2022" ] }, { "cell_type": "markdown", "id": "basic-married", "metadata": {}, "source": [ "## Table of Contents\n", "* [Introduction](#intro)\n", " * [Overview](#overview)\n", " * [Simulated data](#sims)\n", "* [Imports](#imports)\n", "* [Convenience functions](#func)\n", "* [Pipeline processing flag](#flag)\n", "* [Download data](#data)\n", "* [Input data](#inputs)\n", "* [Level 2 association files](#lvl2asn)\n", "* [Run the calwebb_spec2 pipeline](#runspec2)\n", "* [Run individual steps of the spec2 pipeline](#runspec2steps)\n", " * [assign_wcs](#awcs)\n", " * [background](#background)\n", " * [imprint](#imprint)\n", " * [sourcetype](#srctype)\n", " * [flat_field](#flat)\n", " * [pathloss](#pathloss)\n", " * [photom](#photom)\n", " * [cube_build](#cubebuild)\n", " * [extract_1d](#extract1d)\n", "* [Run the calwebb_spec3 pipeline](#runspec3)" ] }, { "cell_type": "markdown", "id": "retained-driving", "metadata": {}, "source": [ "## Introduction " ] }, { "cell_type": "markdown", "id": "current-alexander", "metadata": {}, "source": [ "### Overview " ] }, { "cell_type": "markdown", "id": "express-declaration", "metadata": {}, "source": [ "In this notebook, we will explore the stage 2 and 3 pipelines for NIRSpec IFU data. Here, we will focus on the mechanics of processing \"real\" example data, including how to use associations for background subtraction and multi-exposure combination, what particular steps actually do to the data, and what the primary data products at each stage look like. We will also briefly examine how to interact and work with data models for each product.\n", "\n", "We are using pipeline version 1.1.0 for all data processing in this notebook. Most of the processing runs shown here use the default reference files from the Calibration Reference Data System (CRDS), with one exception at the end to show an example of how to modify/override. Please note that pipeline software development is a continuous process, so results in some cases may be slightly different if using a subsequent version. There are also a few known issues with some of the pipeline steps in this build that are expected to be fixed in the near future, though these do not significantly effect the products you will see here. Finally, some of the calibration reference files used by individual pipeline steps in the current CRDS context are placeholders, as they require calibrations that can only be taken in flight; for this reason, the absolute flux values seen here should not be taken literally." ] }, { "cell_type": "markdown", "id": "small-level", "metadata": {}, "source": [ "### Simulated data " ] }, { "cell_type": "markdown", "id": "brown-tackle", "metadata": {}, "source": [ "We will be using simulated NIRSpec exposures generated by the ESA Instrument Performance Simulator (IPS), using an input scene consisting of a single point source within the IFU aperture described by an emission-line galaxy template spectrum. The observation consists of four dithered exposures taken using the first four positions in the IFU cycling pattern. The instrument was configured with the G140H+F100LP grating/filter combination. The dithers are numbered for convenience in these example file names according to the ordering of the positions in the pattern. These file names are not indicative of actual data product file names you will see in the archive.\n", "\n", "There are a number of caveats to be aware of regarding these simulated data. 1) The IPS does not include a full treatment of all of the effects corrected by the stage 2 pipeline, particularly some of the throughput components. As with the above caveat regarding reference files, the simulations shown here should not be used for any analyses of flux information. 2) The simulated PSF is truncated in order to save on computation time, which results in an artifical drop in apparent count rate in the PSF wings. 3) Spacecraft pointing-related information has to be added by-hand to the headers before the simulated data can be processed by the pipeline. These keywords are used by the stage 3 pipeline when combining exposures, in order to align the spectral traces. Because this has to be a manual process and may be subject to small errors, the quality of the combined products here should not be taken as indicative of actual in-flight performance." ] }, { "cell_type": "markdown", "id": "moral-french", "metadata": {}, "source": [ "## Imports " ] }, { "cell_type": "markdown", "id": "departmental-metropolitan", "metadata": {}, "source": [ "Import packages necessary for this notebook" ] }, { "cell_type": "code", "execution_count": null, "id": "portuguese-harvey", "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "\n", "import glob\n", "\n", "#Modify the path to a directory on your machine\n", "import os\n", "os.environ[\"CRDS_PATH\"] = \"/path/to/my/folder/\"\n", "os.environ[\"CRDS_SERVER_URL\"] = \"https://jwst-crds-pub.stsci.edu\"\n", "\n", "import zipfile\n", "import urllib.request\n", "\n", "import json\n", "\n", "from astropy.io import fits\n", "from astropy.utils.data import download_file\n", "import astropy.units as u\n", "from astropy import wcs\n", "from astropy.wcs import WCS\n", "from astropy.visualization import ImageNormalize, ManualInterval, LogStretch, LinearStretch, AsinhStretch" ] }, { "cell_type": "markdown", "id": "manual-tackle", "metadata": {}, "source": [ "Set up matplotlib for plotting" ] }, { "cell_type": "code", "execution_count": null, "id": "realistic-announcement", "metadata": {}, "outputs": [], "source": [ "import matplotlib.pyplot as plt\n", "import matplotlib as mpl\n", "\n", "# Use this version for non-interactive plots (easier scrolling of the notebook)\n", "%matplotlib inline\n", "\n", "# Use this version (outside of Jupyter Lab) if you want interactive plots\n", "#%matplotlib notebook\n", "\n", "# These gymnastics are needed to make the sizes of the figures\n", "# be the same in both the inline and notebook versions\n", "%config InlineBackend.print_figure_kwargs = {'bbox_inches': None}\n", "\n", "mpl.rcParams['savefig.dpi'] = 80\n", "mpl.rcParams['figure.dpi'] = 80\n", "\n", "plt.rcParams.update({'font.size': 18})" ] }, { "cell_type": "markdown", "id": "opposed-identification", "metadata": {}, "source": [ "Import JWST pipeline modules" ] }, { "cell_type": "code", "execution_count": null, "id": "pregnant-briefing", "metadata": {}, "outputs": [], "source": [ "# The calwebb_spec and spec3 pipelines\n", "from jwst.pipeline import Spec2Pipeline\n", "from jwst.pipeline import Spec3Pipeline\n", "\n", "# individual steps\n", "from jwst.assign_wcs import AssignWcsStep\n", "from jwst.assign_wcs import nirspec\n", "from jwst.background import BackgroundStep\n", "from jwst.imprint import ImprintStep\n", "from jwst.msaflagopen import MSAFlagOpenStep\n", "from jwst.extract_2d import Extract2dStep\n", "from jwst.srctype import SourceTypeStep\n", "from jwst.wavecorr import WavecorrStep\n", "from jwst.flatfield import FlatFieldStep\n", "from jwst.pathloss import PathLossStep\n", "from jwst.photom import PhotomStep\n", "from jwst.cube_build import CubeBuildStep\n", "from jwst.extract_1d import Extract1dStep\n", "\n", "# data models\n", "from jwst import datamodels" ] }, { "cell_type": "markdown", "id": "smart-index", "metadata": {}, "source": [ "## Define convenience functions and parameters " ] }, { "cell_type": "code", "execution_count": null, "id": "passive-berkeley", "metadata": {}, "outputs": [], "source": [ "# All files created by the steps in this notebook have been pre-computed and cached in the following directory\n", "# if you want to actually run everything offline, then comment out this line and specify a desired local directory\n", "#output_dir = '/home/shared/preloaded-fits/nirspec-files/'\n", "output_dir = './nirspec_files/'\n", "if not os.path.exists(output_dir):\n", " os.makedirs(output_dir)" ] }, { "cell_type": "code", "execution_count": null, "id": "registered-adobe", "metadata": {}, "outputs": [], "source": [ "def show_image(data_2d, vmin, vmax, xsize=15, ysize=15, title=None, aspect=1, scale='log', units='MJy/sr'):\n", " \"\"\"Function to generate a 2D, log-scaled image of the data\n", " \n", " Parameters\n", " ----------\n", " data_2d : numpy.ndarray\n", " 2D image to be displayed\n", " \n", " vmin : float\n", " Minimum signal value to use for scaling\n", " \n", " vmax : float\n", " Maximum signal value to use for scaling\n", " \n", " title : str\n", " String to use for the plot title\n", " \n", " scale : str\n", " Specify scaling of the image. Can be 'log' or 'linear'\n", " \n", " units : str\n", " Units of the data. Used for the annotation in the\n", " color bar\n", " \"\"\"\n", " if scale == 'log':\n", " norm = ImageNormalize(data_2d, interval=ManualInterval(vmin=vmin, vmax=vmax),\n", " stretch=LogStretch())\n", " elif scale == 'linear':\n", " norm = ImageNormalize(data_2d, interval=ManualInterval(vmin=vmin, vmax=vmax),\n", " stretch=LinearStretch())\n", " elif scale == 'Asinh':\n", " norm = ImageNormalize(data_2d, interval=ManualInterval(vmin=vmin, vmax=vmax),\n", " stretch=AsinhStretch())\n", " fig = plt.figure(figsize=(xsize, ysize))\n", " ax = fig.add_subplot(1, 1, 1)\n", " im = ax.imshow(data_2d, origin='lower', norm=norm, aspect=aspect, cmap='gist_earth')\n", "\n", " fig.colorbar(im, label=units)\n", " plt.xlabel('Pixel column')\n", " plt.ylabel('Pixel row')\n", " if title:\n", " plt.title(title)" ] }, { "cell_type": "markdown", "id": "cardiovascular-plastic", "metadata": {}, "source": [ "## Pipeline processing flag " ] }, { "cell_type": "markdown", "id": "baking-flash", "metadata": {}, "source": [ "The pipeline and individual steps take too long to run in real time for this demo, so all products shown here have been pre-computed, and the actual pipeline calls will be skipped. Change the following flag to True if you want to run everything offline." ] }, { "cell_type": "code", "execution_count": null, "id": "strong-basic", "metadata": {}, "outputs": [], "source": [ "runflag = True" ] }, { "cell_type": "markdown", "id": "friendly-diamond", "metadata": {}, "source": [ "## Download Data " ] }, { "cell_type": "markdown", "id": "closed-baking", "metadata": {}, "source": [ "The input simulated data and all pipeline products are available on Box. If you want to run the pipeline offline, first download the zip file and unpack in the following cells; just comment out the first line in each." ] }, { "cell_type": "code", "execution_count": null, "id": "useful-bridge", "metadata": {}, "outputs": [], "source": [ "#%%script false --no-raise-error\n", "\n", "# set the Box link and file name\n", "ziplink = 'https://stsci.box.com/shared/static/c8uu9cn0dldkjv86z5rz16rszb39gjzi.zip'\n", "zipfilename = 'nirspec_ifudata.zip'\n", "if not os.path.isfile(os.path.join(output_dir, zipfilename)):\n", " print('Downloading {}...'.format(zipfilename))\n", " demo_file = download_file(ziplink, cache=True)\n", " # Make a symbolic link using a local name for convenience\n", " os.symlink(demo_file, os.path.join(output_dir, zipfilename))\n", "else:\n", " print('{} already exists, skipping download...'.format(zipfilename))" ] }, { "cell_type": "code", "execution_count": null, "id": "pacific-creature", "metadata": {}, "outputs": [], "source": [ "#%%script false --no-raise-error\n", "\n", "# unzip\n", "zf = zipfile.ZipFile(output_dir+'nirspec_ifudata.zip', 'r')\n", "zf.extractall(output_dir)" ] }, { "cell_type": "markdown", "id": "corporate-election", "metadata": {}, "source": [ "## Input simulations " ] }, { "cell_type": "markdown", "id": "labeled-recognition", "metadata": {}, "source": [ "We will be using simulated NIRSpec exposures of a point source observed with the IFU. Because the simulator generates count rate maps, equivalent to level 2a data products, we have to skip stage 1 of the pipeline and instead start the processing with the calwebb_spec2 pipeline.\n", "First, let's take a look at a few of the level 2a images to get familiarized with the inputs. The observation consists of 4 dithered exposures taken with a 4-point dither pattern." ] }, { "cell_type": "code", "execution_count": null, "id": "solved-constant", "metadata": { "scrolled": false }, "outputs": [], "source": [ "# get the data model of dither position 1:\n", "ratefile1 = output_dir+'nirspec_ifusim_d1_nrs1_rate.fits' # for the NRS1 detector\n", "dither1 = datamodels.open(ratefile1)\n", "ratefile2 = output_dir+'nirspec_ifusim_d1_nrs2_rate.fits' # for the NRS2 detector\n", "dither2 = datamodels.open(ratefile2)\n", "\n", "# get the pixel data (the SCI extension of the fits file)\n", "ratesci1 = dither1.data\n", "ratesci2 = dither2.data\n", "\n", "# display the images\n", "show_image(ratesci1, 0, 5.e2, units='DN/s')\n", "show_image(ratesci2, 0, 5.e2, units='DN/s')" ] }, { "cell_type": "markdown", "id": "ranking-climate", "metadata": {}, "source": [ "## Level 2 association files " ] }, { "cell_type": "markdown", "id": "recovered-nowhere", "metadata": {}, "source": [ "The automated pipeline uses association files to specify inputs into the calwebb_spec2 pipeline. The primary use for these is to specify multiple exposures to be used as background for subtraction from a given science exposure. In this demo, we're working with a dither set that is optimized for extended sources, and not used for background subtraction, so we will only use the association file to specify the full set of exposures to be processed one-by-one." ] }, { "cell_type": "code", "execution_count": null, "id": "western-vegetarian", "metadata": {}, "outputs": [], "source": [ "# show the contents of the association file\n", "asn_file = output_dir+\"l2_asn.json\"\n", "with open(asn_file) as f_obj:\n", " asn_data = json.load(f_obj)\n", "asn_data" ] }, { "cell_type": "markdown", "id": "ambient-compact", "metadata": {}, "source": [ "## Run the calwebb_spec2 pipeline " ] }, { "cell_type": "markdown", "id": "circular-basement", "metadata": {}, "source": [ "The following cell demonstrates how to run the calwebb_spec2 pipeline in full, using the above association file as input. Since this takes several hours to run in real time, the cell has been disabled for the webinar. In order to run it offline, simply delete or comment out the top line." ] }, { "cell_type": "code", "execution_count": null, "id": "streaming-rates", "metadata": {}, "outputs": [], "source": [ "if (runflag == True):\n", " spec2 = Spec2Pipeline()\n", " spec2.save_results = True\n", " spec2.output_dir = output_dir\n", " # skip the flat field correction, since the simulations do not include a full treatment of the throughput\n", " spec2.flat_field.skip = True\n", " result = spec2(asn_file)" ] }, { "cell_type": "code", "execution_count": null, "id": "considerable-tournament", "metadata": { "scrolled": true }, "outputs": [], "source": [ "# take a look at the results - open the level 2b files\n", "\n", "callist = [f for f in glob.glob(output_dir+\"*_cal.fits\")]\n", "callist.sort()\n", "for calfile in callist:\n", " # the unrectified 2D calibrated image\n", " print(calfile)\n", " cal = datamodels.MultiSpecModel(calfile) # this contains the calibrated unrectified 2D spectrum\n", " calsci = cal.data\n", " root = calfile[: -9]\n", " \n", " # the exposure-level cube\n", " hdu = fits.open(root+'_s3d.fits')\n", " cube = hdu[1].data\n", " wcs = WCS(hdu[1].header)\n", " hdu.close()\n", " \n", " # 1D spectrum extracted from the exposure-level cube\n", " x1d = datamodels.MultiSlitModel(root+'_x1d.fits') # this contains the aperture-extracted 1D spectrum\n", " \n", " # get wavelength & flux from the x1d data model\n", " x1dwave = x1d.spec[0].spec_table['WAVELENGTH']\n", " x1dflux = x1d.spec[0].spec_table['FLUX']\n", " \n", " # plot the 2D image \n", " show_image(calsci, 0, 1.5e2)\n", "\n", " # plot a slice from the cube\n", " # get slice wavelength values from the WCS transform for plot labels\n", " nslice = 1000\n", " sky, wave = wcs.pixel_to_world(0,0,nslice)\n", " title = np.around(wave.to(u.um),3)\n", " ax = plt.subplot(1, 1, 1, projection=wcs, slices=('x', 'y', nslice))\n", " norm=ImageNormalize(cube[nslice, :, :], vmin=0, vmax=5.e2, stretch=LogStretch())\n", " ax.imshow(cube[nslice, :, :], norm=norm, origin='lower', cmap='gist_earth')\n", " ax.set_xlabel('RA')\n", " ax.set_ylabel('DEC')\n", " ax.grid(color='white', ls='solid')\n", " ax.set_title(title)\n", " \n", " fig = plt.figure(figsize=(19, 8))\n", " plt.plot(x1dwave, x1dflux)\n", " plt.xlabel('wavelength (microns)')\n", " plt.ylabel('flux (Jy)')\n", " plt.show()" ] }, { "cell_type": "markdown", "id": "neither-course", "metadata": {}, "source": [ "## Run individual steps of the spec2 pipeline " ] }, { "cell_type": "markdown", "id": "secret-plain", "metadata": {}, "source": [ "Now let's try running individual pipeline steps on a single exposure and single detector, and take a look at their output to give a flavor for what they are doing to the data." ] }, { "cell_type": "markdown", "id": "supposed-wonder", "metadata": {}, "source": [ "### assign_wcs " ] }, { "cell_type": "markdown", "id": "sealed-broadway", "metadata": {}, "source": [ "This is a complex step that generates tranfsorms between various instrument optical planes and the sky, using the NIRSpec instrument model with multiple reference files." ] }, { "cell_type": "code", "execution_count": null, "id": "sorted-husband", "metadata": {}, "outputs": [], "source": [ "filename = output_dir+'nirspec_ifusim_d1_nrs1_rate.fits' # an input level 2a file\n", "root = filename[:-9] # for later matching\n", "print(root)\n", "\n", "if (runflag == True):\n", " step = AssignWcsStep()\n", " step.save_results = True\n", " step.output_dir = output_dir\n", " result = step(filename)" ] }, { "cell_type": "markdown", "id": "linear-elements", "metadata": {}, "source": [ "Note the verbose output prints the values of various instrument parameters that are relevant for the construction of the instrument transforms, such as the grating wheel tilt and the slit (or slits, in the case of MOS) used for the exposure." ] }, { "cell_type": "code", "execution_count": null, "id": "heavy-vitamin", "metadata": {}, "outputs": [], "source": [ "# The output file has the suffix _assignwcsstep appended:\n", "awcs_output_file = root+'assignwcsstep.fits'\n", "print(awcs_output_file)\n", "\n", "# load the output into a data model container\n", "awcs = datamodels.open(awcs_output_file)\n", "awcs" ] }, { "cell_type": "markdown", "id": "blessed-hawaii", "metadata": {}, "source": [ "The actual pixel data are not changed as a result of this step. Instead, a WCS object is created in the asdf extension of the output file that includes all the various optical transforms and component descriptions encompassed by the instrument model." ] }, { "cell_type": "code", "execution_count": null, "id": "undefined-fabric", "metadata": {}, "outputs": [], "source": [ "# get the WCS information populated by the algorithms of the assign_wcs step\n", "wcsobj = awcs.meta.wcs" ] }, { "cell_type": "code", "execution_count": null, "id": "independent-processor", "metadata": {}, "outputs": [], "source": [ "# list the frames of reference available for transformation\n", "wcsobj.available_frames" ] }, { "cell_type": "code", "execution_count": null, "id": "cognitive-investing", "metadata": {}, "outputs": [], "source": [ "# examples of the units for some of these frames\n", "print(wcsobj.detector.unit) # xy pixel indices\n", "print(wcsobj.slit_frame.unit) # relative xy position in the slit aperture\n", "print(wcsobj.msa_frame.unit) # absolute xy position in the slit aperture\n", "print(wcsobj.world.unit) # RA, Dec sky position" ] }, { "cell_type": "markdown", "id": "alpine-trademark", "metadata": {}, "source": [ "These correspond to image and pupil planes in the instrument, including the detector, the grating wheel, the slit aperture, and various sky planes. One can transform between any two of these using the appropriate transformation. Let's try an example of that to plot the wavelength of an emission line on the 2D spectrum." ] }, { "cell_type": "code", "execution_count": null, "id": "changed-player", "metadata": {}, "outputs": [], "source": [ "# specify the wavelength of a feature of interest\n", "wave0 = 1.2518\n", "\n", "# pick one of the IFU slices\n", "ifuslice = nirspec.nrs_wcs_set_input(awcs, 4) # the source is centered in slice 2 for dither 1\n", "y, x = np.mgrid[:awcs.data.shape[0], : awcs.data.shape[1]]\n", "ra, dec, wave = ifuslice(x, y)\n", "world2det = ifuslice.get_transform('world', 'detector')\n", "\n", "# get bounding box of the slice 2d trace\n", "xstart = int(np.ceil(ifuslice.bounding_box[0][0]))\n", "xstop = int(np.ceil(ifuslice.bounding_box[0][1]))\n", "ystart = int(np.ceil(ifuslice.bounding_box[1][0]))\n", "ystop = int(np.ceil(ifuslice.bounding_box[1][1]))\n", "imslice = awcs.data[ystart:ystop, xstart:xstop]\n", "waveslice = wave[ystart:ystop, xstart:xstop]\n", "\n", "# find the trace of the feature wavelength in 2D space\n", "xm = np.arange(imslice.shape[1])\n", "ym = np.arange(imslice.shape[0])\n", "xm_wave = np.zeros(waveslice.shape[0])\n", "for i in range(len(ym)):\n", " col = waveslice[i,:]\n", "\n", " if not np.all(np.isnan(col)):\n", " xm_wave[i] = np.interp(wave0, col[np.isfinite(col)], xm[np.isfinite(col)], left=0, right=0)\n", "\n", "show_image(imslice, 0., 2.e2, xsize=15, ysize=8, aspect=1., units='DN/s')\n", "plt.plot(xm_wave[xm_wave != 0], ym[xm_wave != 0], color='white', alpha=1, linewidth=2)\n", "plt.xlim(900,1250)\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "stable-brazilian", "metadata": {}, "source": [ "### background " ] }, { "cell_type": "markdown", "id": "operating-matter", "metadata": {}, "source": [ "If there are associated exposures to be used as backgrounds (as in a nod set), this step will average them (if more than one), and then directly subtract from the exposure being processed. In this demo, the example data are part of a cycling dither pattern, which is meant for improving the sampling of extended sources, and the dither separations are not large enough to enable image-to-image subtraction without source overlap. Thus, we will skip running this step (note, the full pipeline automatically skips it if there are no background exposures listed in the association file)." ] }, { "cell_type": "markdown", "id": "joined-alias", "metadata": {}, "source": [ "### imprint " ] }, { "cell_type": "markdown", "id": "affecting-finger", "metadata": {}, "source": [ "The MSA shutters are not perfectly opaque; a small amount of light from the sky leaks through the edges of each shutter. In addition, a very small number of shutters are stuck open and hence always open to sky. This results in a dispersed \"leakage\" signal (also referred to as \"imprint\") that projects onto the same area of the detectors that contain light from the IFU aperture. Dedicated leakage exposures, taken at the same position as the science exposure with the IFU aperture covered and the MSA all closed, can be obtained to remove this signal. The imprint step takes the associated leakage exposures and subtracts them pixel-by-pixel from the science exposure.\n", "\n", "The simulated exposures we are using for this demo do not contain leakage signal, so we will skip running this step (note, the full pipeline automatically skips it if there are no leakage exposures listed in the association file)." ] }, { "cell_type": "markdown", "id": "greatest-artist", "metadata": {}, "source": [ "### sourcetype " ] }, { "cell_type": "markdown", "id": "promising-shelf", "metadata": {}, "source": [ "This step sets the source type (point or extended) based on information passed from PPS (given in the SRCTYAPT header keyword); if no flag was set, by default the type is set to extended for IFU data. This affects subsequent processing steps, which do different things for each type." ] }, { "cell_type": "code", "execution_count": null, "id": "comparative-protest", "metadata": {}, "outputs": [], "source": [ "filename = root+'assignwcsstep.fits' # the output from the previous step\n", "\n", "if (runflag == True):\n", " step = SourceTypeStep()\n", " step.save_results = True\n", " step.output_dir = output_dir\n", " result = step(filename)" ] }, { "cell_type": "markdown", "id": "reliable-measure", "metadata": {}, "source": [ "keyword SRCTYPE should be set to \"POINT\" in this case, because PPS database value SRCTYAPT is \"POINT\"" ] }, { "cell_type": "code", "execution_count": null, "id": "chronic-lambda", "metadata": {}, "outputs": [], "source": [ "# The output file has the suffix _sourcetypestep appended:\n", "stype_output_file = root+'sourcetypestep.fits'\n", "\n", "# load the output into a data model container\n", "stype = datamodels.open(stype_output_file)\n", "\n", "# what is the SRCTYPE keyword value?\n", "stype.find_fits_keyword('SRCTYPE')\n", "print('SRCTYPE=', stype.meta.target.source_type)" ] }, { "cell_type": "markdown", "id": "differential-discipline", "metadata": {}, "source": [ "### flat_field " ] }, { "cell_type": "markdown", "id": "elegant-refrigerator", "metadata": {}, "source": [ "A flat field correction is applied to each pixel in the 2D spectrum. The NIRSpec \"flat\" comprises three components: D flat, which includes the pixel-to-pixel detector response; S flat, which includes the throughput of the spectrograph; F flat, which folds together the throughput of the filter, FORE optics, and OTE." ] }, { "cell_type": "code", "execution_count": null, "id": "intense-green", "metadata": {}, "outputs": [], "source": [ "filename = root+'sourcetypestep.fits' # the output from the previous step\n", "\n", "if (runflag == True):\n", " step = FlatFieldStep()\n", " step.save_interpolated_flat = True # this will save the on-the-fly flat field correction values as an image\n", " step.save_results = True\n", " step.output_dir = output_dir\n", " result = step(filename)" ] }, { "cell_type": "code", "execution_count": null, "id": "demographic-albert", "metadata": {}, "outputs": [], "source": [ "# The output file has the suffix _flatfieldstep appended:\n", "flat_output_file = root+'flatfieldstep.fits'\n", "\n", "# The optional saved flat correction image has the suffix _interpolatedflat appended:\n", "intflat_output_file = root+'interpolatedflat.fits'\n", "\n", "# load the output into a data model container\n", "flat = datamodels.open(flat_output_file)\n", "intflat = datamodels.open(intflat_output_file)" ] }, { "cell_type": "code", "execution_count": null, "id": "arranged-humor", "metadata": {}, "outputs": [], "source": [ "# plot the flat correction image - this is divided into the science exposure to apply the corrections per pixel\n", "show_image(intflat.data, 0.3, 1., scale='Asinh', units='DN/s')" ] }, { "cell_type": "code", "execution_count": null, "id": "federal-things", "metadata": {}, "outputs": [], "source": [ "# plot the corrected science exposure\n", "# note the simulation does not include all of the flat components, so applying the full correction in this case produces incorrect results\n", "show_image(flat.data, 0., 8.e1, scale='Asinh', units='DN/s')" ] }, { "cell_type": "markdown", "id": "underlying-taste", "metadata": {}, "source": [ "### pathloss " ] }, { "cell_type": "markdown", "id": "seeing-centre", "metadata": {}, "source": [ "The pathloss correction scales the 2D spectrum as a function of wavelength to account for diffraction losses incurred by the grating wheel. These will not be measured until flight, so the current reference file is only a placeholder and the correction factors are unity." ] }, { "cell_type": "code", "execution_count": null, "id": "noticed-helen", "metadata": {}, "outputs": [], "source": [ "filename = root+'flatfieldstep.fits' # the output from the previous step\n", "\n", "if (runflag == True):\n", " step = PathLossStep()\n", " step.save_results = True\n", " step.output_dir = output_dir\n", " result = step(filename)" ] }, { "cell_type": "code", "execution_count": null, "id": "interpreted-semester", "metadata": {}, "outputs": [], "source": [ "# The output file has the suffix _pathlossstep appended:\n", "ploss_output_file = root+'pathlossstep.fits'\n", "\n", "# load the output into a data model container\n", "ploss = datamodels.open(ploss_output_file)" ] }, { "cell_type": "code", "execution_count": null, "id": "opening-illinois", "metadata": {}, "outputs": [], "source": [ "# plot the corrected science exposure\n", "show_image(ploss.data, 0., 8.e1, scale='Asinh', units='DN/s')" ] }, { "cell_type": "markdown", "id": "thousand-lotus", "metadata": {}, "source": [ "### photom " ] }, { "cell_type": "markdown", "id": "short-difficulty", "metadata": {}, "source": [ "The photom step applies a scalar conversion factor to convert to physical units. Once on-orbit observations of spectrophotometric standards are obtained, this may also include a wavelength-dependent vector, if necessary, to account for any small discrepancies found in the throughput corrections." ] }, { "cell_type": "code", "execution_count": null, "id": "broad-calgary", "metadata": {}, "outputs": [], "source": [ "filename = root+'pathlossstep.fits' # the output from the previous step\n", "\n", "if (runflag == True):\n", " step = PhotomStep()\n", " step.save_results = True\n", " step.output_dir = output_dir\n", " result = step(filename)" ] }, { "cell_type": "code", "execution_count": null, "id": "pressing-belarus", "metadata": {}, "outputs": [], "source": [ "# The output file has the suffix _photomstep appended:\n", "phot_output_file = root+'photomstep.fits'\n", "\n", "# load the output into a data model container\n", "phot = datamodels.open(phot_output_file)" ] }, { "cell_type": "code", "execution_count": null, "id": "pointed-prairie", "metadata": {}, "outputs": [], "source": [ "# what are the flux calibration-related keywords?\n", "phot.find_fits_keyword('PHOTUJA2')\n", "print('PHOTMJSR=', phot.meta.photometry.conversion_megajanskys)\n", "# note that despite the keyword name, when SRCTYPE=POINT, the units are actually MJy per pixel\n", "# the pre-launch reference file has a placeholder value, so the output fluxes will appear unphysically large\n", "print('PHOTUJA2=', phot.meta.photometry.conversion_microjanskys)\n", "# this is only applicable for extended sources" ] }, { "cell_type": "markdown", "id": "accomplished-camping", "metadata": {}, "source": [ "### cube_build " ] }, { "cell_type": "markdown", "id": "vietnamese-venice", "metadata": {}, "source": [ "The fully calibrated 2D spectrum is resampled onto a rectified cube in sky and wavelength coordinates. This is intended for quick-look purposes only, as all of the calibrated unrectified exposures will be resampled and combined onto a single cube during the level 3 processing." ] }, { "cell_type": "code", "execution_count": null, "id": "adolescent-participant", "metadata": {}, "outputs": [], "source": [ "filename = root+'photomstep.fits' # the output from the previous step\n", "\n", "if (runflag == True):\n", " step = CubeBuildStep()\n", " step.save_results = True\n", " step.output_dir = output_dir\n", " result = step(filename)" ] }, { "cell_type": "code", "execution_count": null, "id": "hawaiian-updating", "metadata": {}, "outputs": [], "source": [ "# The output file has the suffix _photomstep appended:\n", "cube_output_file = root+'photomstep_g140h-f100lp_s3d.fits'\n", "print(cube_output_file)\n", "\n", "# read in the output fits file\n", "hdu = fits.open(cube_output_file)\n", "# get the exposure-level cube\n", "cube = hdu[1].data\n", "# get the WCS info\n", "wcs = WCS(hdu[1].header)\n", "cdelt = hdu[1].header['CDELT1']*3600.\n", "hdu.close()" ] }, { "cell_type": "code", "execution_count": null, "id": "promotional-property", "metadata": {}, "outputs": [], "source": [ "# show a slice\n", "\n", "nslice = 1750 # which wavelength slice of the cube to use for display & fitting\n", "# get slice wavelength values from the WCS transform for plot labels\n", "nslice = 1750\n", "sky, wave = wcs.pixel_to_world(0,0,nslice)\n", "title = np.around(wave.to(u.um),3)\n", "show_image(cube[nslice,:,:], 0., 1.e3, units='MJy')\n", "plt.xlabel('RA')\n", "plt.ylabel('DEC')\n", "plt.grid(color='white', ls='solid')\n", "plt.title(title)" ] }, { "cell_type": "markdown", "id": "rotary-estate", "metadata": {}, "source": [ "### extract_1d " ] }, { "cell_type": "markdown", "id": "acceptable-copyright", "metadata": {}, "source": [ "A 1D spectrum is extracted from the cube, using a circular aperture with radius specified from a reference file. An aperture correction is also applied to account for flux outside of the extraction aperture. As with the level 2b cube, this is intended for quick-look purposes only. An independent extraction of the combined cube will be done during level 3 processing." ] }, { "cell_type": "code", "execution_count": null, "id": "earned-bunny", "metadata": {}, "outputs": [], "source": [ "filename = root+'photomstep_g140h-f100lp_s3d.fits' # the output from the previous step\n", "\n", "if (runflag == True):\n", " step = Extract1dStep()\n", " step.save_results = True\n", " step.output_dir = output_dir\n", " result = step(filename)" ] }, { "cell_type": "code", "execution_count": null, "id": "protecting-onion", "metadata": {}, "outputs": [], "source": [ "# The output file has the suffix _extract1dstep appended:\n", "x1d_output_file = root+'photomstep_g140h-f100lp_extract1dstep.fits'\n", "\n", "# load the output into a data model container\n", "x1d = datamodels.open(x1d_output_file)" ] }, { "cell_type": "code", "execution_count": null, "id": "floral-hartford", "metadata": {}, "outputs": [], "source": [ "# plot the spectrum\n", "x1dwave = x1d.spec[0].spec_table.WAVELENGTH\n", "x1dflux = x1d.spec[0].spec_table.FLUX\n", "\n", "fig = plt.figure(figsize=(19, 8))\n", "plt.plot(x1dwave, x1dflux)\n", "plt.xlabel('wavelength (microns)')\n", "plt.ylabel('flux (Jy)')\n", "plt.show()\n", "\n", "# note that the flux level appears artificially inflated because of the dummy flux conversion factor applied in the photom step" ] }, { "cell_type": "markdown", "id": "clean-christopher", "metadata": {}, "source": [ "## Run the calwebb_spec3 pipeline " ] }, { "cell_type": "markdown", "id": "experimental-underwear", "metadata": {}, "source": [ "For an observation that contains multiple exposures of the same source, such as the dither set in this example, this pipeline will combine all of the exposures into a single output product. Both a 3D combined cube and extracted 1D spectrum are generated. This process also includes an outlier detection step that compares the stack of values at each resampled pixel and flags outliers based on noise threshold parameters; the flagged outliers are not included in the spectral combination. However, the current implementation of this step has not been properly tuned and tends to flag non-outliers associated with the source PSF, so we will skip it for this demo. We will also not demonstrate the optional \"master background\" step, which takes an independent 1D background spectrum to create a 2D background to subtract in the absence of nodded exposures; see the read-the-docs pages for more details." ] }, { "cell_type": "code", "execution_count": null, "id": "hungarian-candidate", "metadata": {}, "outputs": [], "source": [ "# show the contents of the association file\n", "asn_file = output_dir+\"l3_asn.json\"\n", "with open(asn_file) as f_obj:\n", " asn_data = json.load(f_obj)\n", "asn_data" ] }, { "cell_type": "markdown", "id": "handy-clearing", "metadata": {}, "source": [ "run the calwebb_spec3 pipeline using association file of cal.fits files" ] }, { "cell_type": "code", "execution_count": null, "id": "meaning-internship", "metadata": {}, "outputs": [], "source": [ "if (runflag == True):\n", " spec3 = Spec3Pipeline()\n", " spec3.save_results = True\n", " spec3.output_dir = output_dir\n", " spec3.outlier_detection.skip = True # skip this step for now, because the simulations do not include outliers such as CR hits or bad pixels\n", " result = spec3(asn_file)" ] }, { "cell_type": "code", "execution_count": null, "id": "eastern-scale", "metadata": { "scrolled": false }, "outputs": [], "source": [ "# display level 3 products\n", "# combined 3D cube\n", "hdu = fits.open(output_dir+'nirspec_ifusim_comb_1234_g140h-f100lp_s3d.fits')\n", "cube = hdu[1].data\n", "wcs = WCS(hdu[1].header)\n", "cdelt = hdu[1].header['CDELT1']*3600.\n", "hdu.close() \n", "\n", "# plot slices from the cube\n", "fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(17,8))\n", "\n", "# get slice wavelength values from the WCS transform for plot labels\n", "nslice1 = 1000\n", "nslice2 = 3500\n", "sky, wave = wcs.pixel_to_world(0, 0, nslice1)\n", "title1 = np.around(wave.to(u.um), 3)\n", "sky, wave = wcs.pixel_to_world(0, 0, nslice2)\n", "title2 = np.around(wave.to(u.um), 3)\n", "\n", "plt.rcParams.update({'font.size': 16})\n", "ax1 = plt.subplot(1, 2, 1, projection=wcs, slices=('x', 'y', nslice))\n", "norm=ImageNormalize(cube[nslice, :, :], vmin=0, vmax=5.e2, stretch=LogStretch())\n", "slice1 = ax1.imshow(cube[nslice, :, :], norm=norm, origin='lower', cmap='gist_earth')\n", "ax1.set_xlabel('RA')\n", "ax1.set_ylabel('DEC')\n", "ax1.grid(color='white', ls='solid')\n", "ax1.set_title(title1)\n", "\n", "ax2 = plt.subplot(1, 2, 2, projection=wcs, slices=('x', 'y', nslice))\n", "norm=ImageNormalize(cube[nslice, :, :], vmin=0, vmax=5.e2, stretch=LogStretch())\n", "slice2 = ax2.imshow(cube[nslice, :, :], norm=norm, origin='lower', cmap='gist_earth')\n", "ax2.set_xlabel('RA')\n", "ax2.set_ylabel('DEC')\n", "ax2.grid(color='white', ls='solid')\n", "ax2.set_title(title2)" ] }, { "cell_type": "code", "execution_count": null, "id": "responsible-consciousness", "metadata": {}, "outputs": [], "source": [ "# combined 1D extracted spectrum\n", "x1d3 = datamodels.open(output_dir+'nirspec_ifusim_comb_1234_g140h-f100lp_x1d.fits')\n", "\n", "# get wavelength & flux\n", "x1d3wave = x1d3.spec[0].spec_table.WAVELENGTH\n", "x1d3flux = x1d3.spec[0].spec_table.FLUX\n", "\n", "# plot\n", "fig = plt.figure(figsize=(15,9))\n", "plt.plot(x1d3wave,x1d3flux)\n", "plt.xlabel('wavelength (microns)')\n", "plt.ylabel('flux (Jy)')\n", "plt.show()" ] }, { "cell_type": "code", "execution_count": null, "id": "operating-estate", "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "Python 3 (ipykernel)", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.9.13" } }, "nbformat": 4, "nbformat_minor": 5 }