{ "cells": [ { "cell_type": "markdown", "id": "strong-airfare", "metadata": {}, "source": [ "\n", "# Running the spectroscopic pipeline: NIRSpec Fixed Slit" ] }, { "cell_type": "markdown", "id": "covered-hotel", "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", "* [Reference files](#ref)\n", "* [Download the 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", " * [extract_2d](#extract2d)\n", " * [sourcetype](#srctype)\n", " * [wavecorr](#wavecorr)\n", " * [flat_field](#flat)\n", " * [pathloss](#pathloss)\n", " * [photom](#photom)\n", " * [resample_spec](#resample)\n", " * [extract_1d](#extract1d)\n", "* [Run the calwebb_spec3 pipeline](#runspec3)\n", " * [Rerun extract_1d with a different extraction aperture](#rerun)" ] }, { "cell_type": "markdown", "id": "retained-driving", "metadata": {}, "source": [ "## Introduction " ] }, { "cell_type": "markdown", "id": "hundred-consortium", "metadata": {}, "source": [ "### Overview " ] }, { "cell_type": "markdown", "id": "express-declaration", "metadata": {}, "source": [ "In this notebook, we will explore the stage 2 and 3 pipelines for spectroscopic data, using the NIRSpec Fixed Slit mode as an example. Since the previous two notebooks described these pipelines and all their component steps for all spectroscopic modes, 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": "divided-latest", "metadata": {}, "source": [ "### Simulated data " ] }, { "cell_type": "markdown", "id": "removed-prior", "metadata": {}, "source": [ "We will be using simulated NIRSpec exposures generated by the ESA Instrument Performance Simulator (IPS), using an input scene consisting of a point source with an A-star stellar model spectrum. The observation consists of 9 dithered exposures with the S200A1 slit, including a 3-point primary dither pattern with spectral subdithers. The instrument was configured with the PRISM+CLEAR grating/filter combination and a detector subarray of SUBS200A1 (2048x64 pixels on each detector). The dithers are numbered for convenience in these example file names according to the source position in the slit, as shown schematically below. These file names are not indicative of actual data product file names you will see in the archive.\n", "\n", "1 2 3\n", "\n", "4 5 6\n", "\n", "7 8 9\n", "\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. In some cases, you may see 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": { "tags": [] }, "outputs": [], "source": [ "import numpy as np\n", "\n", "import glob\n", "\n", "import os\n", "# Update to a path in your system (see details below at \"Reference files\")\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", "import warnings\n", "\n", "from astropy.io import fits\n", "from astropy.utils.data import download_file\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" ] }, { "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.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.resample import ResampleSpecStep\n", "from jwst.extract_1d import Extract1dStep\n", "\n", "# data models\n", "from jwst import datamodels" ] }, { "cell_type": "markdown", "id": "relevant-austin", "metadata": {}, "source": [ "check the version" ] }, { "cell_type": "code", "execution_count": null, "id": "lonely-point", "metadata": {}, "outputs": [], "source": [ "import jwst\n", "print(jwst.__version__)" ] }, { "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 in this notebook will be saved to the current working directory\n", "# create a subdirectory to organize the 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, 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=(15, 9))\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": "interior-granny", "metadata": {}, "source": [ "## Reference files " ] }, { "cell_type": "markdown", "id": "incoming-cancer", "metadata": {}, "source": [ "Calibration reference files are a collection of FITS and ASDF files that are used to remove instrumental signatures and calibrate JWST data. For example, the dark current reference file contains a multiaccum ramp of dark current signal to be subtracted from the data during the dark current subtraction step.\n", "When running a pipeline or pipeline step, the pipeline will automatically look for any required reference files in a pre-defined local directory. If the required reference files are not present, they will automatically be downloaded from the Calibration Reference Data System (CRDS) at STScI.\n", "\n", "During the JWebbinar, our pre-existing existing environment is set up to correctly use and store calibration reference files, and you do not need to set the environment variables below." ] }, { "cell_type": "markdown", "id": "dynamic-coalition", "metadata": {}, "source": [ "If you wish to run this notebook outside of the JWebbinar, you will have to specify a local directory in which to store [reference files](https://jwst-pipeline.readthedocs.io/en/latest/jwst/introduction.html#crds), along with the server to use to download the reference files from CRDS. To accomplish this, there are two environment variables that should be set prior to importing the pipeline. These are the CRDS_PATH and CRDS_SERVER_URL variables. In the example below, reference files will be downloaded to the \"crds_cache\" directory under the home directory.\n", "\n", "`$ export CRDS_PATH=$HOME/crds_cache`\n", "\n", "`$ export CRDS_SERVER_URL=https://jwst-crds-pub.stsci.edu`\n", "\n", "OR:\n", "\n", "`os.environ[\"CRDS_PATH\"] = \"/user/myself/crds_cache\"`\n", "\n", "`os.environ[\"CRDS_SERVER_URL\"] = \"https://jwst-crds-pub.stsci.edu\"`\n", "\n", "The first time you run the pipeline, the CRDS server should download all of the context and reference files that are needed for that pipeline run, and dump them into the CRDS_PATH directory. Subsequent executions of the pipeline will first look to see if it has what it needs in CRDS_PATH and anything it doesn't have will be downloaded from the STScI cache.\n", "\n", "NOTE: Users at STScI should automatically have access to the Calibration Reference Data System (CRDS) cache for running the pipeline, and can skip setting these environment variables." ] }, { "cell_type": "markdown", "id": "reverse-meaning", "metadata": {}, "source": [ "## Download the data " ] }, { "cell_type": "markdown", "id": "consistent-lodge", "metadata": {}, "source": [ "First, let's grab the data from Box." ] }, { "cell_type": "code", "execution_count": null, "id": "enabling-dress", "metadata": {}, "outputs": [], "source": [ "# set the Box link and file name\n", "ziplink = 'https://stsci.box.com/shared/static/qdj79y7rv99wuui2hot0l3kg5ohq0ah9.zip'\n", "zipfilename = 'nirspec_data.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": "touched-creator", "metadata": {}, "outputs": [], "source": [ "# unzip\n", "zf = zipfile.ZipFile(output_dir+'nirspec_data.zip', 'r')\n", "zf.extractall(output_dir)" ] }, { "cell_type": "markdown", "id": "corporate-election", "metadata": {}, "source": [ "## Input data " ] }, { "cell_type": "markdown", "id": "labeled-recognition", "metadata": {}, "source": [ "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", "\n", "First, let's take a look at a few of the level 2a images to get familiarized with the inputs." ] }, { "cell_type": "code", "execution_count": null, "id": "solved-constant", "metadata": {}, "outputs": [], "source": [ "# get the data model of dither position 1:\n", "ratefile = output_dir+'nirspec_fssim_d1_rate.fits'\n", "dither = datamodels.open(ratefile)\n", "\n", "# get the pixel data (the SCI extension of the fits file)\n", "ratesci = dither.data\n", "\n", "# display the image\n", "show_image(ratesci, 0, 5.e2, units='DN/s')" ] }, { "cell_type": "code", "execution_count": null, "id": "nuclear-cooking", "metadata": {}, "outputs": [], "source": [ "# get the data model of dither position 4:\n", "ratefile = output_dir+'nirspec_fssim_d4_rate.fits'\n", "dither = datamodels.open(ratefile)\n", "\n", "# get the pixel data (the SCI extension of the fits file)\n", "ratesci = dither.data\n", "\n", "# display the image\n", "show_image(ratesci, 0, 5.e2, units='DN/s')" ] }, { "cell_type": "markdown", "id": "under-classics", "metadata": {}, "source": [ "Note how the source spectral trace moves \"down\" in the y (spatial) direction. Because the disperser used here is the PRISM, the length of the trace is much shorter than the ALLSLITS subarray width." ] }, { "cell_type": "markdown", "id": "ranking-climate", "metadata": {}, "source": [ "## Level 2 association files " ] }, { "cell_type": "markdown", "id": "recovered-nowhere", "metadata": {}, "source": [ "One purpose of the primary dither pattern in this case is to enable background subtraction. The spec2 pipeline is set up to handle this by using association files that list the exposures to be used as backgrounds for each input exposure. The background elements can also be specified manually as inputs into the \"background\" step, but we'll be using association files in this example." ] }, { "cell_type": "code", "execution_count": null, "id": "suited-lewis", "metadata": {}, "outputs": [], "source": [ "# get the association files\n", "asnlist = [f for f in glob.glob(output_dir+\"spec2*json\")]\n", "asnlist.sort()\n", "print(asnlist)" ] }, { "cell_type": "code", "execution_count": null, "id": "western-vegetarian", "metadata": {}, "outputs": [], "source": [ "# show the contents of one of the association files\n", "asn_file = asnlist[0]\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": "novel-director", "metadata": {}, "source": [ "The following cell demonstrates how to run the calwebb_spec2 pipeline in full, using the above association files as inputs. Since this takes a bit too long to run in real time during the webinar, the cell has been disabled. In order to run it offline, simply delete or comment out the top line. " ] }, { "cell_type": "code", "execution_count": null, "id": "handed-history", "metadata": {}, "outputs": [], "source": [ "%%script false --no-raise-error\n", "\n", "# loop over the association files (one per exposure)\n", "for asn in asnlist:\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)" ] }, { "cell_type": "code", "execution_count": null, "id": "considerable-tournament", "metadata": { "scrolled": true, "tags": [] }, "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", " print(calfile)\n", " cal = datamodels.open(calfile) # this contains the calibrated unrectified 2D spectrum\n", " root = calfile[:-9]\n", " s2d = datamodels.open(root+'_s2d.fits') # this contains the calibrated *rectified* 2D spectrum\n", " x1d = datamodels.MultiSlitModel(root+'_x1d.fits') # this contains the aperture-extracted 1D spectrum\n", " \n", " for slit in cal.slits:\n", " # this data model is set up to handle multiple slits, which is mainly needed for MOS and WFSS data\n", " # in this case, there is only one slit, S200A1\n", " print(slit.name)\n", " \n", " calsci = slit.data # contains the pixel data from the cal file (SCI extension)\n", " s2dsci = s2d.slits[0].data # contains the pixel data from the s2d file\n", " \n", " # determine the wavelength scale of the s2d data for plotting purposes\n", " # get the data model WCS object\n", " wcsobj = s2d.slits[0].meta.wcs\n", " y, x = np.mgrid[:s2dsci.shape[0], : s2dsci.shape[1]] # grid of pixel x,y indices\n", " det2sky = wcsobj.get_transform('detector', 'world') # the coordinate transform from detector space (pixels) to sky (RA, DEC in degrees)\n", " ra, dec, s2dwave = det2sky(x, y) # RA, Dec, wavelength (microns) for each pixel\n", " s2dwaves = s2dwave[0, :] # only need a single row of values since this is the rectified spectrum\n", " xtint = np.arange(100, s2dsci.shape[1], 100)\n", " xtlab = np.round(s2dwaves[xtint], 2) # wavelength labels for the x-axis\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", " show_image(calsci, -5.e2, 2.5e3, aspect=5., scale='linear', units='Jy')\n", " \n", " show_image(s2dsci, -5.e2, 2.5e3, aspect=5., scale='linear', units='Jy')\n", " plt.xticks(xtint, xtlab)\n", " plt.xlabel('wavelength (microns)')\n", " \n", " fig = plt.figure(figsize=(12, 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 taking 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_fssim_d5_rate.fits' # an input level 2a file\n", "root = filename[:-9] # for later matching\n", "print(root)\n", "\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. In this demo, a fixed slit was used, which is indicated here by \"quadrant 5\"; if this had been a MOS exposure, there would be multiple slits indicated in any or all of quadrants 1-4." ] }, { "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": "focal-switch", "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. We'll see an example of such a transform below, after running the extract_2d step." ] }, { "cell_type": "markdown", "id": "opened-floor", "metadata": {}, "source": [ "### background " ] }, { "cell_type": "markdown", "id": "intermediate-population", "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." ] }, { "cell_type": "code", "execution_count": null, "id": "wireless-queensland", "metadata": {}, "outputs": [], "source": [ "filename = root+'assignwcsstep.fits' # the output from the previous step\n", "\n", "backgrounds = [output_dir+'nirspec_fssim_d2_rate.fits', output_dir+'nirspec_fssim_d8_rate.fits']\n", "# these exposures are for the middle spectral dither positions at primary dither positions 1 and 3\n", "\n", "step = BackgroundStep()\n", "step.save_results = True\n", "step.output_dir = output_dir\n", "result = step(filename, backgrounds)" ] }, { "cell_type": "markdown", "id": "south-vector", "metadata": {}, "source": [ "Look at the output. Since this step does a direct image-to-image subtraction, and the input backgrounds were the \"top\" and \"bottom\" primary dither positions, we should see a positive source trace surrounded by negative traces on either side." ] }, { "cell_type": "code", "execution_count": null, "id": "analyzed-wagner", "metadata": {}, "outputs": [], "source": [ "# The output file has the suffix _backgroundstep appended:\n", "bsubfile = root+'backgroundstep.fits'\n", "\n", "# load the output into a data model container\n", "bsub = datamodels.open(bsubfile)\n", "\n", "# plot the image\n", "show_image(bsub.data, -5.e2, 25.e2, aspect=5., scale='linear', units='DN/s')" ] }, { "cell_type": "markdown", "id": "efficient-change", "metadata": {}, "source": [ "### extract_2d " ] }, { "cell_type": "markdown", "id": "incomplete-dialogue", "metadata": {}, "source": [ "Using the WCS transforms constructed by assign_wcs, this step extracts a 2D \"subwindow\" from the original full frame or (in this case) subarray image. For MOS data, there are multiple 2D extractions corresponding to the multiple sources observed in the exposure." ] }, { "cell_type": "code", "execution_count": null, "id": "alone-amateur", "metadata": {}, "outputs": [], "source": [ "filename = root+'backgroundstep.fits' # the output from the previous step\n", "\n", "step = Extract2dStep()\n", "step.save_results = True\n", "step.output_dir = output_dir\n", "result = step(filename)" ] }, { "cell_type": "markdown", "id": "centered-receiver", "metadata": {}, "source": [ "The bounds of the 2D extraction are printed in the above verbose output as \"Subarray x/y-extents\". They are determined based on the expected wavelength and slit projection limits appropriate for the aperture and grating/filter combination being used." ] }, { "cell_type": "code", "execution_count": null, "id": "configured-journal", "metadata": {}, "outputs": [], "source": [ "# The output file has the suffix _extract2dstep appended:\n", "e2d_output_file = root+'extract2dstep.fits'\n", "print(e2d_output_file)\n", "\n", "# load the output into a data model container\n", "e2d = datamodels.open(e2d_output_file)\n", "\n", "# plot the image\n", "show_image(e2d.slits[0].data, -5.e2, 25.e2, aspect=1., scale='linear', units='DN/s')" ] }, { "cell_type": "markdown", "id": "smaller-proportion", "metadata": {}, "source": [ "Now let's use the WCS object to find the source trace in sky coordinates. This involves some of the transformations that were set up in the assign_wcs step. Various WCS parameters can be calculated per pixel, including wavelength, sky coordinates, and internal instrument optical plane coordinates." ] }, { "cell_type": "code", "execution_count": null, "id": "appropriate-friendly", "metadata": {}, "outputs": [], "source": [ "# get the WCS object\n", "wcsobj = e2d.slits[0].meta.wcs\n", "\n", "# using coordinate transforms embedded in the WCS object,\n", "# we can calculate wavelength and sky coordinates for each pixel in the 2D spectrum\n", "y, x = np.mgrid[:e2d.slits[0].data.shape[0], : e2d.slits[0].data.shape[1]]\n", "det2sky = wcsobj.get_transform('detector', 'world') # this is the transform from detector space (pixel coordinates) to the sky (RA, Dec in degrees)\n", "calra, caldec, calwave = det2sky(x, y)\n", "calwaves = calwave[0, :] # get the set of wavelength values for the first pixel row, for later use\n", "\n", "# get the source RA and Dec coordinates from the metadata (also located in the header of the fits SCI extension)\n", "source_ra = e2d.slits[0].meta.target.ra\n", "source_dec = e2d.slits[0].meta.target.dec\n", "print('catalog RA,DEC:', source_ra, source_dec)\n", " \n", "# find the trace of the source coordinates\n", "calxm = np.arange(calra.shape[1])\n", "calym = np.arange(calra.shape[0])\n", "calym_ra = np.zeros(calra.shape[1])\n", "calym_dec = np.zeros(calra.shape[1])\n", "for j in range(len(calxm)):\n", " col = calra[:, j]\n", " calym_ra[j] = np.interp(source_ra, col[np.isfinite(col)], calym[np.isfinite(col)])\n", " col = caldec[:, j]\n", " calym_dec[j] = np.interp(source_dec, col[np.isfinite(col)], calym[np.isfinite(col)])\n", " \n", "# plot\n", "show_image(e2d.slits[0].data, -5.e2, 25.e2, aspect=1., scale='linear', units='DN/s')\n", "plt.plot(calxm, calym_ra, color='purple', alpha=0.9, linewidth=2)\n", "plt.plot(calxm, calym_dec, color='red', alpha=0.9, linewidth=2)" ] }, { "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 point. 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+'extract2dstep.fits' # the output from the previous step\n", "\n", "step = SourceTypeStep()\n", "step.save_results = True\n", "step.output_dir = output_dir\n", "result = step(filename)" ] }, { "cell_type": "markdown", "id": "little-ability", "metadata": {}, "source": [ "keyword SRCTYPE should be set to \"POINT\" in this case, because PPS database value SRCTYAPT is \"UNKNOWN\"" ] }, { "cell_type": "code", "execution_count": null, "id": "renewable-marking", "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": "refined-confusion", "metadata": {}, "source": [ "### wavecorr " ] }, { "cell_type": "markdown", "id": "straight-blind", "metadata": {}, "source": [ "The wavelength values per pixel need to be corrected when a point source is not located at the center of a slit aperture in the dispersion direction. (Due to a bug in how the source slit position is currently being calculated, the correction is effectively zero for FS data regardless of the dither position.)" ] }, { "cell_type": "code", "execution_count": null, "id": "sustained-spending", "metadata": {}, "outputs": [], "source": [ "filename = root+'sourcetypestep.fits' # the output from the previous step\n", "\n", "step = WavecorrStep()\n", "step.save_results = True\n", "step.output_dir = output_dir\n", "result = step(filename)" ] }, { "cell_type": "code", "execution_count": null, "id": "excess-block", "metadata": {}, "outputs": [], "source": [ "# The output file has the suffix _wavecorrstep appended:\n", "wavecor_output_file = root+'wavecorrstep.fits'\n", "\n", "# load the output into a data model container\n", "wavecor = datamodels.open(wavecor_output_file)" ] }, { "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+'wavecorrstep.fits' # the output from the previous step\n", "\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": "passive-healthcare", "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.slits[0].data, 0., 5., aspect=1., scale='Asinh', units='')" ] }, { "cell_type": "code", "execution_count": null, "id": "velvet-darwin", "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.slits[0].data, -5.e2, 25.e2, aspect=1., scale='linear', 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 geometric and diffraction losses of a non-centered point source incurred by the slit aperture. This is only relative to a centered source, as the absolute pathloss in that case has already been corrected by the F flat component of the flat field correction. (Due to the aforementioned bug in how the source slit position is currently being calculated, this correction is effectively one for FS data regardless of the dither position.)" ] }, { "cell_type": "code", "execution_count": null, "id": "noticed-helen", "metadata": {}, "outputs": [], "source": [ "filename = root+'flatfieldstep.fits' # the output from the previous step\n", "\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": "enhanced-hebrew", "metadata": {}, "outputs": [], "source": [ "# plot the corrected science exposure\n", "show_image(ploss.slits[0].data, -5.e2, 25.e2, aspect=1., scale='linear', 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", "step = PhotomStep()\n", "step.save_results = True\n", "step.output_dir = output_dir\n", "result = step(filename)" ] }, { "cell_type": "code", "execution_count": null, "id": "fantastic-edmonton", "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.slits[0].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.slits[0].meta.photometry.conversion_microjanskys)\n", "# this is only applicable for extended sources" ] }, { "cell_type": "markdown", "id": "accomplished-camping", "metadata": {}, "source": [ "### resample_spec " ] }, { "cell_type": "markdown", "id": "vietnamese-venice", "metadata": {}, "source": [ "The fully calibrated 2D spectrum is resampled onto a rectified grid 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 common grid during the level 3 processing." ] }, { "cell_type": "code", "execution_count": null, "id": "freelance-extra", "metadata": {}, "outputs": [], "source": [ "filename = root+'photomstep.fits' # the output from the previous step\n", "\n", "step = ResampleSpecStep()\n", "step.save_results = True\n", "step.output_dir = output_dir\n", "result = step(filename)" ] }, { "cell_type": "code", "execution_count": null, "id": "soviet-treat", "metadata": {}, "outputs": [], "source": [ "# The output file has the suffix _resamplespecstep appended:\n", "s2d_output_file = root+'resamplespecstep.fits'\n", "\n", "# load the output into a data model container\n", "s2d = datamodels.open(s2d_output_file)" ] }, { "cell_type": "code", "execution_count": null, "id": "opening-grove", "metadata": {}, "outputs": [], "source": [ "# plot the rectified 2D spectrum\n", "show_image(s2d.slits[0].data, -5.e2, 25.e2, aspect=1., scale='linear', units='MJy')" ] }, { "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 resampled 2D product, using an aperture size specified from a reference file. An aperture correction is also applied to account for flux outside of the extraction aperture. As with the resample_spec level 2b product, this is intended for quick-look purposes only. An independent extraction of the combined rectified 2D spectrum will be done during level 3 processing." ] }, { "cell_type": "code", "execution_count": null, "id": "earned-bunny", "metadata": {}, "outputs": [], "source": [ "filename = root+'resamplespecstep.fits' # the output from the previous step\n", "\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+'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": "static-slope", "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 nod set in this example, this pipeline will combine all of the exposures into a single output product. Both a rectified 2D and extracted 1D spectrum are generated. The pipeline includes 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, we will skip this step because the simulations are noiseless, which would result in spurious results. 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": "meaning-internship", "metadata": {}, "outputs": [], "source": [ "# run the calwebb_spec3 pipeline using association file of cal.fits files\n", "\n", "# get the association file\n", "asn = output_dir+\"spec3_all_asn.json\"\n", "\n", "# loop over the association files (one per exposure)\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 noise\n", "result = spec3(asn)" ] }, { "cell_type": "code", "execution_count": null, "id": "eastern-scale", "metadata": {}, "outputs": [], "source": [ "# display level 3 products\n", "# combined 2D rectified spectrum\n", "s2d3 = datamodels.open(output_dir+'nirspec_fssim_s00001_nirspec_clear-prism-s200a1-subs200a1_s2d.fits')\n", "s2d3sci = s2d3.data\n", "\n", "# get the wavelength scale\n", "wcsobj = s2d3.meta.wcs\n", "y, x = np.mgrid[:s2d3sci.shape[0], : s2d3sci.shape[1]]\n", "det2sky = wcsobj.get_transform('detector', 'world')\n", "s2d3ra, s2d3dec, s2d3wave = det2sky(x, y)\n", "s2d3waves = s2d3wave[0, :]\n", "xtint = np.arange(100, s2d3sci.shape[1], 100)\n", "xtlab = np.round(s2d3waves[xtint], 2)\n", "\n", "# show the combined s2d image\n", "show_image(s2d3sci, -5.e2, 2.5e3, aspect=5., scale='linear')\n", "plt.xticks(xtint, xtlab)\n", "plt.xlabel('wavelength (microns)')" ] }, { "cell_type": "code", "execution_count": null, "id": "related-great", "metadata": {}, "outputs": [], "source": [ "# show the overlap image (\"CON\" fits extension)\n", "s2d3con = s2d3.con\n", "show_image(s2d3con, 0, 512, aspect=5., scale='linear', units='overlap index')" ] }, { "cell_type": "markdown", "id": "different-novel", "metadata": {}, "source": [ "This indicates which input exposures (identified by a binary flag) contributed to the flux in each output pixel. (Note that the values at rows >~ 27 are incorrect as a result of known a bug in this pipeline version.) " ] }, { "cell_type": "code", "execution_count": null, "id": "responsible-consciousness", "metadata": {}, "outputs": [], "source": [ "# combined 1D extracted spectrum\n", "x1d3 = datamodels.open(output_dir+'nirspec_fssim_s00001_nirspec_clear-prism-s200a1-subs200a1_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=(12, 8))\n", "plt.plot(x1d3wave, x1d3flux)\n", "plt.xlabel('wavelength (microns)')\n", "plt.ylabel('flux (Jy)')\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "hindu-feature", "metadata": {}, "source": [ "### Rerun extract_1d with a different extraction aperture " ] }, { "cell_type": "code", "execution_count": null, "id": "executive-twelve", "metadata": {}, "outputs": [], "source": [ "# show the contents of the extract_1d reference file\n", "reffile = output_dir+'jwst_nirspec_extract1d_0004.json'\n", "with open(reffile) as f_obj:\n", " e1dref = json.load(f_obj)\n", "e1dref" ] }, { "cell_type": "code", "execution_count": null, "id": "arctic-summer", "metadata": {}, "outputs": [], "source": [ "# edit the extraction aperture width\n", "for params in e1dref['apertures']:\n", " if (params['id'] == 'S200A1'):\n", " params['extract_width'] = 3\n", "\n", "# create new reference file\n", "newfile = open(output_dir+'jwst_nirspec_extract1d_news200a1.json', 'w')\n", "json.dump(e1dref, newfile)\n", "newfile.close()" ] }, { "cell_type": "code", "execution_count": null, "id": "jewish-commons", "metadata": {}, "outputs": [], "source": [ "# rerun extract_1d step with the new reference file\n", "root = output_dir+'nirspec_fssim_s00001_nirspec_clear-prism-s200a1-subs200a1'\n", "filename = root+'_s2d.fits' # the level 3 2D spectrum\n", "\n", "step = Extract1dStep()\n", "step.override_extract1d = output_dir+'jwst_nirspec_extract1d_news200a1.json'\n", "step.save_results = True\n", "step.output_dir = output_dir\n", "result = step(filename)" ] }, { "cell_type": "code", "execution_count": null, "id": "accepted-awareness", "metadata": {}, "outputs": [], "source": [ "# The output file has the suffix _extract1dstep appended:\n", "x1d3_output_file = root+'_extract1dstep.fits'\n", "\n", "# load the output into a data model container\n", "x1d3new = datamodels.open(x1d3_output_file)" ] }, { "cell_type": "code", "execution_count": null, "id": "primary-static", "metadata": {}, "outputs": [], "source": [ "# get wavelength & flux\n", "x1d3newwave = x1d3new.spec[0].spec_table.WAVELENGTH\n", "x1d3newflux = x1d3new.spec[0].spec_table.FLUX\n", "\n", "# plot & compare with the original version\n", "fig = plt.figure(figsize=(12, 8))\n", "plt.plot(x1d3newwave, x1d3newflux, color='blue')\n", "plt.plot(x1d3wave, x1d3flux, color='gray')\n", "plt.xlabel('wavelength (microns)')\n", "plt.ylabel('flux (Jy)')\n", "plt.show()" ] }, { "cell_type": "markdown", "id": "aware-criterion", "metadata": {}, "source": [ "Using a narrower extraction aperture in this case results in lower flux because the default aperture correction reference file is only a placeholder. For real data obtained in flight, the two spectra should overlap, though of course the narrow-aperture version would have larger uncertainties." ] }, { "cell_type": "markdown", "id": "38d6a63c-6626-409f-b4aa-27f5be3d29a5", "metadata": {}, "source": [ "\"Space" ] } ], "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 }