{ "cells": [ { "cell_type": "markdown", "id": "4097f313", "metadata": {}, "source": [ "\"stsci_logo\" " ] }, { "cell_type": "markdown", "id": "56018571", "metadata": {}, "source": [ "\n", "# MIRI MRS Calibration Notebook #1 (Detailed Spec2/Spec3 Walkthrough) #" ] }, { "cell_type": "markdown", "id": "54a14522", "metadata": {}, "source": [ "**Author**: David Law, AURA Associate Astronomer, MIRI branch\n", "
\n", "**Last Updated**: November 29, 2021\n", "
\n", "**Pipeline Version**: 1.3.3 \n", "**Updated to**: 1.5.2 on June 10, 2022" ] }, { "cell_type": "markdown", "id": "c0ef1e21", "metadata": {}, "source": [ "## Table of contents\n", "1. [Introduction](#intro)
\n", " 1.1 [Purpose of this Notebook](#purpose)
\n", " 1.2 [Input Simulations](#inputs)
\n", " 1.3 [Caveats for Simulated Data](#mirisim)
\n", "2. [Setup](#setup)
\n", " 2.1 [CRDS Context](#crds)
\n", " 2.2 [Python Imports](#imports)
\n", " 2.3 [Data I/O Directories](#iodir)
\n", " 2.4 [Reprocessing Flag](#redo)
\n", "3. [Detector1 Pipeline](#det1)
\n", "4. [Spec2 Pipeline](#spec2)
\n", " 4.1 [Assign WCS](#wcs)
\n", " 4.2 [Pixel-wise Background Subtraction](#spec2bg)
\n", " 4.3 [Imprint Subtraction and MSA Failed-Open Flagging (NIRSpec-only)](#imprintmsa)
\n", " 4.4 [Flatfield Correction](#flat)
\n", " 4.5 [Source Type Identification](#stype)
\n", " 4.6 [Straylight Subtraction (MIRI-only)](#stray)
\n", " 4.7 [Fringe Flat (MIRI-only)](#fringe)
\n", " 4.8 [Path-Loss Correction (NIRSpec-only)](#pathloss)
\n", " 4.9 [Flux Calibration](#fluxcal)
\n", " 4.10 [Quicklook Cube Building](#spec2cube)
\n", " 4.11 [Quicklook Spectral Extraction](#spec2x1d)
\n", "5. [Spec3 Pipeline](#spec3)
\n", " 5.1 [Association Files](#l3assoc)
\n", " 5.2 [Moving Target WCS](#spec3wcs)
\n", " 5.3 [Master Background Subtraction](#masterbg)
\n", " 5.4 [Residual Background Matching (MIRI-only)](#residbg)
\n", " 5.5 [Outlier Detection](#orej)
\n", " 5.6 [Cube Creation](#spec3cube)
\n", " 5.7 [Spectral Extraction](#spec3x1d)
" ] }, { "cell_type": "markdown", "id": "77fbc6d7", "metadata": {}, "source": [ "
" ] }, { "cell_type": "markdown", "id": "dea72d3b", "metadata": {}, "source": [ "1.-Introduction \n", "------------------" ] }, { "cell_type": "markdown", "id": "9ac76d41", "metadata": {}, "source": [ "### 1.1-Purpose of this Notebook ###\n", "\n", "In this notebook we provide an introduction to JWST spectroscopic mode data reduction, focusing on the specific use case of the integral-field unit (IFU) observing mode. As a demonstration case we use IFU data for the MIRI Medium Resolution Spectrometer (MRS) IFU; while many of the pipeline steps are similar for both MIRI and NIRSpec IFUs we note differences where they arise.\n", "\n", "We will start with a simple simulated MRS point source observation (created using mirisim: https://wiki.miricle.org/Public/MIRISim_Public), process the data through the Detector1 pipeline (which turns raw detector counts into uncalibrated rate images), the Spec2 pipeline (which turns uncalibrated rate images into calibrated rate images), and the Spec3 pipeline (which turns calibrated rate images into composite data cubes and extracted 1d spectra).\n", "\n", "Since the Detector1 pipeline has been discussed extensively in previous JWebbinars we will not dig into that stage of the pipeline in detail, and focus instead on the Spec2 and Spec3 stages. We will step individually through each step in these two pipeline stages, discuss how they work, and examine some sample outputs. Since we'll be examining each step individually, this notebook is thus not a good template to use for designing your own notebook to process/inspect large quantities of observational data. This use case is more directly addressed by MRS Calibration Notebooks #2 and #3 (point sources and extended sources respectively).\n", "\n", "A few additional caveats:\n", "- This notebook covers the v1.3.2 baseline pipeline as it existed in October 2021. The pipeline is under continuous development and there are therefore some changes in the latest pipeline build that will not be reflected here.\n", "- Likewise, there are some advanced algorithms slated for development prior to cycle 1 observations that will not be discussed here." ] }, { "cell_type": "markdown", "id": "4bc6e500", "metadata": {}, "source": [ "### 1.2-Input Simulations ###\n", "\n", "As input to this notebook, we'll be using a 4-pt dithered observation of a point source created using mirisim that covers the 1C and 2C bands (i.e., just one of the two MIRI detectors, with the LONG grating configuration). The point source spectrum is chosen to be astrophysically realistic, in this case taken from the galaxy NGC 5728.\n", "\n", "\"nb1input\" " ] }, { "cell_type": "markdown", "id": "d0534256", "metadata": {}, "source": [ "### 1.3-Caveats for Simulated Data ###\n", "\n", "As noted above, in this notebook we will be processing simulated data created with the 'mirisim' tool. Like the pipeline, mirisim is also an evolving piece of software and there are multiple known issues that can cause problems. A few of the most important such mirisim issues include:\n", "\n", "- FAST mode has incorrect noise properties, rendering FAST mode data processed by the pipeline unreliable. The simulations in this notebook therefore use simulated SLOW mode data.\n", "\n", "- Extended sources are not simulated unless they meet a minimum size that varies with each band.\n", "\n", "- Point sources are not simulated properly, with the PSF profile being simulated incorrectly in the cores. Mirisim simulations therefore should not be used to study the PSF shape.\n", "\n", "- Reference pixels are not treated consistently, the refpix step of detector1 must therefore be turned off to process mirisim data without artifacts.\n", "\n", "- Channel 4 flux calibration is incorrect in mirisim. No workaround is currently available- channel 4 fluxes provided by the pipeline from simulated data will be incorrect.\n", "\n", "- WCS alignment is incorrect in mirisim, causing sources to jump in location by a couple of pixels between channels. No workaround is available- do not use mirisim data to test spatial alignment.\n", "\n", "- Flux conservation is not perfect within mirisim. Likewise, the aperture correction factors in use by the pipeline correspond to the expected performance in flight (to be udpated during on-orbit commissioning) and are not well matched to mirisim data. No workaround available, do not use mirisim data to test flux conservation.\n", "\n", "- mirisim does not add all of the necessary header keywords for the pipeline to know how to do background subtraction, identify source type, etc. In order to get these APT-derived keywords correct they will need to be set manually." ] }, { "cell_type": "markdown", "id": "274b4d23", "metadata": {}, "source": [ "
" ] }, { "cell_type": "markdown", "id": "eeb4f687", "metadata": {}, "source": [ "2.-Setup \n", "------------------\n", "\n", "In this section we set things up a number of necessary things in order for the pipeline to run successfully.\n", "\n", "First we'll import the various python packages that we're actually going to use in this notebook, including both generic utility functions and the actual pipeline modules themselves.\n", "\n", "Next, we'll specify the data directory structure that we want to use. In order to keep our filesystem clean we'll separate simulated inputs and outputs from each pipeline stage into their own folders.\n", "\n", "Finally, for convenience in this JWebbinar we'll define a flag that sets whether or not to actually run some of the longer pipeline steps in this notebook or just to rely upon cached reductions provided ahead of time. This is because some steps can take quite a long time to run, and in a short Webbinar we don't want to be waiting for them to all run in real time. This flag is set to False by default for use in the live Webbinar; if you want to experiment with running all steps yourself ahead of time just set this flag to True. Total notebook runtime with True/False values is about 15 minutes vs 6 minutes.\n", "\n", "Note that this runtime reflects substantial gains in the speed of 3d cube building in the build 1.3.2 JWST pipeline, which now uses C routines to perform cube building instead of pure python. Previously, total notebook runtime was approximately 1.5 hours." ] }, { "cell_type": "markdown", "id": "fe2d4497", "metadata": {}, "source": [ "### 2.1-CRDS Context ###" ] }, { "cell_type": "markdown", "id": "3b7677f0", "metadata": {}, "source": [ "If we wanted to set a specific CRDS reference file context (see https://jwst-crds.stsci.edu/) we would do it here. However, it is not currently necessary for version 1.3.2 of the pipeline." ] }, { "cell_type": "code", "execution_count": null, "id": "4c6722ad", "metadata": {}, "outputs": [], "source": [ "# Comment out this line if you want to use the latest reference files tagged for a specific pipeline version\n", "#%env CRDS_CONTEXT jwst_0771.pmap" ] }, { "cell_type": "markdown", "id": "18258a5b", "metadata": {}, "source": [ "### 2.2-Python Imports ###" ] }, { "cell_type": "code", "execution_count": null, "id": "7b632113", "metadata": { "scrolled": false }, "outputs": [], "source": [ "# First let's use the entire available screen width for the notebook\n", "from IPython.core.display import display, HTML\n", "display(HTML(\"\"))" ] }, { "cell_type": "code", "execution_count": null, "id": "285c5c06", "metadata": {}, "outputs": [], "source": [ "# Basic system utilities for interacting with files\n", "import glob\n", "import sys\n", "import os\n", "# Update to a path in your system\n", "os.environ[\"CRDS_PATH\"] = \"/path/to/my/folder/\"\n", "os.environ[\"CRDS_SERVER_URL\"] = \"https://jwst-crds-pub.stsci.edu\"\n", "import time\n", "import shutil\n", "import warnings\n", "import zipfile\n", "import urllib.request\n", "\n", "# Astropy utilities for opening FITS and ASCII files\n", "from astropy.io import fits\n", "from astropy.io import ascii\n", "from astropy.utils.data import download_file\n", "# Astropy utilities for making plots\n", "from astropy.visualization import LinearStretch, LogStretch, ImageNormalize, ZScaleInterval\n", "\n", "# Numpy for doing calculations\n", "import numpy as np\n", "\n", "# Matplotlib for making plots\n", "import matplotlib.pyplot as plt\n", "from matplotlib import rc" ] }, { "cell_type": "code", "execution_count": null, "id": "a342207b", "metadata": {}, "outputs": [], "source": [ "# Import the base JWST package and warn if not the expected version\n", "import jwst\n", "\n", "print('JWST pipeline version', jwst.__version__)" ] }, { "cell_type": "code", "execution_count": null, "id": "ec3e18e7", "metadata": { "scrolled": false }, "outputs": [], "source": [ "# JWST pipelines (encompassing many steps)\n", "from jwst.pipeline import Detector1Pipeline\n", "from jwst.pipeline import Spec2Pipeline\n", "from jwst.pipeline import Spec3Pipeline\n", "\n", "# Individual JWST pipeline steps\n", "from jwst.assign_wcs import AssignWcsStep\n", "from jwst.background import BackgroundStep\n", "from jwst.flatfield import FlatFieldStep\n", "from jwst.srctype import SourceTypeStep\n", "from jwst.straylight import StraylightStep\n", "from jwst.fringe import FringeStep\n", "from jwst.photom import PhotomStep\n", "from jwst.cube_build import CubeBuildStep\n", "from jwst.extract_1d import Extract1dStep\n", "from jwst.cube_skymatch import CubeSkyMatchStep\n", "from jwst.master_background import MasterBackgroundStep\n", "from jwst.outlier_detection import OutlierDetectionStep\n", "\n", "# JWST pipeline utilities\n", "from jwst import datamodels # JWST datamodels\n", "from jwst.associations import asn_from_list as afl # Tools for creating association files\n", "from jwst.associations.lib.rules_level2_base import DMSLevel2bBase # Definition of a Lvl2 association file\n", "from jwst.associations.lib.rules_level3_base import DMS_Level3_Base # Definition of a Lvl3 association file\n", "\n", "from stcal import dqflags # Utilities for working with the data quality (DQ) arrays" ] }, { "cell_type": "markdown", "id": "dcecf94b", "metadata": {}, "source": [ "### 2.3-Data I/O Directories ###" ] }, { "cell_type": "markdown", "id": "51e8a8a9", "metadata": {}, "source": [ "Data for this notebook comes in two parts, a stage0/ directory that contains the mirisim inputs to the pipeline, and a cache/ directory containing pre-reduced stage1, stage2, and stage3 pipeline results. This cache directory is intended to provide reference examples against which new reductions can be compared, or to serve as the source of intermediate-stage pipeline inputs if skipping lengthy processing steps in this notebook.\n", "\n", "These directories are zipped together and can be downloaded from Box automatically using the code below. However, since the zip file is quite large (~ 1 GB) it can also be downloaded through any other means and extracted into the target data directory." ] }, { "cell_type": "code", "execution_count": null, "id": "5a9b164a", "metadata": {}, "outputs": [], "source": [ "# input_dir should point to where you installed the data from Box\n", "input_dir = './MRS_data1/'\n", "\n", "ziplink = 'https://stsci.box.com/shared/static/gcphs4m12j7l8zmqj4cgbmw0be12dzz3.zip'\n", "zipfilename = 'mrs_data1.zip'\n", "\n", "# Look for a stage0/ directory, if not found then download the data\n", "if os.path.isdir(os.path.join(input_dir, 'stage0/')):\n", " print('stage0/ input directory already exists, skipping download.')\n", "else:\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(input_dir, zipfilename))\n", " zf = zipfile.ZipFile(os.path.join(input_dir, zipfilename), 'r')\n", " print('Extracting zip file...')\n", " zf.extractall(input_dir)\n", " print('Cleaning up.')\n", " os.remove(os.path.join(input_dir, zipfilename))" ] }, { "cell_type": "code", "execution_count": null, "id": "51b43a48", "metadata": {}, "outputs": [], "source": [ "# Set up mirisim input path\n", "mirisim_dir = os.path.join(input_dir, 'stage0/') # Simulated inputs are here\n", "\n", "# Set up cache paths\n", "det1_cachedir = os.path.join(input_dir, 'cache/stage1/') # Cached Detector1 pipeline outputs are here\n", "spec2_cachedir = os.path.join(input_dir, 'cache/stage2/') # Cached Spec2 pipeline outputs are here\n", "spec3_cachedir = os.path.join(input_dir, 'cache/stage3/') # Cached Spec3 pipeline outputs are here" ] }, { "cell_type": "code", "execution_count": null, "id": "6d2cb178", "metadata": {}, "outputs": [], "source": [ "# output_dir should point to where you want the output results to go\n", "output_dir = './MRS_output1/'\n", "\n", "# Specify output directories to keep data products organized\n", "det1_dir = os.path.join(output_dir, 'stage1/') # Detector1 pipeline outputs will go here\n", "spec2_dir = os.path.join(output_dir, 'stage2/') # Spec2 pipeline outputs will go here\n", "spec3_dir = os.path.join(output_dir, 'stage3/') # Spec3 pipeline outputs will go here\n", "\n", "# We need to check that the desired output directories exist, and if not create them\n", "if not os.path.exists(det1_dir):\n", " os.makedirs(det1_dir)\n", "if not os.path.exists(spec2_dir):\n", " os.makedirs(spec2_dir)\n", "if not os.path.exists(spec3_dir):\n", " os.makedirs(spec3_dir)" ] }, { "cell_type": "markdown", "id": "0ad4c945", "metadata": {}, "source": [ "### 2.4-Reprocessing Flag ###" ] }, { "cell_type": "markdown", "id": "e618da8a", "metadata": {}, "source": [ "Since some parts of the pipeline take a long time to run, for a first use of this notebook we will disable those steps and simply copy results out of the cache for informational purposes. In order to run the full pipeline on the data, or to run this notebook on your own simulated data not downloaded from the Box link above, this reprocessing flag will need to be enabled." ] }, { "cell_type": "code", "execution_count": null, "id": "a9140935", "metadata": {}, "outputs": [], "source": [ "# To rerun all steps use:\n", "# redolong = True\n", "# To skip lengthy steps and copy results from the cache use:\n", "redolong = False\n", "\n", "# If skipping lengthy steps, check that the cache has been installed properly\n", "if (redolong is False):\n", " if not os.path.exists(det1_cachedir):\n", " print('Could not find cache, running all steps.')\n", " redolong = True" ] }, { "cell_type": "markdown", "id": "74103475", "metadata": {}, "source": [ "
" ] }, { "cell_type": "markdown", "id": "1353b88a", "metadata": {}, "source": [ "3.-Detector1 Pipeline \n", "------------------\n", "
\n", "\n", "In this section we process our simulated data through the Detector1 pipeline to create Lvl2a data products (i.e., uncalibrated slope images). We won't go into detail however, as these steps are treated in detail elsewhere.\n", "\n", "See https://jwst-pipeline.readthedocs.io/en/latest/jwst/pipeline/calwebb_detector1.html\n", " \n", "
" ] }, { "cell_type": "code", "execution_count": null, "id": "658a880f", "metadata": {}, "outputs": [], "source": [ "# Start a timer to keep track of runtime\n", "time0 = time.perf_counter()" ] }, { "cell_type": "code", "execution_count": null, "id": "b4f854ba", "metadata": { "scrolled": false }, "outputs": [], "source": [ "# First we'll define a function that will call the detector1 pipeline with our desired set of parameters:\n", "def rundet1(filenames):\n", " det1 = Detector1Pipeline() # Instantiate the pipeline\n", " det1.output_dir = det1_dir # Specify where the output should go\n", " det1.refpix.skip = True # Skip the reference pixel subtraction (as it doesn't interact well with simulated data)\n", " det1.save_results = True # Save the final resulting _rate.fits files\n", " det1(filenames) # Run the pipeline on an input list of files" ] }, { "cell_type": "code", "execution_count": null, "id": "87f5c381", "metadata": { "scrolled": false }, "outputs": [], "source": [ "# Now let's look for input files in our (cached) mirisim simulation directory\n", "sstring = mirisim_dir + 'det*exp1.fits'\n", "simfiles = sorted(glob.glob(sstring))\n", "print('Found ' + str(len(simfiles)) + ' input files to process')" ] }, { "cell_type": "code", "execution_count": null, "id": "a41f942f", "metadata": { "scrolled": true }, "outputs": [], "source": [ "# Run the pipeline on these input files by a simple loop over our pipeline function\n", "\n", "# If rerunning long pipeline steps, actually run the step\n", "if (redolong is True):\n", " for file in simfiles:\n", " rundet1(file)\n", "# Otherwise, just copy cached outputs into our output directory structure\n", "else:\n", " sstring = os.path.join(det1_cachedir, 'det*rate.fits')\n", " files = sorted(glob.glob(sstring))\n", " for file in files:\n", " outfile = str.replace(file, det1_cachedir, det1_dir)\n", " shutil.copy(file, outfile)" ] }, { "cell_type": "markdown", "id": "6e5f71ba", "metadata": {}, "source": [ "Let's take a look at the output data products (i.e., uncalibrated slope data) to get an idea what they look like." ] }, { "cell_type": "code", "execution_count": null, "id": "f3eaffa8", "metadata": {}, "outputs": [], "source": [ "# Look for our _rate.fits files produced by the Detector1 pipeline\n", "sstring = det1_dir + 'det*rate.fits'\n", "ratefiles = sorted(glob.glob(sstring))\n", "# And print them out so that we can see them\n", "ratefiles" ] }, { "cell_type": "code", "execution_count": null, "id": "6b9b4e26", "metadata": {}, "outputs": [], "source": [ "# We'll open the first image in the list and take a look at its contents\n", "hdu1 = fits.open(ratefiles[0])\n", "hdu1.info()" ] }, { "cell_type": "code", "execution_count": null, "id": "ec1a12e6", "metadata": {}, "outputs": [], "source": [ "# Let's show what the SCI extension of the first two files looks like\n", "image1 = hdu1['SCI'].data\n", "hdu2 = fits.open(ratefiles[1])\n", "image2 = hdu2['SCI'].data\n", "\n", "# Use a classic ZScale normalization\n", "norm = ImageNormalize(image1, interval=ZScaleInterval(), stretch=LinearStretch())\n", "\n", "rc('axes', linewidth=2) \n", "fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(7, 7), dpi=100)\n", "\n", "# And plot the data\n", "ax1.imshow(image1, cmap='gray', norm=norm, origin='lower')\n", "ax1.set_title('Exposure 1')\n", "ax1.set_xlabel('X pixel')\n", "ax1.set_ylabel('Y pixel')\n", "\n", "ax2.imshow(image2, cmap='gray', norm=norm, origin='lower')\n", "ax2.set_title('Exposure 2')\n", "ax2.set_xlabel('X pixel')" ] }, { "cell_type": "markdown", "id": "36b4a626", "metadata": {}, "source": [ "Figure 1: Each exposure above shows the slope image for a MIRI MRS exposure. X pixels 1-500 (roughly) correspond to data from Channel 1, and pixels 500-1024 correspond to data from Channel 2. Data are dispersed along the spectral Y axis in each of multiple slices. In this example of a bright point source, the point source is visible in many different slices- note how different slices are illuminated in Exposure 1 vs Exposure 2 though as we have dithered the source location between the two." ] }, { "cell_type": "code", "execution_count": null, "id": "8d5b8f37", "metadata": {}, "outputs": [], "source": [ "# Close our files behind us\n", "hdu1.close()\n", "hdu2.close()" ] }, { "cell_type": "code", "execution_count": null, "id": "7dbc97a4", "metadata": {}, "outputs": [], "source": [ "# Print out the time benchmark\n", "time1 = time.perf_counter()\n", "print(f\"Runtime so far: {time1 - time0:0.4f} seconds\")" ] }, { "cell_type": "markdown", "id": "11f21a4f", "metadata": {}, "source": [ "We're now done with the Detector1 pipeline!" ] }, { "cell_type": "markdown", "id": "edfb7d82", "metadata": {}, "source": [ "
" ] }, { "cell_type": "markdown", "id": "d20f10d1", "metadata": {}, "source": [ "4.-Spec2 Pipeline \n", "------------------\n", "\n", "
\n", "In this section we process our simulated data through the Spec2 pipeline in order to produce Lvl2b data products (i.e., calibrated slope images and quick-look data cubes). Since we're going to run each step of this pipeline individually we won't set up a master spec2 pipeline function- see Notebook #2 for a guide to doing that.\n", "\n", "Note that the individual steps ordinarily wouldn't produce output files if we were running Spec2 as a pipeline; they would just feed from one step directly into the next. However, we'll break them out individual here so that we can inspect them.\n", "\n", "See https://jwst-pipeline.readthedocs.io/en/latest/jwst/pipeline/calwebb_spec2.html\n", "
" ] }, { "cell_type": "markdown", "id": "4fec624c", "metadata": {}, "source": [ "### 4.1-Assign WCS ###\n", "\n", "
\n", "This first step is critical for many later stages of the pipeline; this is where the distortion solution of the instrument is put together with the various telescope pointing keywords to determine the overall world coordinate solution mapping detector pixels to their RA, DEC, WAVELENGTH on the sky. This transformation is extremely complex for the MRS, which has many slices for each of two channels that both map to overlapping regions of the sky simultaneously at different wavelengths. This information will be stored in an ASDF extension to the FITS file which can be queried using the JWST datamodels code, or both other ASDF parsers.\n", "\n", "See https://jwst-pipeline.readthedocs.io/en/latest/jwst/assign_wcs/index.html\n", "
" ] }, { "cell_type": "code", "execution_count": null, "id": "d117db09", "metadata": { "scrolled": true }, "outputs": [], "source": [ "# Call the step, specifying that we want results saved into the spec2_dir directory\n", "for file in ratefiles:\n", " AssignWcsStep.call(file, save_results=True, output_dir=spec2_dir)" ] }, { "cell_type": "code", "execution_count": null, "id": "88baaa9f", "metadata": { "scrolled": false }, "outputs": [], "source": [ "# Look for the _assignwcsstep.fits files produced by this step\n", "sstring = spec2_dir + 'det*assignwcsstep.fits'\n", "wcsfiles = sorted(glob.glob(sstring))\n", "# And print them out so that we can see them\n", "wcsfiles" ] }, { "cell_type": "markdown", "id": "13915742", "metadata": {}, "source": [ "The assign_wcs step doesn't modify the science data, so this will look identical to the \\_rate.fits files above. However, now the ASDF extension contains all of the information about the distortion transforms between detector pixel values and world coordinates. A detailed treatment of how to interact with this WCS is beyond the scope of this notebook, but we can show an example quickly:" ] }, { "cell_type": "code", "execution_count": null, "id": "ae459278", "metadata": { "scrolled": false }, "outputs": [], "source": [ "# Open the first file using the JWST datamodels\n", "image = datamodels.ImageModel(wcsfiles[0])\n", "\n", "# And show the available transforms\n", "print(image.meta.wcs)" ] }, { "cell_type": "code", "execution_count": null, "id": "21fc1911", "metadata": { "scrolled": false }, "outputs": [], "source": [ "# Here's an example transformation giving the world coordinates of pixel x=376,y=589\n", "ra, dec, wave = image.meta.wcs.transform(\"detector\", \"world\", 376, 589)\n", "print('RA = ', ra, ' deg')\n", "print('DEC = ', dec, ' deg')\n", "print('Wavelength = ', wave, ' micron')" ] }, { "cell_type": "markdown", "id": "945564d2", "metadata": {}, "source": [ "### 4.2-Pixel-wise background subtraction ###\n", "\n", "
\n", "The MIRI MRS is expected to see a substantial background signal, especially in Channel 4 where the telescope thermal background becomes significant. It will therefore be necessary to subtract this background (strictly, a foreground) signal from the data prior to analysis. Ideally, this background should be sufficiently smooth and uniform that it is possible to model it and subtract the model from the science data (see Spec3- Master Background). However, it is possible that the signal has sufficient spatial/spectral variability that it may prove necessary to do background subtraction via direct subtraction of on/off pointings.\n", "\n", "This background step in the Spec2 pipeline is largely a placeholder for if it becomes necessary to do such a direct subtraction at the cost of increased total noise. As such, we do not recommend running this step at the present time- additional guidance will be provided after on-orbit commissioning. If we *were* to run the step, here's one way in which in could be done.\n", "\n", "\n", "See https://jwst-pipeline.readthedocs.io/en/latest/jwst/background_step/\n", "
" ] }, { "cell_type": "code", "execution_count": null, "id": "330d62e1", "metadata": {}, "outputs": [], "source": [ "# Look for our assignwcsstep.fits files produced by the assign_wcs step\n", "sstring = spec2_dir + 'det*assignwcsstep.fits'\n", "wcsfiles = sorted(glob.glob(sstring))\n", "# And print them out so that we can see them\n", "wcsfiles" ] }, { "cell_type": "markdown", "id": "d1d8a6be", "metadata": {}, "source": [ "Define which exposures we're going to call background pointings for each science exposure. If we were observing an extended source with a dedicated background this would be easy- the background pointings are any dedicated background observations. Since we're simulating dithered point source observations however this is a little trickier. With a four-point dither pattern we have A-B-A'-B' pointings; where the A-B offset distance is large and the A' and B' pointings are close to the A and B pointings respectively. We'll use A/B and A'/B' as background pairs." ] }, { "cell_type": "code", "execution_count": null, "id": "9c99a50d", "metadata": {}, "outputs": [], "source": [ "# Alternate background assignments\n", "bgfiles = wcsfiles.copy()\n", "bgfiles[0] = wcsfiles[1]\n", "bgfiles[1] = wcsfiles[0]\n", "bgfiles[2] = wcsfiles[3]\n", "bgfiles[3] = wcsfiles[2]" ] }, { "cell_type": "code", "execution_count": null, "id": "0f5233c4", "metadata": { "scrolled": true }, "outputs": [], "source": [ "# Loop over the input files, specifying background to subtract from each\n", "for ii in range(0, 4):\n", " BackgroundStep.call(wcsfiles[ii], [bgfiles[ii]], save_results=True, output_dir=spec2_dir)" ] }, { "cell_type": "code", "execution_count": null, "id": "772ad768", "metadata": {}, "outputs": [], "source": [ "# Look for our backgroundstep.fits files produced by the background step\n", "sstring = spec2_dir + 'det*backgroundstep.fits'\n", "\n", "bgfiles = sorted(glob.glob(sstring))\n", "# And print them out so that we can see them\n", "bgfiles" ] }, { "cell_type": "code", "execution_count": null, "id": "8c5b67dd", "metadata": {}, "outputs": [], "source": [ "# Let's take a quick look at what it did by reading in the original \n", "# exposures 1 and 2, and the modified exposure 1\n", "hdu1 = fits.open(wcsfiles[0])\n", "image1 = hdu1['SCI'].data\n", "hdu2 = fits.open(wcsfiles[1])\n", "image2 = hdu2['SCI'].data\n", "hdu3 = fits.open(bgfiles[0])\n", "image3 = hdu3['SCI'].data\n", "\n", "# Use a classic ZScale normalization\n", "norm = ImageNormalize(image3, interval=ZScaleInterval(), stretch=LinearStretch())\n", "\n", "rc('axes', linewidth=2) \n", "fig, (ax1, ax2, ax3) = plt.subplots(1, 3, figsize=(9, 7), dpi=100)\n", "\n", "# And plot the data\n", "ax1.imshow(image1, cmap='gray', origin='lower', norm=norm)\n", "ax1.set_title('Frame 1')\n", "ax1.set_xlabel('X pixel')\n", "ax1.set_ylabel('Y pixel')\n", "\n", "ax2.imshow(image2, cmap='gray', origin='lower', norm=norm)\n", "ax2.set_title('Frame 2')\n", "ax2.set_xlabel('X pixel')\n", "\n", "ax3.imshow(image3, cmap='gray', origin='lower', norm=norm)\n", "ax3.set_title('Frame 1 - Frame 2')\n", "ax3.set_xlabel('X pixel')" ] }, { "cell_type": "markdown", "id": "40db2f52", "metadata": {}, "source": [ "Figure 2: Data before and after Spec2 pixel-wise background subtraction. Frame 1 and Frame 2 show two dithered frames of point source data, the right-hand column shows the difference between them computed by this step when Frame 2 is used as the background for Frame 1." ] }, { "cell_type": "code", "execution_count": null, "id": "8af236da", "metadata": {}, "outputs": [], "source": [ "# Close our files behind us\n", "hdu1.close()\n", "hdu2.close()\n", "hdu3.close()" ] }, { "cell_type": "markdown", "id": "bf7b36ef", "metadata": {}, "source": [ "### 4.3-Imprint Subtraction and MSA Failed-Open Flagging (NIRSpec-only) ###\n", "\n", "
\n", "These NIRSpec-specific steps account for the contamination of the detector by artifacts introduced by the microshutter array. The imprint subtraction step subtracts a dedicated imprint exposure to remove patterns created by the microshutter array. Likewise, the MSA failed-open flagging step identifies and flags detector pixels that are affected by stuck-open shutters in the microshutter array based on a master reference list.\n", "\n", "We do not treat these steps in detail in the present notebook.\n", "\n", "See https://jwst-pipeline.readthedocs.io/en/latest/jwst/imprint/index.html and https://jwst-pipeline.readthedocs.io/en/latest/jwst/msaflagopen/index.html\n", "
" ] }, { "cell_type": "markdown", "id": "d3c39d55", "metadata": {}, "source": [ "### 4.4-Flatfield correction ###\n", "\n", "
\n", "This step divides the array data by the pixel flatfield reference file. Ordinarily, this would therefore normalize out both differences in the pixel-to-pixel response function and differences in the pixel solid angle. However, for MIRI MRS the reference flatfield is currently unity everywhere that light falls on the detector, so this step has no practical effect. These corrections are instead taken care of in the PHOTOM step.\n", "\n", "See https://jwst-pipeline.readthedocs.io/en/latest/jwst/flatfield/index.html\n", "
" ] }, { "cell_type": "code", "execution_count": null, "id": "f9247b62", "metadata": {}, "outputs": [], "source": [ "# Look for our assignwcsstep.fits files produced by the assign_wcs step\n", "# (since we're going to skip background subtraction)\n", "sstring = spec2_dir + 'det*assignwcsstep.fits'\n", "\n", "wcsfiles = sorted(glob.glob(sstring))\n", "# And print them out so that we can see them\n", "wcsfiles" ] }, { "cell_type": "code", "execution_count": null, "id": "45dc83e3", "metadata": { "scrolled": true }, "outputs": [], "source": [ "# Call the step, specifying that we want results saved into the spec2_dir directory\n", "for file in wcsfiles:\n", " FlatFieldStep.call(file, save_results=True, output_dir=spec2_dir)" ] }, { "cell_type": "markdown", "id": "aefcefac", "metadata": {}, "source": [ "### 4.5-Source Type Identification ###\n", "\n", "
\n", "This step attempts to determine whether the observation should be considered as a POINT or EXTENDED source, and set the SRCTYPE keyword in the science extension header accordingly. There are no user-facing arguments that modify this step.\n", "\n", "Each instrument uses a slightly different logic to determine which kind of source to use. At present the logic is roughly as follows:\n", "- If a source type was specified in APT (keyword SRCTYAPT), use that.\n", "- If the dither pattern in use is a point-source specific pattern, assume this is a POINT source.\n", "- If no other guidance is provided, assume the target is an EXTENDED source (for MIRI MRS) or a POINT source (for NIRSpec IFU)\n", "\n", "Note that these are only relevant in that they set the default method by which the pipeline automatically extracts 1d spectra from the data cubes.\n", "\n", "See https://jwst-pipeline.readthedocs.io/en/latest/jwst/srctype/index.html\n", "
" ] }, { "cell_type": "code", "execution_count": null, "id": "2b058f32", "metadata": {}, "outputs": [], "source": [ "# Look for our assignwcsstep.fits files produced by the assign_wcs step\n", "sstring = spec2_dir + 'det*flatfieldstep.fits'\n", "flatfiles = sorted(glob.glob(sstring))\n", "# And print them out so that we can see them\n", "flatfiles" ] }, { "cell_type": "code", "execution_count": null, "id": "0ff50aed", "metadata": { "scrolled": true }, "outputs": [], "source": [ "# Call the step, specifying that we want results saved into the spec2_dir directory\n", "for file in flatfiles:\n", " SourceTypeStep.call(file, save_results=True, output_dir=spec2_dir)" ] }, { "cell_type": "code", "execution_count": null, "id": "010d03e1", "metadata": {}, "outputs": [], "source": [ "# Look for our sourcetypestep.fits files produced by the source type step\n", "sstring = spec2_dir + 'det*sourcetypestep.fits'\n", "srcfiles = sorted(glob.glob(sstring))\n", "# And print them out so that we can see them\n", "srcfiles" ] }, { "cell_type": "code", "execution_count": null, "id": "be32c911", "metadata": {}, "outputs": [], "source": [ "# We can look at the keyword in the first exposure\n", "hdu = fits.open(srcfiles[0])\n", "print('SRCTYPE = ', hdu['SCI'].header['SRCTYPE'])\n", "hdu.close()" ] }, { "cell_type": "markdown", "id": "c6f01042", "metadata": {}, "source": [ "In this case it was set to EXTENDED. This can't be kludged in the pipeline step, but we'll do so directly in the files in order to tell them that we're working with a POINT source:" ] }, { "cell_type": "code", "execution_count": null, "id": "54eb8e3a", "metadata": {}, "outputs": [], "source": [ "# Loop over files kludging the source type to POINT\n", "for file in srcfiles:\n", " hdu = fits.open(file)\n", " hdu['SCI'].header['SRCTYPE'] = 'POINT'\n", " hdu.writeto(file, overwrite=True)\n", " hdu.close()" ] }, { "cell_type": "code", "execution_count": null, "id": "beb5e8f8", "metadata": {}, "outputs": [], "source": [ "# Now the source type is set to POINT in the headers!\n", "hdu = fits.open(srcfiles[0])\n", "print('SRCTYPE = ', hdu['SCI'].header['SRCTYPE'])\n", "hdu.close()" ] }, { "cell_type": "markdown", "id": "115aeef6", "metadata": {}, "source": [ "### 4.6-Straylight Subtraction (MIRI-only) ###\n", "\n", "
\n", "The MIRI MRS has been observed to have appreciable straylight at short wavelengths in ground-test data, and this step is therefore designed to model and subtract this component from the detector data using the small interstitial regions of detector pixels between the illuminated slices. However, a limitation of the test data set that we're using is that mirisim does *not* add straylight to the data; as such, running this step on simulated data can only introduce artifacts rather than remove them.\n", "\n", "As such, we usually recommend skipping this step when working with simulated data. However, we'll run it here just to see what it looks like.\n", "\n", "See https://jwst-pipeline.readthedocs.io/en/latest/jwst/straylight/index.html\n", "
" ] }, { "cell_type": "code", "execution_count": null, "id": "2993bd4a", "metadata": { "scrolled": true }, "outputs": [], "source": [ "# Call the straylight step, specifying that we want results saved into the spec2_dir directory\n", "for file in srcfiles:\n", " StraylightStep.call(file, save_results=True, output_dir=spec2_dir)" ] }, { "cell_type": "code", "execution_count": null, "id": "71745c17", "metadata": {}, "outputs": [], "source": [ "# Look for our straylightstep.fits files produced by the straylight step\n", "sstring = spec2_dir + 'det*straylightstep.fits'\n", "strayfiles = sorted(glob.glob(sstring))\n", "# And print them out so that we can see them\n", "strayfiles" ] }, { "cell_type": "code", "execution_count": null, "id": "02612854", "metadata": {}, "outputs": [], "source": [ "# Let's show what the SCI extension of the first file before/after straylight subtraction looks like\n", "hdu1 = fits.open(srcfiles[0])\n", "image1 = hdu1['SCI'].data\n", "hdu2 = fits.open(strayfiles[0])\n", "image2 = hdu2['SCI'].data\n", "\n", "# Use a classic ZScale normalization\n", "norm = ImageNormalize(image1, interval=ZScaleInterval(), stretch=LinearStretch())\n", "\n", "rc('axes', linewidth=2) \n", "fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(7, 7), dpi=100)\n", "\n", "# And plot the data\n", "ax1.imshow(image1, cmap='gray', origin='lower', norm=norm)\n", "ax1.set_title('Before Straylight')\n", "ax1.set_xlabel('X pixel')\n", "ax1.set_ylabel('Y pixel')\n", "\n", "ax2.imshow(image2, cmap='gray', origin='lower', norm=norm)\n", "ax2.set_title('After Straylight')\n", "ax2.set_xlabel('X pixel')" ] }, { "cell_type": "markdown", "id": "953bae61", "metadata": {}, "source": [ "Figure 3: 2d detector images of MIRI data before and after background subtraction. Differences are not noticeable to the eye." ] }, { "cell_type": "code", "execution_count": null, "id": "785a3c45", "metadata": {}, "outputs": [], "source": [ "# It's pretty hard to see the difference by eye, so let's plot a horizontal cut across the images\n", "plt.plot(image1[512, 0:400], label='Original')\n", "plt.plot(image2[512, 0:400], label='Straylight Corrected')\n", "plt.plot(image1[512, 0:400] - image2[512, 0:400], label='Straylight Signal')\n", "plt.yscale('log')\n", "plt.xlabel('Detector X pixel')\n", "plt.ylabel('Flux')\n", "plt.legend()" ] }, { "cell_type": "markdown", "id": "9c3c7a78", "metadata": {}, "source": [ "Figure 4: Horizontal cut through 2d MIRI data before and after straylight correction. Note that the fitted straylight signal shows some peaks around the slice edges that are probably unphysical. More work on this step will be taking place during commissioning as we learn more about what the real straylight signal looks like on orbit." ] }, { "cell_type": "code", "execution_count": null, "id": "81562c8b", "metadata": {}, "outputs": [], "source": [ "# Close out our files\n", "hdu1.close()\n", "hdu2.close()" ] }, { "cell_type": "markdown", "id": "cff79275", "metadata": {}, "source": [ "### 4.7-Fringe Flat (MIRI-only) ###\n", "\n", "
\n", "This crucial step is the first pipeline correction for the strong periodic amplitude modulation (i.e., fringing) that occurs in the MIRI detectors due to internal reflections within the detectors. In this step, the pipeline simply divides by a reference fringe flatfield to make a first-order correction to the data. In detail, the fringing signal will depend on the geometry of sources within the scene, and thus there will be a residual fringe correction later in the pipeline too.\n", "\n", "See https://jwst-pipeline.readthedocs.io/en/latest/jwst/fringe/index.html\n", "
" ] }, { "cell_type": "code", "execution_count": null, "id": "6ee10116", "metadata": { "scrolled": true }, "outputs": [], "source": [ "# Call the step, specifying that we want results saved into the spec2_dir directory\n", "for file in strayfiles:\n", " FringeStep.call(file, save_results=True, output_dir=spec2_dir)" ] }, { "cell_type": "code", "execution_count": null, "id": "a4b78db0", "metadata": {}, "outputs": [], "source": [ "# Look for our fringestep.fits files produced by the fringe flat step\n", "sstring = spec2_dir + 'det*fringestep.fits'\n", "fringefiles = sorted(glob.glob(sstring))\n", "# And print them out so that we can see them\n", "fringefiles" ] }, { "cell_type": "code", "execution_count": null, "id": "666c567a", "metadata": {}, "outputs": [], "source": [ "# Let's show what the SCI extension of the first file before/after application of the fringe flat looks like\n", "# We'll zoom in on a region of the detector to make the results more clear\n", "hdu1 = fits.open(strayfiles[0])\n", "image1 = hdu1['SCI'].data\n", "hdu2 = fits.open(fringefiles[0])\n", "image2 = hdu2['SCI'].data\n", "\n", "# Use a classic ZScale normalization\n", "norm = ImageNormalize(image1, interval=ZScaleInterval(), stretch=LinearStretch())\n", "\n", "rc('axes', linewidth=2) \n", "fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(7, 7), dpi=100)\n", "\n", "# And plot the data\n", "ax1.imshow(image1, cmap='gray', origin='lower', norm=norm)\n", "ax1.set_title('Before Fringe Flat')\n", "ax1.set_xlabel('X pixel')\n", "ax1.set_ylabel('Y pixel')\n", "ax1.set_xlim(0, 250)\n", "ax1.set_ylim(0, 250)\n", "\n", "ax2.imshow(image2, cmap='gray', origin='lower', norm=norm)\n", "ax2.set_title('After Fringe Flat')\n", "ax2.set_xlabel('X pixel')\n", "ax2.set_xlim(0, 250)\n", "ax2.set_ylim(0, 250)" ] }, { "cell_type": "markdown", "id": "320d6479", "metadata": {}, "source": [ "Figure 5: MIRI 2d data before and after application of the static reference fringe flat; note that the periodic modulation in the traces has been removed (by construction perfectly for these simulated data)." ] }, { "cell_type": "code", "execution_count": null, "id": "52229b9d", "metadata": {}, "outputs": [], "source": [ "# Close our files behind us\n", "hdu1.close()\n", "hdu2.close()" ] }, { "cell_type": "markdown", "id": "89f7f576", "metadata": {}, "source": [ "### 4.8-Path-Loss Correction (NIRSpec-only) ###\n", "\n", "
\n", "This is a NIRSpec and NIRISS-only step that estimates flux losses from the data from such effects as light scattering outside the grating. The relevant correction factor is determine for both point sources and uniformly extended sources by interpolation from a reference file given a known point source location.\n", "\n", "For MIRI these corrections are folded into the flux calibration step for typical scenes, and a detector-based point source optimal extraction method that will be available sometime in Cycle 1.\n", "\n", "See https://jwst-pipeline.readthedocs.io/en/latest/jwst/pathloss/index.html\n", "
" ] }, { "cell_type": "markdown", "id": "0be15f65", "metadata": {}, "source": [ "### 4.9-Flux Calibration ###\n", "\n", "
\n", "This step is a crucial part of the pipeline in which count-rate pixel slopes are converted to physical radiometric units by multiplying by reference files derived from observations of spectrophotometric standard sources. In the case of the MRS, the differences in effective pixel solid angle are also accounted for in this step.\n", "\n", "See https://jwst-pipeline.readthedocs.io/en/latest/jwst/photom/index.html\n", "
" ] }, { "cell_type": "code", "execution_count": null, "id": "ae4c30c6", "metadata": { "scrolled": true }, "outputs": [], "source": [ "# Call the step, specifying that we want results saved into the spec2_dir directory\n", "for file in fringefiles:\n", " PhotomStep.call(file, save_results=True, output_dir=spec2_dir)" ] }, { "cell_type": "code", "execution_count": null, "id": "7aab1eae", "metadata": {}, "outputs": [], "source": [ "# Look for our photomstep.fits files produced by the photometric calibration step\n", "sstring = spec2_dir + 'det*photomstep.fits'\n", "photomfiles = sorted(glob.glob(sstring))\n", "# And print them out so that we can see them\n", "photomfiles" ] }, { "cell_type": "markdown", "id": "7a19a31d", "metadata": {}, "source": [ "Ordinarily we'd run the Spec2 pipeline as a pipeline rather than individual steps though, in which\n", "case the final outputs would have the extension _cal.fits rather than _photomstep.fits
\n", "Just to make the rest of this notebook more typical, we'll rename the _photomstep.fits fits to _cal.fits files" ] }, { "cell_type": "code", "execution_count": null, "id": "1e52be37", "metadata": {}, "outputs": [], "source": [ "# Rename photomstep to cal files\n", "calfiles = photomfiles.copy()\n", "for ii in range(0, len(photomfiles)):\n", " calfiles[ii] = str.replace(photomfiles[ii], 'photomstep', 'cal')\n", " thisphotomfile = photomfiles[ii]\n", " thiscalfile = calfiles[ii]\n", " %mv -f $thisphotomfile $thiscalfile" ] }, { "cell_type": "code", "execution_count": null, "id": "80f28087", "metadata": {}, "outputs": [], "source": [ "# Let's show what the SCI extension of the first file before/after application of the flux calibration looks like\n", "# We'll zoom in on a region of the detector to make the results more clear\n", "hdu1 = fits.open(fringefiles[0])\n", "image1 = hdu1['SCI'].data\n", "hdu2 = fits.open(calfiles[0])\n", "image2 = hdu2['SCI'].data\n", "\n", "# Use a classic ZScale normalization\n", "norm1 = ImageNormalize(image1, interval=ZScaleInterval(), stretch=LinearStretch())\n", "norm2 = ImageNormalize(image2, interval=ZScaleInterval(), stretch=LinearStretch())\n", "\n", "rc('axes', linewidth=2) \n", "fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(7, 7), dpi=100)\n", "\n", "# And plot the data\n", "ax1.imshow(image1, cmap='gray', origin='lower', norm=norm1)\n", "ax1.set_title('Before Photom')\n", "ax1.set_xlabel('X pixel')\n", "ax1.set_ylabel('Y pixel')\n", "ax1.set_xlim(0, 250)\n", "ax1.set_ylim(0, 250)\n", "\n", "ax2.imshow(image2, cmap='gray', origin='lower', norm=norm2)\n", "ax2.set_title('After Photom')\n", "ax2.set_xlabel('X pixel')\n", "ax2.set_xlim(0, 250)\n", "ax2.set_ylim(0, 250)" ] }, { "cell_type": "markdown", "id": "a0a07cbc", "metadata": {}, "source": [ "Figure 6: MIRI 2d data before and after application of the photometric calibration reference file." ] }, { "cell_type": "code", "execution_count": null, "id": "a77669e7", "metadata": {}, "outputs": [], "source": [ "# Close our files behind us\n", "hdu1.close()\n", "hdu2.close()" ] }, { "cell_type": "markdown", "id": "01246568", "metadata": {}, "source": [ "### 4.10-Quicklook Cube Building ###\n", "\n", "
\n", "Now that we've got fully calibrated 2d detector-level data (\\_photom.fits if we've run the pipeline step by step, \\_cal.fits if we simply ran the spec2 pipeline as a whole) the pipeline builds some simple data cubes from each exposure. These are multi-band data cubes that combine data from both channels on a given detector into a single data cube with a large blank region in the middle corresponding to the gap between the bands (e.g., it would create a 1A+2A data cube from the SHORT wavelength detector with the SHORT grating setting, and have an empty gap corresponding to the 1B+1C wavelength range).\n", "\n", "We won't generally want to interact with these cubes much; they're big, unwieldy, and since they're for single exposures they aren't well sampled. However, at the moment they're used to populate the background flux information that will be used in the Master Background step in the Spec3 pipeline and are thus necessary to produce.\n", "\n", "See https://jwst-pipeline.readthedocs.io/en/latest/jwst/cube_build/index.html\n", "
" ] }, { "cell_type": "code", "execution_count": null, "id": "701f0492", "metadata": {}, "outputs": [], "source": [ "# Look for our _cal.fits files\n", "sstring = spec2_dir + 'det*cal.fits'\n", "calfiles = sorted(glob.glob(sstring))\n", "# And print them out so that we can see them\n", "calfiles" ] }, { "cell_type": "code", "execution_count": null, "id": "616bf91a", "metadata": { "scrolled": true }, "outputs": [], "source": [ "# Call the step, specifying that we want results saved into the spec2_dir directory\n", "\n", "# If rerunning long pipeline steps, actually run the step\n", "if (redolong is True):\n", " for file in calfiles:\n", " # Note that the 'multi' makes this function like in the spec2 pipeline\n", " CubeBuildStep.call(file, output_type='multi', save_results=True, output_dir=spec2_dir)\n", " \n", "# Otherwise, just copy cached outputs into our output directory structure\n", "else:\n", " sstring = os.path.join(spec2_cachedir, 'det*s3d.fits')\n", " files = sorted(glob.glob(sstring))\n", " for file in files:\n", " outfile = str.replace(file, spec2_cachedir, spec2_dir)\n", " shutil.copy(file, outfile)" ] }, { "cell_type": "code", "execution_count": null, "id": "c1668df9", "metadata": { "scrolled": true }, "outputs": [], "source": [ "# Look for our intermediate cubes produced by the cube building step\n", "sstring = spec2_dir + 'det*s3d.fits'\n", "cubefiles = sorted(glob.glob(sstring))\n", "# And print them out so that we can see them\n", "cubefiles" ] }, { "cell_type": "code", "execution_count": null, "id": "e7026c7e", "metadata": {}, "outputs": [], "source": [ "# Now let's display a couple of these cubes\n", "hdu1 = fits.open(cubefiles[0])\n", "cube1 = hdu1['SCI'].data\n", "hdu2 = fits.open(cubefiles[1])\n", "cube2 = hdu2['SCI'].data\n", "\n", "# Use a linear stretch\n", "norm = ImageNormalize(cube1[0, :, :], vmin=-200, vmax=1e3, stretch=LinearStretch())\n", "\n", "rc('axes', linewidth=2) \n", "fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(7, 7), dpi=100)\n", "\n", "# And plot the data\n", "ax1.imshow(cube1[0, :, :], cmap='gray', norm=norm, origin='lower')\n", "ax1.set_title('Dither #1')\n", "ax1.set_xlabel('X pixel')\n", "ax1.set_ylabel('Y pixel')\n", "\n", "ax2.imshow(cube2[0, :, :], cmap='gray', norm=norm, origin='lower')\n", "ax2.set_title('Dither #2')\n", "ax2.set_xlabel('X pixel')" ] }, { "cell_type": "markdown", "id": "17d1e6f1", "metadata": {}, "source": [ "Figure 7: MIRI rectified 3d data cubes for two different dither positions" ] }, { "cell_type": "markdown", "id": "524eed1f", "metadata": {}, "source": [ "### 4.11-Quicklook Spectral Extraction ###\n", "\n", "
\n", "Now that we've got rough per-exposure data cubes, the pipeline will do a 1d spectral extraction from those cubes. If it thinks the target is an extended source (see Source Type Identification above) it will sum up the flux from the entire field of view; if it thinks the target is a point source it will perform aperture photometry at the expected location of the source (set by the TARG_RA and TARG_DEC keywords in the primary FITS header). At the moment there's no way to override this location short of hacking the keyword values in the FITS header, but this capability should be available soon.\n", "\n", "See https://jwst-pipeline.readthedocs.io/en/latest/jwst/extract_1d/index.html\n", "
" ] }, { "cell_type": "code", "execution_count": null, "id": "dd2c2c4f", "metadata": {}, "outputs": [], "source": [ "Extract1dStep.spec" ] }, { "cell_type": "code", "execution_count": null, "id": "0361a94b", "metadata": { "scrolled": true }, "outputs": [], "source": [ "# Call the step, specifying that we want results saved into the spec2_dir directory\n", "for file in cubefiles:\n", " Extract1dStep.call(file, save_results=True, output_dir=spec2_dir)" ] }, { "cell_type": "code", "execution_count": null, "id": "12df3458", "metadata": {}, "outputs": [], "source": [ "# Look for our intermediate 1d spectra\n", "sstring = spec2_dir + 'det*extract1dstep.fits'\n", "specfiles = sorted(glob.glob(sstring))\n", "# And print them out so that we can see them\n", "specfiles" ] }, { "cell_type": "code", "execution_count": null, "id": "05d06db8", "metadata": {}, "outputs": [], "source": [ "# Let's look at one of them\n", "hdu = fits.open(specfiles[0])\n", "spec = hdu['EXTRACT1D']\n", "\n", "plt.plot(spec.data['WAVELENGTH'], spec.data['FLUX'])\n", "plt.xlabel('Wavelength (micron)')\n", "plt.ylabel('Flux (Jy)')" ] }, { "cell_type": "markdown", "id": "286ca790", "metadata": {}, "source": [ "Figure 8: Extracted 1d spectra from the quick-look per-exposure data cube. Note that in this case we've got spectra in 1A and 2A together in the same spectrum with zeros in between where there is no coverage." ] }, { "cell_type": "code", "execution_count": null, "id": "e0ce9c2e", "metadata": {}, "outputs": [], "source": [ "# In case you forget the units, they're in the FITS header\n", "print(spec.header['TTYPE1'], spec.header['TUNIT1'])\n", "print(spec.header['TTYPE2'], spec.header['TUNIT2'])" ] }, { "cell_type": "code", "execution_count": null, "id": "0ad0e51b", "metadata": {}, "outputs": [], "source": [ "# Close our files behind us\n", "hdu.close()" ] }, { "cell_type": "code", "execution_count": null, "id": "30d5f133", "metadata": {}, "outputs": [], "source": [ "# Print out the time benchmark\n", "time1 = time.perf_counter()\n", "print(f\"Runtime so far: {time1 - time0:0.4f} seconds\")" ] }, { "cell_type": "markdown", "id": "67b385b8", "metadata": {}, "source": [ "We're now done with the Spec2 pipeline!" ] }, { "cell_type": "markdown", "id": "e69364c0", "metadata": {}, "source": [ "
" ] }, { "cell_type": "markdown", "id": "a6059697", "metadata": {}, "source": [ "5.-Spec3 Pipeline \n", "------------------\n", "\n", "
\n", "In this section we process our simulated data through the Spec3 pipeline in order to produce Lvl3 data products (i.e., calibrated combined data cubes and extracted 1d spectra).\n", "\n", "Instead of calling each step in Spec3 individually, we will instead call the overall Spec3 pipeline with various steps turned on/off in order to demonstrate the operation of each of these steps. In part this is to reduce the number of intermediate files that get produced- in pipeline mode the output from one step feeds transparently into the next without needing to produce additional output files. In addition, most of the steps in Spec3 are difficult to evaluate on their own merits and their effects only become clearly visible once the resulting data have been fed through cube building to produce a final 3d data product.\n", "\n", "Spec3 can also become significantly more complex than previous stages in that it can process data from multiple different wavelength bands together. In the present notebook we focus only on the basic operation of the Spec3 steps for a single wavelength band- the interaction between many different bands is treated in greater detail in MRS Calibration Notebook #2 (JWebbinar Series #5 on IFU spectroscopy).\n", "\n", "See https://jwst-pipeline.readthedocs.io/en/latest/jwst/pipeline/calwebb_spec3.html\n", "
" ] }, { "cell_type": "markdown", "id": "ece8213e", "metadata": {}, "source": [ "### 5.1-Association Files ###\n", "\n", "
\n", "The Spec3 pipeline is the first place where we really have to deal with many files at the same time and how they interact with each other (i.e., background observations, dithered observations, etc). As such, we need to create an 'Association File' describing these files and how they should be treated by the pipeline. These files are still under some development, and the recommended methods of using them will likely change before Cycle 1. We'll therefore define a function to create a very simple association file that will treat all exposures it is given as science exposures to be combined together.\n", "\n", "See https://jwst-pipeline.readthedocs.io/en/latest/jwst/associations/index.html\n", "
" ] }, { "cell_type": "code", "execution_count": null, "id": "629dd5dc", "metadata": { "scrolled": false }, "outputs": [], "source": [ "# Define a useful function to write out a Lvl3 association file from an input list\n", "def writel3asn(files, asnfile, prodname, **kwargs):\n", " # Define the basic association of science files\n", " asn = afl.asn_from_list(files, rule=DMS_Level3_Base, product_name=prodname)\n", " # Add any background files to the association\n", " if ('bg' in kwargs):\n", " for bgfile in kwargs['bg']:\n", " asn['products'][0]['members'].append({'expname': bgfile, 'exptype': 'background'})\n", " # Write the association to a json file\n", " _, serialized = asn.dump()\n", " with open(asnfile, 'w') as outfile:\n", " outfile.write(serialized)" ] }, { "cell_type": "markdown", "id": "b1e57612", "metadata": {}, "source": [ "### 5.2-Moving Target WCS ###\n", "\n", "
\n", "This step modifies the WCS embedded in the individual calibrated frames to account for motion of the target between individual frames (e.g., for solar system objects) such that the WCS for all frames is centered at the average location of the target within the association.\n", "\n", "Since this step represents a rather specialized use case we will skip it for purposes of this notebook.\n", "\n", "See https://jwst-pipeline.readthedocs.io/en/latest/jwst/assign_mtwcs/main.html\n", "
" ] }, { "cell_type": "markdown", "id": "5b04b907", "metadata": {}, "source": [ "### 5.3-Master Background Subtraction ###\n", "\n", "
\n", "As described above, the MIRI MRS is expected to see a significant thermal and zodiacal background that needs to be subtracted from the science data. In the Master Background step, we construct a model of the background and subtract it from the data. This can produce significantly higher SNR results than simply doing a pixel-by-pixel subtraction of the background reference data, but it will not be known until flight if the results from this approach will be satisfactory. As such, this step of the pipeline is still undergoing some development and is focused around processing of extended source observations with dedicated background pointings.\n", "\n", "Since this is not the use case that we're exploring in this notebook, we will skip this step for now (an annular background will be subtracted in the Extract1D stage).\n", "\n", "See https://jwst-pipeline.readthedocs.io/en/latest/jwst/master_background/index.html\n", "
" ] }, { "cell_type": "markdown", "id": "85b7e316", "metadata": {}, "source": [ "### 5.4-Residual Background Matching (MIRI-only) ###\n", "\n", "
\n", "Since the thermal background at MIRI wavelength is significant, there is the possibility that it could vary in undesirable ways between exposures. As such, a single background image would be insufficient to bring all exposures to the same level. The residual background matching step therefore ensures that the low-order background is consistent across all exposures.\n", "\n", "See https://jwst-pipeline.readthedocs.io/en/latest/jwst/mrs_imatch/index.html\n", "
" ] }, { "cell_type": "code", "execution_count": null, "id": "1dc590cf", "metadata": { "scrolled": false }, "outputs": [], "source": [ "# Read in two calibrated frames\n", "# Look for our _rate.fits files produced by the Detector1 pipeline\n", "sstring = spec2_dir + 'det*cal.fits'\n", "calfiles = sorted(glob.glob(sstring))\n", "calfiles" ] }, { "cell_type": "markdown", "id": "b5b30085", "metadata": {}, "source": [ "In these simulations the background doesn't vary, but let's pretend that it did.\n", "We'll crudely hack one of the files to mimic a background level shift." ] }, { "cell_type": "code", "execution_count": null, "id": "2e8db49e", "metadata": { "scrolled": false }, "outputs": [], "source": [ "hdu = fits.open(calfiles[0])\n", "hdu['SCI'].data += 300 # Hack this first file to add a pedestal offset\n", "hdu.writeto(spec2_dir+'rbm_test.fits', overwrite=True)\n", "hdu.close()" ] }, { "cell_type": "code", "execution_count": null, "id": "056bc1ec", "metadata": { "scrolled": true }, "outputs": [], "source": [ "# Now we'll create an association file including this hacked exposure\n", "testfiles = calfiles.copy()\n", "testfiles[0] = spec2_dir + 'rbm_test.fits'\n", "writel3asn(testfiles, 'rbm.json', 'rbm')\n", "\n", "# And run it through cube building (we'll just build a cube for the Ch2 data as an example to save time), calling the output\n", "# file 'rbm_before'\n", "cb = CubeBuildStep()\n", "\n", "# If rerunning long pipeline steps, actually run the step\n", "if (redolong is True):\n", " cb.call('rbm.json', channel='2', save_results=True, output_dir=spec3_dir, output_file='rbm_before')\n", "# Otherwise, just copy cached outputs into our output directory structure\n", "else:\n", " sstring = os.path.join(spec3_cachedir, 'rbm_before*s3d.fits')\n", " files = sorted(glob.glob(sstring))\n", " for file in files:\n", " outfile = str.replace(file, spec3_cachedir, spec3_dir)\n", " shutil.copy(file, outfile)" ] }, { "cell_type": "code", "execution_count": null, "id": "bcb80f72", "metadata": { "scrolled": true }, "outputs": [], "source": [ "# And for comparison we'll run this association through the Spec3 pipeline with just Residual Background Matching and Cube Build,\n", "# calling the output 'rbm_after'\n", "spec3 = Spec3Pipeline()\n", "spec3.output_dir = spec3_dir\n", "spec3.save_results = True\n", "spec3.master_background.skip = True\n", "spec3.outlier_detection.skip = True\n", "spec3.cube_build.channel = '2'\n", "spec3.cube_build.output_file = 'rbm_after'\n", "spec3.extract_1d.skip = True\n", "\n", "# If rerunning long pipeline steps, actually run the step\n", "if (redolong is True):\n", " spec3('rbm.json')\n", "# Otherwise, just copy cached outputs into our output directory structure\n", "else:\n", " sstring = os.path.join(spec3_cachedir, 'rbm_after*s3d.fits')\n", " files = sorted(glob.glob(sstring))\n", " for file in files:\n", " outfile = str.replace(file, spec3_cachedir, spec3_dir)\n", " shutil.copy(file, outfile)" ] }, { "cell_type": "code", "execution_count": null, "id": "cbf7f162", "metadata": { "scrolled": false }, "outputs": [], "source": [ "# Now let's compare what the SCI data of the two resulting cubes looks like\n", "hdu1 = fits.open(spec3_dir + 'rbm_before_ch2-long_s3d.fits')\n", "cube1 = hdu1['SCI'].data\n", "hdu2 = fits.open(spec3_dir + 'rbm_after_ch2-long_s3d.fits')\n", "cube2 = hdu2['SCI'].data\n", "\n", "# Use a classic ZScale normalization\n", "norm = ImageNormalize(cube1[0, :, :], interval=ZScaleInterval(), stretch=LinearStretch())\n", "\n", "rc('axes', linewidth=2) \n", "fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(7, 7), dpi=100)\n", "\n", "# And plot the data\n", "ax1.imshow(cube1[0, :, :], cmap='gray', norm=norm, origin='lower')\n", "ax1.set_title('No BG matching')\n", "ax1.set_xlabel('X pixel')\n", "ax1.set_ylabel('Y pixel')\n", "\n", "ax2.imshow(cube2[0, :, :], cmap='gray', norm=norm, origin='lower')\n", "ax2.set_title('With BG matching')\n", "ax2.set_xlabel('X pixel')" ] }, { "cell_type": "markdown", "id": "729ee8e0", "metadata": {}, "source": [ "Figure 9: Data cubes built with and without residual background matching. The example with background matching clear looks much better! We should therefore be sure to use this step if such patterns are visible in the data. Note that the pipeline doesn't know what the 'true' background level should be, so it brings everything to the same average level in spatially overlapping regions (i.e., if you're mosaicing a nebula it won't try to make all parts of the nebular have the same flux, just all exposures of a given point in the nebula should be the same)." ] }, { "cell_type": "code", "execution_count": null, "id": "9bfa7899", "metadata": {}, "outputs": [], "source": [ "# Close our files behind us\n", "hdu1.close()\n", "hdu2.close()" ] }, { "cell_type": "markdown", "id": "ee0e86b4", "metadata": {}, "source": [ "### 5.5-Outlier Detection ###\n", "\n", "
\n", "In general, cosmic rays should have been flagged in the Detector1 pipeline as jumps in the readout ramp and corrected when determining the final slope image. However, sometimes fainter cosmic rays slip through Detector1, and other outliers (due to bad pixels, etc) can be present that cause artifacts in the final data cubes. The Outlier Detection routine is designed to identify and flag these artifacts.\n", "\n", "In brief, it builds a data cube from each individual exposure on a common data cube grid, and median combines the cubes to construct a cleaned cube. This cleaned cube is mapped back to the 2d calibrated detector data of each exposure, and outliers in individual frames with respect to this mapped-back data are flagged. When we next do cube building, it will use this updated flagging to construct our high-quality combined data cube.\n", "\n", "Since the IFUs are significantly undersampled however, we need to be careful not to accidentally flag point source as outliers since the fluxes can change significantly between different pointings. Note also that outlier detection cannot detect outliers when there are just 2 or fewer frames; therefore even in a four-point dither pattern it will only be able to work in the central regions where all four exposures overlap.\n", "\n", "See https://jwst-pipeline.readthedocs.io/en/latest/jwst/outlier_detection/index.html\n", "
" ] }, { "cell_type": "markdown", "id": "0bfdb188", "metadata": {}, "source": [ "As for residual background matching, there aren't generally many outliers in mirisim simulated data for us to demonstrate this algorithm on.\n", "Therefore, we'll add some outliers into a copy of the simulated data and show that they get successfully flagged and rejected." ] }, { "cell_type": "code", "execution_count": null, "id": "e315c171", "metadata": {}, "outputs": [], "source": [ "# Read in four calibrated frames\n", "# Look for our _rate.fits files produced by the Detector1 pipeline\n", "sstring = spec2_dir + 'det*cal.fits'\n", "calfiles = sorted(glob.glob(sstring))\n", "# And print them out so that we can see them\n", "calfiles" ] }, { "cell_type": "code", "execution_count": null, "id": "b8e9937e", "metadata": {}, "outputs": [], "source": [ "# We'll crudely hack the first files to introduce a artifact on column 925 and write it back out to an _odhack.fits files\n", "hdu = fits.open(calfiles[0])\n", "data = hdu['SCI'].data\n", "data[:, 925] = 13000\n", "hdu['SCI'].data = data\n", "hdu.writeto(str.replace(calfiles[0], 'cal', 'od_test'), overwrite=True)\n", "hdu.close()" ] }, { "cell_type": "code", "execution_count": null, "id": "89aee526", "metadata": {}, "outputs": [], "source": [ "# Now we'll create an association file including these hacked exposures\n", "testfiles = calfiles.copy()\n", "testfiles[0] = str.replace(calfiles[0], 'cal', 'od_test')\n", "writel3asn(testfiles, 'od.json', 'od')\n", "\n", "# And run it through cube building (we'll just build a cube for the Ch2 data as an example to save time),\n", "# calling the result 'od_before'\n", "cb = CubeBuildStep()\n", "\n", "# If rerunning long pipeline steps, actually run the step\n", "if (redolong is True):\n", " cb.call('od.json', channel='2', save_results=True, output_dir=spec3_dir, output_file='od_before')\n", "# Otherwise, just copy cached outputs into our output directory structure\n", "else:\n", " sstring = os.path.join(spec3_cachedir, 'od_before*s3d.fits')\n", " files = sorted(glob.glob(sstring))\n", " for file in files:\n", " outfile = str.replace(file, spec3_cachedir, spec3_dir)\n", " shutil.copy(file, outfile)" ] }, { "cell_type": "code", "execution_count": null, "id": "9e272ccb", "metadata": {}, "outputs": [], "source": [ "# And for comparison we'll run this association through the Spec3 pipeline with just \n", "# Outlier Detection and Cube Build, saving our final cube as 'od_after'\n", "\n", "# This initial setup is just to make sure that we get the latest parameter reference files\n", "# pulled in for our files. This is a temporary workaround to get around an issue with\n", "# how this pipeline calling method works.\n", "crds_config = Spec3Pipeline.get_config_from_reference('od.json')\n", "spec3 = Spec3Pipeline.from_config_section(crds_config)\n", "\n", "# Now make any modifications to the pipeline specification\n", "spec3.output_dir = spec3_dir\n", "spec3.save_results = True\n", "spec3.master_background.skip = True\n", "spec3.mrs_imatch.skip = True\n", "spec3.outlier_detection.save_intermediate_results = True # We'll write out intermediate files to explore what they look like\n", "spec3.cube_build.channel = '2'\n", "spec3.cube_build.output_file = 'od_after'\n", "spec3.extract_1d.skip = True\n", "\n", "# If rerunning long pipeline steps, actually run the step\n", "if (redolong is True):\n", " spec3('od.json')\n", "# Otherwise, just copy cached outputs into our output directory structure\n", "else:\n", " sstring = os.path.join(spec3_cachedir, 'od_after*s3d.fits')\n", " files = sorted(glob.glob(sstring))\n", " for file in files:\n", " outfile = str.replace(file, spec3_cachedir, spec3_dir)\n", " shutil.copy(file, outfile)\n", " sstring = os.path.join(spec3_cachedir, '*crf.fits')\n", " files = sorted(glob.glob(sstring))\n", " for file in files:\n", " outfile = str.replace(file, spec3_cachedir, spec3_dir)\n", " shutil.copy(file, outfile)" ] }, { "cell_type": "code", "execution_count": null, "id": "398ff580", "metadata": {}, "outputs": [], "source": [ "# Show the image and mask\n", "hdu = fits.open(spec3_dir + 'det_image_seq1_MIRIFUSHORT_12LONGexp1_od_test_a3001_crf.fits')\n", "flux = hdu['SCI'].data\n", "dq = hdu['DQ'].data\n", "\n", "# Use a classic ZScale normalization\n", "norm = ImageNormalize(flux, interval=ZScaleInterval(), stretch=LinearStretch())\n", "\n", "rc('axes', linewidth=2) \n", "fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(7, 7), dpi=100)\n", "\n", "# And plot the data. Highlight a pixel in the bad column with a red X\n", "ax1.imshow(flux, cmap='gray', norm=norm, origin='lower')\n", "ax1.set_title('2d SCI array')\n", "ax1.set_xlabel('X pixel')\n", "ax1.set_ylabel('Y pixel')\n", "ax1.set_xlim(850, 1000)\n", "ax1.set_ylim(0, 150)\n", "ax1.plot([925], [60], 'X', color='red')\n", "\n", "ax2.imshow(dq, cmap='gray', vmin=0, vmax=1, origin='lower')\n", "ax2.set_title('2d DQ array')\n", "ax2.set_xlabel('X pixel')\n", "ax2.set_xlim(850, 1000)\n", "ax2.set_ylim(0, 150)\n", "ax2.plot([925], [60], 'X', color='red')" ] }, { "cell_type": "markdown", "id": "ddfe0c9d", "metadata": {}, "source": [ "Figure 10: SCI and ERR extensions for a MIRI 2d calibrated detector image in which one column has been deliberately set to an erroneous value. Note that the DQ array identifies this erroneous column after the outlier detection step has been run." ] }, { "cell_type": "markdown", "id": "0acfd5dc", "metadata": {}, "source": [ "That DQ array is pretty messy! What is it telling us?\n", "We can see that the column is flagged in our DQ array; let's examine what the DQ flag value is" ] }, { "cell_type": "code", "execution_count": null, "id": "3e257f42", "metadata": {}, "outputs": [], "source": [ "print(dq[60, 925])" ] }, { "cell_type": "markdown", "id": "e30cb3ed", "metadata": {}, "source": [ "That's not very useful on its own, but we can ask the pipeline to tell us what this DQ flag value actually means:" ] }, { "cell_type": "code", "execution_count": null, "id": "194cbd86", "metadata": {}, "outputs": [], "source": [ "dqflags.dqflags_to_mnemonics(dq[60, 925], mnemonic_map=datamodels.dqflags.pixel)" ] }, { "cell_type": "markdown", "id": "432380e7", "metadata": {}, "source": [ "That is, it's now been flagged as an OUTLIER and warned the pipeline not to use it! Now let's compare what the 3d cubes with and without outlier detection look like." ] }, { "cell_type": "code", "execution_count": null, "id": "86661e51", "metadata": {}, "outputs": [], "source": [ "# Close our files behind us\n", "hdu.close()" ] }, { "cell_type": "code", "execution_count": null, "id": "4f9ab6b2", "metadata": {}, "outputs": [], "source": [ "# Cube without outlier rejection\n", "hdu1 = fits.open(spec3_dir + 'od_before_ch2-long_s3d.fits')\n", "flux1 = hdu1['SCI'].data\n", "\n", "# Cube with outlier rejection\n", "hdu2 = fits.open(spec3_dir + 'od_after_ch2-long_s3d.fits')\n", "flux2 = hdu2['SCI'].data\n", "\n", "# Use a classic ZScale normalization\n", "norm = ImageNormalize(flux1, interval=ZScaleInterval(), stretch=LinearStretch())\n", "\n", "rc('axes', linewidth=2) \n", "fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(7, 7), dpi=100)\n", "\n", "# And plot the data. Highlight a pixel in the bad column with a red X\n", "ax1.imshow(flux1[0, :, :], cmap='gray', norm=norm, origin='lower')\n", "ax1.set_title('No Outlier Det.')\n", "ax1.scatter([14], [22], s=580, facecolors='none', edgecolors='r')\n", "\n", "ax2.imshow(flux2[0, :, :], cmap='gray', norm=norm, origin='lower')\n", "ax2.set_title('With Outlier Det.')\n", "ax2.scatter([14], [22], s=580, facecolors='none', edgecolors='r')" ] }, { "cell_type": "markdown", "id": "0bfb2f5f", "metadata": {}, "source": [ "Figure 11: 3d data cubes constructed with and without the outlier detection step. The red circle denotes the location of the outlier that we injected into the data." ] }, { "cell_type": "code", "execution_count": null, "id": "47487310", "metadata": {}, "outputs": [], "source": [ "# Close our files behind us\n", "hdu1.close()\n", "hdu2.close()" ] }, { "cell_type": "markdown", "id": "6fe6742a", "metadata": {}, "source": [ "### 5.6-Cube Creation ###\n", "\n", "
\n", "Now that all of our individual exposures have been flux calibrated, background subtracted/matched, and any outliers flagged for removal, we can finally build a composite data cube from our dithered data! Under the hood we've already been using this code for the background matching and outlier detection steps, but this step is the one that produces a 'final' 3d output product from the pipeline.\n", "\n", "Note, however, that in the example given below we're building cubes directly from the 2d calibration data (\\_cal.fits files), and thus not incorporating any changes from the outlier detection or residual background matching steps. This is why it's much simpler to run Spec3 as a pipeline (as demonstrated in MRS Notebook 2) as the relevant information is transparently passed from step to step.\n", "\n", "A final word of caution: the data cubes created by the JWST pipeline are in SURFACE BRIGHTNESS units (MJy/steradian), not flux units. What that means is that if you intend to sum spectra within an aperture you need to be sure to multiply by the pixel area in steradians first in order to get a spectrum in flux units (the PIXAR_SR keyword can be found in the SCI extension header). This correction is already build into the pipeline Extract1D algorithm.\n", "\n", "See https://jwst-pipeline.readthedocs.io/en/latest/jwst/cube_build/index.html\n", "
" ] }, { "cell_type": "code", "execution_count": null, "id": "4a837a1d", "metadata": { "scrolled": true }, "outputs": [], "source": [ "# Let's create an association file with the calibrated 2d data\n", "writel3asn(calfiles, 'l3.json', 'l3')\n", "\n", "# And run it through cube building (we'll just build a cube for the Ch2 data as an example to save time)\n", "cb = CubeBuildStep()\n", "\n", "# If rerunning long pipeline steps, actually run the step\n", "if (redolong is True):\n", " cb.call('l3.json', channel='2', save_results=True, output_dir=spec3_dir)\n", "# Otherwise, just copy cached outputs into our output directory structure\n", "else:\n", " sstring = os.path.join(spec3_cachedir, 'l3*s3d.fits')\n", " files = sorted(glob.glob(sstring))\n", " for file in files:\n", " outfile = str.replace(file, spec3_cachedir, spec3_dir)\n", " shutil.copy(file, outfile)" ] }, { "cell_type": "markdown", "id": "e61ab0a9", "metadata": {}, "source": [ "We can also build data cubes with a different rotation convention. By default, cubes are built in the 'skyalign' system, in which N is up and E is left. However, we can also build them in the 'ifualign' system, in which the axes of the cube are tied to the physical axes of the IFU footprint; let's call it 'rotated'." ] }, { "cell_type": "code", "execution_count": null, "id": "3a00ed53", "metadata": { "scrolled": true }, "outputs": [], "source": [ "# Build a rotated-frame cube\n", "\n", "cb = CubeBuildStep()\n", "\n", "# If rerunning long pipeline steps, actually run the step\n", "if (redolong is True):\n", " cb.call('l3.json', channel='2', save_results=True, output_dir=spec3_dir, coord_system='ifualign', output_file='rotated')\n", "# Otherwise, just copy cached outputs into our output directory structure\n", "else:\n", " sstring = os.path.join(spec3_cachedir, 'rotated*s3d.fits')\n", " files = sorted(glob.glob(sstring))\n", " for file in files:\n", " outfile = str.replace(file, spec3_cachedir, spec3_dir)\n", " shutil.copy(file, outfile)" ] }, { "cell_type": "code", "execution_count": null, "id": "f01c294b", "metadata": {}, "outputs": [], "source": [ "# Let's take a quick look at these two cubes\n", "# Cube without outlier rejection\n", "hdu1 = fits.open(spec3_dir + 'l3_ch2-long_s3d.fits')\n", "flux1 = hdu1['SCI'].data\n", "\n", "# Cube with outlier rejection\n", "hdu2 = fits.open(spec3_dir + 'rotated_ch2-long_s3d.fits')\n", "flux2 = hdu2['SCI'].data\n", "\n", "# Use a classic ZScale normalization with a logarithmic stretch to make sure that\n", "# we can see the actual cube footprint well\n", "norm = ImageNormalize(flux1, interval=ZScaleInterval(), stretch=LogStretch())\n", "\n", "rc('axes', linewidth=2) \n", "fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(7, 7), dpi=100)\n", "\n", "# And plot the data. Highlight a pixel in the bad column with a red X\n", "ax1.imshow(flux1[0, :, :], cmap='gray', norm=norm, origin='lower')\n", "ax1.set_title('SKYALIGN')\n", "\n", "ax2.imshow(flux2[0, :, :], cmap='gray', norm=norm, origin='lower')\n", "ax2.set_title('IFUALIGN')" ] }, { "cell_type": "markdown", "id": "bceae952", "metadata": {}, "source": [ "Figure 12: MIRI data cubes constructed using the 'skyalign' and 'ifualign' coordinate reference frames. Both have full WCS information embedded in the data so that they can be (e.g.) displayed similarly using tools such as ds9." ] }, { "cell_type": "code", "execution_count": null, "id": "a503135a", "metadata": {}, "outputs": [], "source": [ "# Close our files behind us\n", "hdu1.close()\n", "hdu2.close()" ] }, { "cell_type": "markdown", "id": "c59297ef", "metadata": {}, "source": [ "### 5.7-Spectral Extraction ###\n", "\n", "
\n", "This step extracts a 1-dimensional spectrum from the final composite 3d data cube. As for the spectral extraction on a per-exposure basis at the end of the Spec2 pipeline, it will sum the entire field of view for EXTENDED sources, and perform aperture photometry at the nominal source location for POINT sources (which includes both background subtraction using an annular ring to define the background signal, and aperture correction).\n", "\n", "In practice, for many science cases the 3d data cubes will be much to complex to be distilled down to a single 1-dimensional spectrum. With the exception of genuine point sources, the 1d extracted spectra should thus be viewed as a convenience rather than the end-stage product to be used for scientific analysis.\n", "\n", "See https://jwst-pipeline.readthedocs.io/en/latest/jwst/extract_1d/index.html\n", "
" ] }, { "cell_type": "code", "execution_count": null, "id": "ce744390", "metadata": {}, "outputs": [], "source": [ "# Run the Extract1D step on the final 3d cube\n", "cubefile = spec3_dir + 'l3_ch2-long_s3d.fits'\n", "Extract1dStep.call(cubefile, save_results=True, output_dir=spec3_dir)" ] }, { "cell_type": "code", "execution_count": null, "id": "1e84b034", "metadata": {}, "outputs": [], "source": [ "# Let's look at the result\n", "specfile = spec3_dir + 'l3_ch2-long_extract1dstep.fits'\n", "\n", "# Let's look at one of them\n", "hdu = fits.open(specfile)\n", "spec = hdu['EXTRACT1D']\n", "\n", "plt.plot(spec.data['WAVELENGTH'], spec.data['FLUX'])\n", "plt.xlabel('Wavelength (micron)')\n", "plt.ylabel('Flux (Jy)')" ] }, { "cell_type": "markdown", "id": "a49b1ecb", "metadata": {}, "source": [ "Figure 13: Extracted 1d spectrum for band 2C" ] }, { "cell_type": "code", "execution_count": null, "id": "aa826037", "metadata": {}, "outputs": [], "source": [ "# Close our files behind us\n", "hdu.close()" ] }, { "cell_type": "code", "execution_count": null, "id": "0d63ab12", "metadata": {}, "outputs": [], "source": [ "# Print out the time benchmark\n", "time1 = time.perf_counter()\n", "print(f\"Runtime so far: {time1 - time0:0.4f} seconds\")" ] }, { "cell_type": "markdown", "id": "00c8f014", "metadata": {}, "source": [ "We're now done with the Spec3 pipeline!" ] }, { "cell_type": "markdown", "id": "1435c29c", "metadata": {}, "source": [ "
" ] }, { "cell_type": "markdown", "id": "483f94a9", "metadata": {}, "source": [ "\"stsci_logo\" " ] } ], "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 }