{ "cells": [ { "cell_type": "markdown", "id": "77ca1b9a", "metadata": {}, "source": [ "\"stsci_logo\" " ] }, { "cell_type": "markdown", "id": "7e783dc5", "metadata": {}, "source": [ "\n", "# MIRI MRS Calibration Notebook #2 (Multiband and Parallel Processing Demo; Point Sources) #" ] }, { "cell_type": "markdown", "id": "cc39cee4", "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": "dd8a38c1", "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", " 2.5 [Multithreading Function](#mthread)
\n", "3. [Detector1 Pipeline](#det1)
\n", "4. [Spec2 Pipeline](#spec2)
\n", "5. [Spec3 Pipeline: Default configuration (12 per-band cubes)](#spec3)
\n", "6. [Spec3 Pipeline: Channel configuration (4 per-channel cubes)](#spec3_v2)
\n", "7. [Spec3 Pipeline: Uber configuration (1 all-wavelength cube)](#spec3_v3)
\n", "8. [Examine the output cubes, and compare spectra to the simulation inputs](#examine)
" ] }, { "cell_type": "markdown", "id": "862e6477", "metadata": {}, "source": [ "
" ] }, { "cell_type": "markdown", "id": "c7f7fd72", "metadata": {}, "source": [ "1.-Introduction \n", "------------------" ] }, { "cell_type": "markdown", "id": "aed5986a", "metadata": {}, "source": [ "### 1.1-Purpose of this Notebook ###\n", "\n", "In this notebook we provide a series of realistic example for running the JWST pipeline on large quantities of MIRI MRS data similar to those that will be obtained by many Cycle 1 observing programs. In particular, we focus on handling issues related to multi-band ditheredobservations spanning the MRS wavelength range, both as regards efficient reduction times and data display issues. We therefore do not discuss the purpose and input/output of each individual pipeline step, but focus on ways to modify the final pipeline output for the best scientific utility. For a walkthrough of individual steps, see Notebook #1. Similarly, this will cover the baseline pipeline as it exists in October 2021 and does not touch advanced algorithms or steps still under development.\n", "\n", "This notebook is specifically for POINT-SOURCE observations; extended source observations are somewhat different and will be presented in MRS Notebook #3.\n", "\n", "We will start with a simple simulated MRS observation (created using mirisim: https://wiki.miricle.org/Public/MIRISim_Public) with a dedicated background, 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", "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": "20c382e7", "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 entire MRS waveband from 5-28 micron (i.e., both MIRI detectors, with the SHORT/MEDIUM/LONG grating configurations). 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": "de1a66d6", "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 updated 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": "fddc9773", "metadata": {}, "source": [ "
" ] }, { "cell_type": "markdown", "id": "d0dbbd9f", "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. This includes a variety of multiprocessing functions that allow us to parallelize pipeline reductions of many exposures using the multiple cores now standard in many computers.\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 40 minutes (8 CPUs; vs 1.5 hours with one CPU) vs 10 seconds.\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 6 hours with 8 CPUs." ] }, { "cell_type": "markdown", "id": "17e81629", "metadata": {}, "source": [ "### 2.1-CRDS Context ###" ] }, { "cell_type": "markdown", "id": "acc46ca5", "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": "115ca6c7", "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": "d9cdcfba", "metadata": {}, "source": [ "### 2.2-Python Imports ###" ] }, { "cell_type": "code", "execution_count": null, "id": "f5a01ba3", "metadata": {}, "outputs": [], "source": [ "# Import packages for multiprocessing. These won't be used on the online demo, but can be\n", "# very useful for local data processing unless/until they get integrated natively into\n", "# the cube building code. These need to be imported before anything else.\n", "\n", "import multiprocessing\n", "multiprocessing.set_start_method('fork')\n", "from multiprocessing import Pool\n", "#Modify the path to a directory on your machine\n", "import os\n", "os.environ[\"CRDS_PATH\"] = \"/path/to/my/folder/\"\n", "os.environ[\"CRDS_SERVER_URL\"] = \"https://jwst-crds-pub.stsci.edu\"\n", "\n", "# Set the maximum number of processes to spawn based on available cores\n", "usage = 'none' # Either 'none' (single thread), 'quarter', 'half', or 'all' available cores" ] }, { "cell_type": "code", "execution_count": null, "id": "1e5c73b7", "metadata": {}, "outputs": [], "source": [ "# Multiprocessing can snarl up the Jupyter notebook display, so if we're doing it we need to\n", "# tell the pipeline to send its output to a log file instead of the screen BEFORE we import\n", "# the pipeline.\n", "\n", "# If not multiprocessing, remove the local .cfg file that the pipeline will look for\n", "if (usage == 'none'):\n", " if os.path.exists('stpipe-log.cfg'):\n", " os.remove('stpipe-log.cfg')\n", "\n", "# If multiprocessing, create the .cfg file that the pipeline will look for\n", "else:\n", " # Create a .cfg file in the working directory to specify where to send pipeline output to\n", " print('[*]', file=open('stpipe-log.cfg', \"w\"))\n", " print('handler = file:pipeline.log', file=open('stpipe-log.cfg', \"a\"))\n", " print('level = INFO', file=open('stpipe-log.cfg', \"a\"))" ] }, { "cell_type": "code", "execution_count": null, "id": "01c552cc", "metadata": {}, "outputs": [], "source": [ "num_cores = multiprocessing.cpu_count()\n", "if usage == 'quarter':\n", " maxp = num_cores // 4 or 1\n", "elif usage == 'half':\n", " maxp = num_cores // 2 or 1\n", "elif usage == 'all':\n", " maxp = num_cores\n", "else:\n", " maxp = 1\n", "\n", "print('We will use '+str(maxp)+' CPUs in this run')" ] }, { "cell_type": "code", "execution_count": null, "id": "7fecce2a", "metadata": {}, "outputs": [], "source": [ "# Now 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": "a5a835ac", "metadata": {}, "outputs": [], "source": [ "# Basic system utilities for interacting with files\n", "import glob\n", "import sys\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": "8fb78d51", "metadata": {}, "outputs": [], "source": [ "# Import the base JWST package\n", "import jwst" ] }, { "cell_type": "code", "execution_count": null, "id": "c853939c", "metadata": {}, "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", "# 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\n", "\n", "import stcal.ramp_fitting.utils as utils # Utilities for handling multiprocessing" ] }, { "cell_type": "markdown", "id": "f6c32ec1", "metadata": {}, "source": [ "### 2.3-Data I/O Directories ###" ] }, { "cell_type": "markdown", "id": "d941c6a8", "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 (~ 3 GB) it can also take a few minutes or more to download, and can be downloaded through any other means and extracted into the target data directory." ] }, { "cell_type": "code", "execution_count": null, "id": "597d707a", "metadata": {}, "outputs": [], "source": [ "# input_dir should point to where you installed the data from Box\n", "input_dir = './MRS_data2/'\n", "\n", "ziplink = 'https://stsci.box.com/shared/static/94hv0t7jmourchlhpf3ysxoq8p3rxlj7.zip'\n", "zipfilename = 'mrs_data2.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": "1c40ca4b", "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": "50bbdfb4", "metadata": {}, "outputs": [], "source": [ "# output_dir should point to where you want the output results to go\n", "output_dir = './MRS_output2/'\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": "8ef8c2c9", "metadata": {}, "source": [ "### 2.4-Reprocessing Flag ###" ] }, { "cell_type": "markdown", "id": "ad7ff1e3", "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": "a9cddd63", "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 = True\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": "02b61432", "metadata": {}, "source": [ "### 2.5-Multithreading Function ###" ] }, { "cell_type": "code", "execution_count": null, "id": "51271da5", "metadata": {}, "outputs": [], "source": [ "# Define a function that can oversee the multiprocessing\n", "# We'll just pass this a function name and a set of inputs, and it'll figure out the\n", "# multithreading accordingly.\n", "def runmany(step, filenames):\n", " if __name__ == '__main__':\n", " p = Pool(maxp)\n", " res = p.map(step, filenames)\n", " p.close()\n", " p.join()" ] }, { "cell_type": "markdown", "id": "266a1b0f", "metadata": {}, "source": [ "
" ] }, { "cell_type": "markdown", "id": "202aface", "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 will turn off the reference pixel subtraction step as it does not handle mirisim data well.\n", "\n", "See https://jwst-pipeline.readthedocs.io/en/latest/jwst/pipeline/calwebb_detector1.html\n", "
" ] }, { "cell_type": "code", "execution_count": null, "id": "805cb62e", "metadata": {}, "outputs": [], "source": [ "# Start a timer to keep track of runtime\n", "time0 = time.perf_counter()" ] }, { "cell_type": "code", "execution_count": null, "id": "9f3660af", "metadata": {}, "outputs": [], "source": [ "# First we'll define a function that will call the detector1 pipeline with our desired set of parameters\n", "# We won't enumerate the individual steps\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": "c1383d60", "metadata": {}, "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": "d44991e3", "metadata": {}, "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", " runmany(rundet1, simfiles)\n", " \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": "code", "execution_count": null, "id": "e3e79df2", "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": "9f1151ba", "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 and 1d spectra). We will not look into each individual step in detail- see MRS Notebook #1 for a guide to that.\n", "\n", "We will skip the 'background' step as this is largely a placeholder for in case pixel-by-pixel differencing of science and background pointings turns out to be necessary on orbit.\n", "\n", "Likewise, we will skip the 'straylight' step as this is unnecessary for simulated data and can sometimes introduce artifacts at present.\n", "\n", "See https://jwst-pipeline.readthedocs.io/en/latest/jwst/pipeline/calwebb_spec2.html\n", "
" ] }, { "cell_type": "code", "execution_count": null, "id": "af61d821", "metadata": {}, "outputs": [], "source": [ "# Define a function that will call the spec2 pipeline with our desired set of parameters\n", "# We'll list the individual steps just to make it clear what's running\n", "def runspec2(filename):\n", " spec2 = Spec2Pipeline()\n", " spec2.output_dir = spec2_dir\n", " \n", " spec2.assign_wcs.skip = False # Derives the world coordinate solution- never skip this!\n", " spec2.bkg_subtract.skip = True # Performs pixel-wise background subtraction- placeholder in case necessary\n", " spec2.flat_field.skip = False # Applies a pixel flatfield (currently all unity)\n", " spec2.srctype.skip = False # Guesses at the source type (point/extended) based on APT inputs and dither information\n", " spec2.straylight.skip = True # Model and subtraction straylight (may want to skip)\n", " spec2.fringe.skip = False # Applies static fringe flat to remove fringes from the data\n", " spec2.photom.skip = False # Applies photometric calibration flat\n", " spec2.cube_build.skip = False # Build quicklook data cubes (needed for Master background in stage 3)\n", " spec2.extract_1d.skip = False # Extract quicklook 1d spectra (needed for Master background in stage 3)\n", " \n", " spec2.save_results = True\n", " spec2(filename)" ] }, { "cell_type": "code", "execution_count": null, "id": "10ff3a2a", "metadata": {}, "outputs": [], "source": [ "# Look for input uncalibrated slope files from the Detector1 pipeline\n", "sstring = det1_dir + 'det*rate.fits'\n", "ratefiles = sorted(glob.glob(sstring))\n", "print('Found ' + str(len(ratefiles)) + ' input files to process')" ] }, { "cell_type": "code", "execution_count": null, "id": "e4ab8d2f", "metadata": {}, "outputs": [], "source": [ "# Simulated data doesn't have the right keywords to tell the pipeline what kind of data is being processed.\n", "# Overwrite rate file header info to specify that these are point sources.\n", "for ii in range(0, len(ratefiles)):\n", " hdu = fits.open(ratefiles[ii])\n", " hdu[1].header['SRCTYPE'] = 'POINT'\n", " hdu.writeto(ratefiles[ii], overwrite=True)\n", " hdu.close()" ] }, { "cell_type": "code", "execution_count": null, "id": "fa5fd149", "metadata": {}, "outputs": [], "source": [ "# If rerunning long pipeline steps, actually run the step\n", "if (redolong is True):\n", " runmany(runspec2, ratefiles)\n", " \n", "# Otherwise, just copy cached outputs into our output directory structure\n", "else:\n", " sstring = os.path.join(spec2_cachedir, 'det*.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": "ec795648", "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": "20ccb089", "metadata": {}, "source": [ "5.-Spec3 Pipeline: Default configuration (12 per-band cubes)\n", "------------------\n", "\n", "
\n", "Here we'll run the Spec3 pipeline in close to 'default' configuration, where we produce a composite data cube from all dithered exposures for each of the 12 MRS bands.\n", "\n", " \n", "Note that we will not be using the 'Master Background' step as we're working with point sources. At the present time background subtraction for point sources is performed during final spectral extraction from the 3d data cubes using an annulus to define the background region. If we were working with extended sources, we would also need to populate our Level3 Association file with the appropriate Spec2 X1d filenames of the dedicated background observations.\n", "\n", "A 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/pipeline/calwebb_spec3.html\n", " \n", "
" ] }, { "cell_type": "code", "execution_count": null, "id": "127897fc", "metadata": {}, "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": "code", "execution_count": null, "id": "bcadac5d", "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": "code", "execution_count": null, "id": "5bef9d09", "metadata": {}, "outputs": [], "source": [ "# Define a convenience function that will split the input _cal.fits files into their corresponding bands\n", "def sort_calfiles(files):\n", " channel = []\n", " band = []\n", "\n", " for file in files:\n", " hdr = (fits.open(file))[0].header\n", " channel.append(hdr['CHANNEL'])\n", " band.append(hdr['BAND'])\n", " channel = np.array(channel)\n", " band = np.array(band)\n", "\n", " indx = np.where((channel == '12') & (band == 'SHORT'))\n", " files12A = files[indx]\n", " indx = np.where((channel == '12') & (band == 'MEDIUM'))\n", " files12B = files[indx]\n", " indx = np.where((channel == '12') & (band == 'LONG'))\n", " files12C = files[indx]\n", " indx = np.where((channel == '34') & (band == 'SHORT'))\n", " files34A = files[indx]\n", " indx = np.where((channel == '34') & (band == 'MEDIUM'))\n", " files34B = files[indx]\n", " indx = np.where((channel == '34') & (band == 'LONG'))\n", " files34C = files[indx]\n", " \n", " return files12A, files12B, files12C, files34A, files34B, files34C" ] }, { "cell_type": "code", "execution_count": null, "id": "6a9aa3d7", "metadata": {}, "outputs": [], "source": [ "# Find and sort all of the input files\n", "sstring = spec2_dir + 'det*cal.fits'\n", "calfiles = np.array(sorted(glob.glob(sstring)))\n", "sortfiles = sort_calfiles(calfiles) # Split them up into bands\n", "print('Found ' + str(len(calfiles)) + ' input files to process')" ] }, { "cell_type": "code", "execution_count": null, "id": "69e17a38", "metadata": {}, "outputs": [], "source": [ "# Create association files for each band with exposures\n", "# Write out the files and add the filenames into our list\n", "\n", "# We could also write all files to a single .json file and cube_build would be smart enough\n", "# to know how to split this up, but then we wouldn't be able to parallelize the runs.\n", "asnlist = []\n", "names = ['12A', '12B', '12C', '34A', '34B', '34C']\n", "for ii in range(0, len(sortfiles)):\n", " thesefiles = sortfiles[ii]\n", " ninband = len(thesefiles)\n", " if (ninband > 0):\n", " filename = 'l3asn-' + names[ii] + '.json'\n", " asnlist.append(filename)\n", " writel3asn(thesefiles, filename, 'Level3')" ] }, { "cell_type": "code", "execution_count": null, "id": "389cc285", "metadata": {}, "outputs": [], "source": [ "# Define a function that will call the spec3 pipeline with our desired set of parameters\n", "def runspec3(filename):\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('l3asn-12A.json')# The exact asn file used doesn't matter\n", " spec3 = Spec3Pipeline.from_config_section(crds_config)\n", " \n", " spec3.output_dir = spec3_dir\n", " spec3.save_results = True\n", " \n", " spec3.master_background.skip = True # Computes and subtracts a master background signal\n", " spec3.outlier_detection.skip = False # Identifies and flags any pixels with values that produce outliers in overlapping regions of cube space\n", " spec3.mrs_imatch.skip = False # Ensure that there are no jumps in the background between individual exposures\n", " spec3.cube_build.skip = False # Build the composite data cubes\n", " spec3.extract_1d.skip = False # Extract 1d spectra from the composite data cubes\n", "\n", " spec3(filename)" ] }, { "cell_type": "code", "execution_count": null, "id": "5653e513", "metadata": {}, "outputs": [], "source": [ "# If rerunning long pipeline steps, actually run the step\n", "if (redolong is True):\n", " runmany(runspec3, asnlist)\n", " \n", "# Otherwise, just copy cached outputs into our output directory structure\n", "else:\n", " sstring = os.path.join(spec3_cachedir, 'Level3*.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": "53f4efda", "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": "670edaa0", "metadata": {}, "source": [ "6.-Spec3 Pipeline: Channel configuration (4 per-channel cubes)\n", "------------------\n", "\n", "
\n", "Previously we used the Spec3 pipeline to make a composite data cube for each of the 12 MRS bands, but the pipeline can also be configured to build cubes combining many bands together into the same data cube. For instance, it can produce 4 per-channel data cubes instead (where, e.g., 1A, 1B, and 1C wavelengths are all combined into a single 'Ch1' cube). Note that this doesn't rely on sampling the individual band cubes- the 2d calibrated data are simply built directly onto this new output grid.\n", "\n", "Let's run Spec3 again, but this time making cubes that cover a full channel (composed of three individual bands).\n", "
" ] }, { "cell_type": "code", "execution_count": null, "id": "fc2e9401", "metadata": {}, "outputs": [], "source": [ "# First we'll make an association file that includes all of the different band exposures\n", "sstring = spec2_dir + 'det*cal.fits'\n", "calfiles = np.array(sorted(glob.glob(sstring)))\n", "writel3asn(calfiles, 'l3asn.json', 'Level3')\n", "print('Found ' + str(len(calfiles)) + ' input files to process')" ] }, { "cell_type": "code", "execution_count": null, "id": "c75c7392", "metadata": {}, "outputs": [], "source": [ "# And we'll define a new Spec3 instance with different cube building parameters\n", "# We list a few of these options here (for more extensive options see\n", "# https://jwst-pipeline.readthedocs.io/en/latest/jwst/cube_build/arguments.html)\n", "# We'll skip master_background, outlier detection, and mrs imatch this time to save runtime\n", "def runspec3_ch(filename):\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('l3asn.json')# The exact asn file used doesn't matter\n", " spec3 = Spec3Pipeline.from_config_section(crds_config)\n", " \n", " spec3.output_dir = spec3_dir\n", " spec3.save_results = True\n", " \n", " # Whether to skip individual pipeline steps or not\n", " spec3.master_background.skip = True # Computes and subtracts a master background signal\n", " spec3.outlier_detection.skip = True # Identifies and flags any pixels with values that produce outliers in overlapping regions of cube space\n", " spec3.mrs_imatch.skip = True # Ensure that there are no jumps in the background between individual exposures\n", " spec3.cube_build.skip = False # Build the composite data cubes\n", " spec3.extract_1d.skip = False # Extract 1d spectra from the composite data cubes\n", " \n", " # Some of the available cube building options\n", " spec3.cube_build.output_file = 'chancube' # Custom output name\n", " spec3.cube_build.output_type = 'channel' # We want to create per-channel cubes\n", " # spec3.cube_build.channel = 'ALL' # Build everything from just channel 1 into a single cube (we could also choose '2','3','4', or 'ALL')\n", " # spec3.cube_build.weighting = 'emsm' # Current default cube build method uses a radial exponential Modified Shepard algorithm, but a 'drizzle' 3d Drizzle option is coming soon\n", " # spec3.cube_build.coord_system = 'ifualign' # The default is 'skyalign'; they differ in whether cube is aligned with North/East or the IFU field angle\n", " # spec3.cube_build.scale1 = 0.5 # Output cube spaxel scale (arcsec) in dimension 1 if setting it by hand\n", " # spec3.cube_build.scale2 = 0.5 # Output cube spaxel scale (arcsec) in dimension 2 if setting it by hand\n", " # spec3.cube_build.scalew = 0.002 # Output cube spaxel size (microns) in dimension 3 if setting it by hand\n", " \n", " spec3(filename)" ] }, { "cell_type": "code", "execution_count": null, "id": "1b2ba992", "metadata": {}, "outputs": [], "source": [ "# If rerunning long pipeline steps, actually run the step\n", "if (redolong is True):\n", " runspec3_ch('l3asn.json')\n", " \n", "# Otherwise, just copy cached outputs into our output directory structure\n", "else:\n", " sstring = os.path.join(spec3_cachedir, 'chancube*.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": "e0a31ad3", "metadata": {}, "source": [ "7.-Spec3 Pipeline: Uber configuration (1 all-wavelength cube)\n", "------------------\n", "\n", "
\n", "Now we can also tell the pipeline to combine **everything** into a single cube spanning the wavelength range from 5-28 microns.\n", "\n", "These cubes are quite unwieldy; they're about 9000 spectral planes in length and nearly 1GB in size each. However, it can be a convenience to have all of the bands together in one place. There are significant drawbacks to these cubes as well though. Since both the PSF and field size change substantially across the MIRI wavelength range the ideal cube spaxel sizes also change. Since cubes can only use a spaxel size in the spatial dimension, a single cube combining both Channel 1 and Channel 4 must therefore either dramatically undersample the data at Ch1 or dramatically oversample the data at Ch4. Likewise, the spectral sampling must change too, although we have more freedom here since the Cube Build algorithm defines a non-linear wavelength solution that changes step size with wavelength in order to more optimally sample spectral features at all wavelengths simultaneously.\n", "
" ] }, { "cell_type": "code", "execution_count": null, "id": "8d035604", "metadata": {}, "outputs": [], "source": [ "# And we'll define a new Spec3 instance with different cube building parameters\n", "# We list a few of these options here (for more extensive options see\n", "# https://jwst-pipeline.readthedocs.io/en/latest/jwst/cube_build/arguments.html)\n", "# We'll skip master_background, outlier detection, and mrs imatch this time to save runtime\n", "def runspec3_all(filename):\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('l3asn-12A.json')# The exact asn file used doesn't matter\n", " spec3 = Spec3Pipeline.from_config_section(crds_config)\n", " \n", " spec3.output_dir = spec3_dir\n", " spec3.save_results = True\n", " \n", " # Whether to skip individual pipeline steps or not\n", " spec3.master_background.skip = True # Computes and subtracts a master background signal\n", " spec3.outlier_detection.skip = True # Identifies and flags any pixels with values that produce outliers in overlapping regions of cube space\n", " spec3.mrs_imatch.skip = True # Ensure that there are no jumps in the background between individual exposures\n", " spec3.cube_build.skip = False # Build the composite data cubes\n", " spec3.extract_1d.skip = False # Extract 1d spectra from the composite data cubes\n", " \n", " # Some of the available cube building options\n", " spec3.cube_build.output_file = 'allcube' # Custom output name\n", " spec3.cube_build.output_type = 'multi'\n", " # spec3.cube_build.weighting = 'emsm' # Current default cube build method uses a radial exponential Modified Shepard algorithm, but a 'drizzle' 3d Drizzle option is coming soon\n", " # spec3.cube_build.coord_system = 'ifualign' # The default is 'skyalign'; they differ in whether cube is aligned with North/East or the IFU field angle\n", " # spec3.cube_build.channel = 'ALL' # Build everything from just channel 1 into a single cube (we could also choose '2','3','4', or 'ALL')\n", " # spec3.cube_build.scale1 = 0.5 # Output cube spaxel scale (arcsec) in dimension 1 if setting it by hand\n", " # spec3.cube_build.scale2 = 0.5 # Output cube spaxel scale (arcsec) in dimension 2 if setting it by hand\n", " # spec3.cube_build.scalew = 0.002 # Output cube spaxel size (microns) in dimension 3 if setting it by hand\n", " \n", " spec3(filename)" ] }, { "cell_type": "code", "execution_count": null, "id": "1608a4e6", "metadata": {}, "outputs": [], "source": [ "# If rerunning long pipeline steps, actually run the step\n", "if (redolong is True):\n", " runspec3_all('l3asn.json')\n", " \n", "# Otherwise, just copy cached outputs into our output directory structure\n", "else:\n", " sstring = os.path.join(spec3_cachedir, 'allcube*.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": "9b1c46c5", "metadata": {}, "source": [ "8.-Examine the output cubes, and compare spectra to the simulation inputs\n", "------------------\n", "\n", "
\n", "Now let's look at the final products (cubes and 1d spectra) created by the pipeline, and compare the different kinds of data cubes.\n", "\n", "As above, remember that the data cubes created by the JWST pipeline are in SURFACE BRIGHTNESS units (MJy/steradian), not flux units.\n", "
" ] }, { "cell_type": "code", "execution_count": null, "id": "fc9babac", "metadata": {}, "outputs": [], "source": [ "# Show an image of the Ch1A, Ch1, and ALL cubes\n", "hdu1A = fits.open(spec3_dir + 'Level3_ch1-short_s3d.fits')\n", "data1A = hdu1A['SCI'].data\n", "hdr1A = hdu1A['SCI'].header\n", "# Linear wavelength solution in per-band cubes\n", "wave1A = np.arange(hdr1A['NAXIS3']) * hdr1A['CDELT3'] + hdr1A['CRVAL3']\n", "\n", "cubefile_ch1 = glob.glob(spec3_dir + 'chancube_ch1*s3d.fits')\n", "hdu1 = fits.open(cubefile_ch1[0])\n", "data1 = hdu1['SCI'].data\n", "hdr1 = hdu1['SCI'].header\n", "# Linear wavelength solution in per-channel cubes\n", "wave1 = np.arange(hdr1['NAXIS3']) * hdr1['CDELT3'] + hdr1['CRVAL3']\n", "\n", "cubefile_all = glob.glob(spec3_dir + 'allcube_*s3d.fits')\n", "hduALL = fits.open(cubefile_all[0])\n", "dataALL = hduALL['SCI'].data\n", "hdrALL = hduALL['SCI'].header\n", "# Reference table of wavelengths for the ALL cube\n", "waveALL = hduALL['WCS-TABLE'].data['wavelength'][0]\n", "\n", "# Use a logarithmic stretch to make sure that\n", "# we can see the actual cube footprint well\n", "norm = ImageNormalize(data1A[0, :, :], stretch=LogStretch(), vmin=-30, vmax=1e4)\n", "\n", "rc('axes', linewidth=2) \n", "fig, (ax1, ax2, ax3) = plt.subplots(1, 3, figsize=(10, 7), dpi=100)\n", "\n", "# And plot the data. Highlight a pixel in the bad column with a red X\n", "ax1.imshow(data1A[0, :, :], cmap='gray', norm=norm, origin='lower')\n", "ax1.set_title('Ch1A Cube')\n", "\n", "ax2.imshow(data1[0, :, :], cmap='gray', norm=norm, origin='lower')\n", "ax2.set_title('Ch1 Cube')\n", "\n", "ax3.imshow(dataALL[0, :, :], cmap='gray', norm=norm, origin='lower')\n", "ax3.set_title('ALL Cube')" ] }, { "cell_type": "markdown", "id": "52c026cd", "metadata": {}, "source": [ "Figure 1: 5 micron plane for data cubes containing Ch1A (4.9 - 5.8 micron), Ch1 (4.9 - 7.5 micron), and all wavelengths (4.9 - 28.3 micron). Note that the ALL-wavelength cube has a larger field of view as it must accommodate the larger Ch4 field in the cube as well." ] }, { "cell_type": "markdown", "id": "a61b976e", "metadata": {}, "source": [ "We can also compare the spectra of these three kinds of cubes to demonstrate their different wavelength coverage." ] }, { "cell_type": "code", "execution_count": null, "id": "01c914db", "metadata": {}, "outputs": [], "source": [ "# Plot a spectrum of the source from the Ch1A, Ch1 and ALL cubes\n", "spec1A = fits.open(spec3_dir + 'Level3_ch1-short_x1d.fits')\n", "specfile_ch1 = glob.glob(spec3_dir + 'chancube_ch1*x1d.fits')\n", "spec1 = fits.open(specfile_ch1[0])\n", "specfile_all = glob.glob(spec3_dir + 'allcube_ch1*x1d.fits')\n", "specALL = fits.open(specfile_all[0])\n", "\n", "rc('axes', linewidth=2) \n", "fig, ax = plt.subplots(1, 1, figsize=(8, 4), dpi=100)\n", "\n", "plt.plot(wave1A, spec1A['EXTRACT1D'].data['FLUX'], label='Band 1A', zorder=2, color='red')\n", "plt.plot(wave1, spec1['EXTRACT1D'].data['FLUX'], label='Ch 1', zorder=1)\n", "plt.plot(waveALL, specALL['EXTRACT1D'].data['FLUX'], label='All Wavelength', zorder=0)\n", "plt.legend()\n", "plt.xlabel('Wavelength (micron)')\n", "plt.ylabel('Fnu (Jy)')\n", "plt.tight_layout()\n", "\n", "spec1A.close()\n", "spec1.close()\n", "specALL.close()" ] }, { "cell_type": "markdown", "id": "ad34c1f4", "metadata": {}, "source": [ "Figure 2: Spectra of point sources extracted from data cubes containing Ch1A (4.9 - 5.8 micron), Ch1 (4.9 - 7.5 micron), and all wavelengths (4.9 - 28.3 micron)." ] }, { "cell_type": "markdown", "id": "2d101461", "metadata": {}, "source": [ "Finally, we can compare the JWST pipeline output spectrum to the input spectrum used to create the mirisim simulations." ] }, { "cell_type": "code", "execution_count": null, "id": "0979572c", "metadata": {}, "outputs": [], "source": [ "# Find the mirisim input spectrum\n", "inputsim = ascii.read(mirisim_dir + 'ngc5728_mirisim.txt')\n", "inputsim['fnu'] /= 1e6 # Mirisim inputs are in units of uJy; convert to Jy to match pipeline outputs\n", "\n", "# Find the 12-band pipeline output 1d spectra\n", "sstring = spec3_dir + 'Level3*x1d.fits'\n", "x1dfiles = np.array(sorted(glob.glob(sstring)))" ] }, { "cell_type": "code", "execution_count": null, "id": "e455243d", "metadata": {}, "outputs": [], "source": [ "# Loop over input spectra reading them into a big array\n", "x1d_wave = []\n", "x1d_flux = []\n", "medwaves = []\n", "for ii in range(0, len(x1dfiles)):\n", " hdu = fits.open(x1dfiles[ii])\n", " specdata = hdu['EXTRACT1D'].data\n", " x1d_wave.append(specdata['WAVELENGTH'])\n", " x1d_flux.append(specdata['FLUX'])\n", " medwaves.append(np.median(specdata['WAVELENGTH']))\n", " hdu.close()\n", "\n", "x1d_wave = np.array(x1d_wave)\n", "x1d_flux = np.array(x1d_flux) \n", "# For convenience, sort according to increasing wavelength\n", "indx = np.argsort(medwaves)\n", "x1d_wave = x1d_wave[indx]\n", "x1d_flux = x1d_flux[indx]\n", "\n", "# Introduce a 10% kludge factor to account for the fact that the mirisim PSF is oversized and thus loses too much flux beyond the aperture radius\n", "x1d_flux *= 1.1" ] }, { "cell_type": "code", "execution_count": null, "id": "99b9ebc5", "metadata": {}, "outputs": [], "source": [ "# Make the comparison plot\n", "rc('axes', linewidth=2) \n", "fig, ax = plt.subplots(figsize=(8, 4), dpi=100)\n", "\n", "plt.plot(inputsim['wave'], inputsim['fnu'], label='Input Spectrum', color='black')\n", "plt.xlabel('Wavelength (micron)')\n", "plt.ylabel('Fnu (Jy)')\n", "plt.tight_layout()\n", "\n", "for ii in range(0, len(x1d_flux)):\n", " if (ii == 0):\n", " plt.plot(x1d_wave[ii], x1d_flux[ii], color='red', label='Pipeline Result')\n", " else:\n", " plt.plot(x1d_wave[ii], x1d_flux[ii], color='red')\n", " \n", "plt.legend()\n" ] }, { "cell_type": "markdown", "id": "e79e5dac", "metadata": {}, "source": [ "Figure 3: Comparison between simulated input spectrum and output spectrum produced by the pipeline. A 10% kludge has been applied to correct for known deficiencies of the simulated PSF compared to the pipeline reference files; deviations at long wavelengths are due to calibration problems in the simulator." ] }, { "cell_type": "code", "execution_count": null, "id": "3cb79b54", "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": "c186839d", "metadata": {}, "source": [ "
" ] }, { "cell_type": "markdown", "id": "6199b024", "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 }