diff --git a/notebooks/NIRSPEC/FSlit/JWPipeNB-NIRSpec-FS.ipynb b/notebooks/NIRSPEC/FSlit/JWPipeNB-NIRSpec-FS.ipynb new file mode 100644 index 0000000..7362ce2 --- /dev/null +++ b/notebooks/NIRSPEC/FSlit/JWPipeNB-NIRSpec-FS.ipynb @@ -0,0 +1,2490 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "590b784c-96e4-4089-907e-ecb4a2652599", + "metadata": {}, + "source": [ + "\"stsci_logo\" " + ] + }, + { + "cell_type": "markdown", + "id": "79f5e34a-377a-44be-8d25-a7d03dfe4aab", + "metadata": {}, + "source": [ + "# NIRSpec FS Pipeline Notebook" + ] + }, + { + "cell_type": "markdown", + "id": "a70d2e41-e4c9-47c0-b693-74b69df76b8e", + "metadata": {}, + "source": [ + "**Authors**: Elena Manjavacas (emanjavacas@stsci.edu), building on the work of Peter Zeidler (zeidler@stsci.edu), Kayli Glidic (kglidic@stsci.edu), and James Muzerolle (muzerol@stsci.edu); NIRSpec branch
\n", + "**Last Updated**: January 6, 2024
\n", + "**Pipeline Version**: 1.16.0 (Build 11.1, Context jwst_1298.pmap)\n", + "\n", + "**Purpose**:
\n", + "This notebook provides a framework for processing generic Near-Infrared Spectrograph (NIRSpec) fixed slit (FS) data through the three stages of the JWST pipeline. It includes how to use associations for multi-exposure observations and how to interact and work with JWST datamodels. Data is assumed to be organized into three folders: science, background, and associations, as specified in the paths set up below. In most cases, editing cells outside the [Configuration](#1.-Configuration) section is unnecessary unless the standard pipeline processing options or plot parameters need to be modified.\n", + "\n", + "**[Data](#3.-Demo-Mode-Setup-(ignore-if-not-using-demo-data))**:
\n", + "This notebook is set up to use observations of HD1808347 A3V standard star (point source) with the G235M grism obtained by Proposal ID (PID) 1128, Observation 6. The demo data will be automatically downloaded in the `demo_mode` unless disabled (i.e., to use local files instead).\n", + "\n", + "**[JWST pipeline version and CRDS context](#Set-CRDS-Context-and-Server)**: 
\n", + "This notebook was written for the calibration pipeline version given above and uses the context associated with this version of the JWST Calibration Pipeline. Information about this an other contexts can be found in the JWST Calibration Reference Data System (CRDS) \n", + "[server](https://jwst-crds.stsci.edu/). If you use different pipeline versions, please refer to the table [here](https://jwst-crds.stsci.edu/display_build_contexts/) to determine what context to use. To learn more about the differences for the pipeline, read the relevant [documentation](https://jwst-docs.stsci.edu/jwst-science-calibration-pipeline/jwst-operations-pipeline-build-information#references)\n", + "\n", + "Please note that pipeline software development is a continuous process, so results in some cases may be slightly different if a subsequent version is used. **For optimal results, users are strongly encouraged to reprocess their data using the most recent pipeline version and [associated CRDS context](https://jwst-crds.stsci.edu/display_build_contexts/), taking advantage of bug fixes and algorithm improvements.**\n", + "Any [known issues](https://jwst-docs.stsci.edu/known-issues-with-jwst-data/nirspec-known-issues/nirspec-fs-known-issues#gsc.tab=0) for this build are noted in the notebook. \n", + "\n", + "**Updates**:
\n", + "This notebook is regularly updated to incorporate the latest pipeline improvements. Find the most up-to-date version of this notebook [here](https://github.com/spacetelescope/jwst-pipeline-notebooks/). \n", + "\n", + "**Recent Changes**:
\n", + "* October 15, 2024: Converted notebook to follow standard template (kglidic@stsci.edu).
\n", + "* November 4, 2024: Notebook updated to JWST pipeline version 1.16.0 (Build 11.1).\n", + "* January 6, 2024: Updated formatting and added examples for creating association files." + ] + }, + { + "cell_type": "markdown", + "id": "15c418e0-437a-48bf-822a-a2262965666e", + "metadata": { + "jp-MarkdownHeadingCollapsed": true + }, + "source": [ + "## Table of Contents\n", + "\n", + "* [1. Configuration](#1.-Configuration)\n", + "* [2. Package Imports](#2.-Package-Imports)\n", + "* [3. Demo Mode Setup](#3.-Demo-Mode-Setup-(ignore-if-not-using-demo-data))\n", + "* [4. Directory Setup](#4.-Directory-Setup)\n", + "* [5. Stage 1: `Detector1Pipeline` (`calwebb_detector1`)](#5.-Stage-1:-Detector1Pipeline-(calwebb_detector1))\n", + " * [5.1 Configure `Detector1Pipeline`](#5.1-Configure-Detector1Pipeline)\n", + " * [5.2 Run `Detector1Pipeline`](#5.2-Run-Detector1Pipeline)\n", + "* [6. Stage 2: `Spec2Pipeline` (`calwebb_spec2`)](#6.-Stage-2:-Spec2Pipeline-(calwebb_spec2))\n", + " * [6.1 Configure `Spec2Pipeline`](#6.1-Configure-Spec2Pipeline)\n", + " * [6.2 Create `Spec2Pipeline` Association Files](#6.2-Create-Spec2Pipeline-Association-Files)\n", + " * [6.3 Run `Spec2Pipeline`](#6.3-Run-Spec2Pipeline)\n", + "* [7. Stage 3: `Spec3Pipeline` (`calwebb_spec3`)](#7.-Stage-3:-Spec3Pipeline-(calwebb_spec3))\n", + " * [7.1 Configure `Spec3Pipeline`](#7.1-Configure-Spec3Pipeline)\n", + " * [7.2 Create `Spec3Pipeline` Association Files](#7.2-Create-Spec3Pipeline-Association-Files)\n", + " * [7.3 Run `Spec3Pipeline`](#7.3-Run-Spec3Pipeline)\n", + "* [8. Visualize the Data](#8.-Visualize-the-Data)\n", + " * [8.1 Display `Detector1Pipeline` Products](#8.1-Display-Detector1Pipeline-Products)\n", + " * [8.2 Display `Spec2Pipeline` Products](#8.2-Display-Spec2Pipeline-Products)\n", + " * [8.3 Display `Spec3Pipeline` Products](#8.3-Display-Spec3Pipeline-Products)\n", + "* [9. Modifying the EXTRACT1D Reference File (as needed)](#9.-Modifying-the-EXTRACT1D-Reference-File-(as-needed))\n", + "---" + ] + }, + { + "cell_type": "markdown", + "id": "2cc95671-0291-4c81-a253-87c1fb970ad1", + "metadata": {}, + "source": [ + "## 1. Configuration\n", + "\n", + "### Install dependencies and parameters\n", + "\n", + "To make sure that the pipeline version is compatible with the steps discussed below and that the required dependencies and packages get installed, you can create a fresh conda environment and install the provided requirements.txt file before starting this notebook:\n", + "\n", + " conda create -n nirspec_fs_pipeline python=3.12\n", + " conda activate nirspec_fs_pipeline\n", + " pip install -r requirements.txt\n", + "\n", + "### Set the basic parameters to configure the notebook\n", + "\n", + "These parameters determine what data gets used, where data is located (if already on disk), and the type of background subtraction (if any). The list of parameters includes:\n", + "\n", + "* `demo_mode`:\n", + " * `True`: Downloads example data from the [Barbara A. Mikulski Archive for Space Telescopes (MAST)](https://archive.stsci.edu/) and processes it through the pipeline. All processing will occur in a local directory unless modified in [Section 3](#3.-Demo-Mode-Setup-(ignore-if-not-using-demo-data)) below.\n", + " * `False`: Process your own downloaded data; provide its location.

\n", + " \n", + "* **Directories with data**:\n", + " * `sci_dir`: Directory where science observation data is stored.\n", + " * `bg_dir`: Directory where dedicated background observation data is stored.

\n", + " \n", + "* **[Backgroud subtraction methods](https://jwst-pipeline.readthedocs.io/en/latest/jwst/background_subtraction/main.html#spectroscopic-modes:~:text=the%20calwebb_image3%20pipeline.-,Spectroscopic%20Modes,%EF%83%81,-Spectroscopic%20observations%20allow)**:
\n", + " * `master_bg` = `True`: Apply master-background subtraction in `calwebb_spec3`. For dedicated background observations.\n", + " * `pixel_bg` = `True`: Apply pixel-to-pixel background subtraction in `calwebb_spec2`. This is the default pipeline setting. Uses noded observations.

" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "25f22f8a-300e-4cb5-b5fc-d9505cfe0bdb", + "metadata": {}, + "outputs": [], + "source": [ + "# Basic import necessary for configuration.\n", + "\n", + "import os\n", + "import crds\n", + "import warnings\n", + "\n", + "# Control logging level: INFO, WARNING, ERROR.\n", + "# Set up if want to hide RuntimeWarning messages.\n", + "#import logging\n", + "#logging.disable(logging.ERROR)\n", + "warnings.simplefilter(\"ignore\", RuntimeWarning)" + ] + }, + { + "cell_type": "markdown", + "id": "91e0b06e-a03b-42e6-81b0-125a39ca62e1", + "metadata": {}, + "source": [ + "
\n", + " \n", + "Note that `demo_mode` must be set appropriately below.\n", + "
" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "0d08eac8-a2e1-4237-92eb-5c7eba9f115e", + "metadata": {}, + "outputs": [], + "source": [ + "# Set parameters for demo_mode, data mode directories, and processing steps.\n", + "\n", + "# -------------------------------DEMO MODE-----------------------------------\n", + "demo_mode = True\n", + "\n", + "if demo_mode:\n", + " print('Running in demonstration mode using online example data!')\n", + "\n", + "# ----------------------------User Mode Directories--------------------------\n", + "else: # If demo_mode = False, look for user data in these paths.\n", + "\n", + " # Set directory paths for processing specific data; adjust to your local\n", + " # directory setup (examples provided below).\n", + " basedir = os.path.join(os.getcwd(), '')\n", + "\n", + " # Directory to science observation data; expects uncalibrated data in\n", + " # sci_dir/uncal/ and results in stage1, stage2, and stage3 directories.\n", + " sci_dir = os.path.join(basedir, 'fs_data_01128/Obs006', '')\n", + "\n", + " # Directory to background observation data; expects uncalibrated data in\n", + " # bg_dir/uncal/ and results in stage1, stage2, and stage3 directories.\n", + " #bg_dir = os.path.join(basedir, 'fs_data_02288/Obs002', '')\n", + " bg_dir = '' # If no background observation, use an empty string.\n", + "\n", + " # Directory to stage 2/3 association files.\n", + " asn_dir = os.path.join(basedir, 'asn', '')\n", + "\n", + "# ---------------------------Set Processing Steps----------------------------\n", + "# Individual pipeline stages can be turned on/off here. Note that a later\n", + "# stage won't be able to run unless data products have already been\n", + "# produced from the prior stage.\n", + "\n", + "# Science processing.\n", + "dodet1 = True # calwebb_detector1\n", + "dospec2 = True # calwebb_image2\n", + "dospec3 = True # calwebb_image3\n", + "doviz = True # Visualize calwebb outputs\n", + "\n", + "# How should background subtraction be done?\n", + "# If none are selected, data will not be background subtracted.\n", + "# Note: Master-background subtraction is for dedicated background observations.\n", + "# Dedicated backgrounds must be processed through spec2 first.\n", + "master_bg = False # Master-background subtraction in spec3.\n", + "pixel_bg = True # Pixel-based background subtraction in spec2." + ] + }, + { + "cell_type": "markdown", + "id": "66489df1-6a4e-427f-adf9-59813012f391", + "metadata": {}, + "source": [ + "---\n", + "\n", + "### Set CRDS Context and Server\n", + "\n", + "Before importing `CRDS` and `JWST` modules, we need to configure our environment. This includes defining a CRDS cache directory in which to keep the reference files that will be used by the calibration pipeline. If the local CRDS cache directory has not been set, it will automatically be created in the home directory." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "ef3d8a2a-9da6-401b-a4d5-825cf1d41d98", + "metadata": {}, + "outputs": [], + "source": [ + "# ------------------------Set CRDS context and paths------------------------\n", + "# Each version of the calibration pipeline is associated with a specific CRDS\n", + "# context file. The pipeline will select the appropriate context file behind\n", + "# the scenes while running. However, if you wish to override the default context\n", + "# file and run the pipeline with a different context, you can set that using\n", + "# the CRDS_CONTEXT environment variable. Here, we show how this is done,\n", + "# although we leave the line commented out in order to use the default context.\n", + "# If you wish to specify a different context, uncomment the line below.\n", + "# os.environ['CRDS_CONTEXT'] = 'jwst_1298.pmap' # CRDS context for 1.16.0\n", + "\n", + "# Set CRDS cache directory to user home dir if not already set.\n", + "if os.getenv('CRDS_PATH') is None:\n", + " os.environ['CRDS_PATH'] = os.path.join(os.path.expanduser('~'), 'crds_cache')\n", + "\n", + "# Check whether the CRDS server URL has been set. If not, set it.\n", + "if os.getenv('CRDS_SERVER_URL') is None:\n", + " os.environ['CRDS_SERVER_URL'] = 'https://jwst-crds.stsci.edu'\n", + "\n", + "# Output the current CRDS path and server URL in use.\n", + "print(f\"CRDS local filepath: {os.environ.get('CRDS_PATH', 'Not Set')}\")\n", + "print(f\"CRDS file server: {os.environ.get('CRDS_SERVER_URL', 'Not Set')}\")\n", + "print(f\"Default CRDS Context for JWST Build: \"\n", + " f\"{crds.get_default_context('jwst', state='build')}\")\n", + "print(f\"Currently Active CRDS Context: {os.environ.get('CRDS_CONTEXT', 'Not Set')}\")" + ] + }, + { + "cell_type": "markdown", + "id": "a5cc8ea9-a9a6-456c-8377-d9294fbff5b0", + "metadata": {}, + "source": [ + "\n", + "---\n", + "\n", + "## 2. Package Imports\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "9efc0021-ddc4-44f1-a4ca-42f429a75603", + "metadata": {}, + "outputs": [], + "source": [ + "# Use the entire available screen width for this notebook.\n", + "from IPython.display import display, HTML, JSON\n", + "display(HTML(\"\"))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "44a12eb7-b697-4cb3-9c8e-df7f8a936fc6", + "metadata": {}, + "outputs": [], + "source": [ + "# ----------------------General Imports----------------------\n", + "import time\n", + "import glob\n", + "import json\n", + "import itertools\n", + "import numpy as np\n", + "\n", + "# ----------------------Astropy Imports----------------------\n", + "# Astropy utilities for opening FITS files, downloading demo files, etc.\n", + "from astropy.io import fits\n", + "from astropy.stats import sigma_clip\n", + "from astropy.visualization import ImageNormalize, ManualInterval, LogStretch\n", + "from astropy.visualization import LinearStretch, AsinhStretch, simple_norm\n", + "\n", + "# ----------------------Astroquery Import----------------------\n", + "from astroquery.mast import Observations\n", + "\n", + "# ----------------------Plotting Imports---------------------\n", + "import matplotlib.pyplot as plt\n", + "from matplotlib.patches import Rectangle\n", + "from matplotlib.collections import PatchCollection" + ] + }, + { + "cell_type": "markdown", + "id": "694021fe-c418-427e-8bf8-93624565968c", + "metadata": {}, + "source": [ + "
\n", + "\n", + "Installation instructions for the JWST pipeline found in [JWST User Documentation](https://jwst-docs.stsci.edu/jwst-science-calibration-pipeline-overview) and\n", + "[ReadtheDocs](https://jwst-pipeline.readthedocs.io) • \n", + "[Github](https://github.com/spacetelescope/jwst)\n", + "\n", + "
" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "e5f34a5e-7372-4468-a89a-35e51f833adc", + "metadata": {}, + "outputs": [], + "source": [ + "# ----------------------JWST Calibration Pipeline Imports----------------------\n", + "import jwst # Import the base JWST and CRDS packages.\n", + "from crds.client import api\n", + "from stpipe import crds_client\n", + "\n", + "# JWST pipelines (each encompassing many steps).\n", + "from jwst.pipeline import Detector1Pipeline # calwebb_detector1\n", + "from jwst.pipeline import Spec2Pipeline # calwebb_spec2\n", + "from jwst.pipeline import Spec3Pipeline # calwebb_spec3\n", + "from jwst.extract_1d import Extract1dStep # Extract1D Step\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 ASN files.\n", + "from jwst.associations.lib.rules_level2_base import DMSLevel2bBase # Lvl2 ASN file.\n", + "from jwst.associations.lib.rules_level3_base import DMS_Level3_Base # Lvl3 ASN file.\n", + "\n", + "print(\"JWST Calibration Pipeline Version = {}\".format(jwst.__version__))" + ] + }, + { + "cell_type": "markdown", + "id": "b1ad4b68-ef9c-4879-98cb-8d5a74861fd5", + "metadata": {}, + "source": [ + "---\n", + "### Define Convience Functions\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "02c2b6eb-9458-438b-9c20-08fb34e4d6a8", + "metadata": {}, + "outputs": [], + "source": [ + "def update_asn_paths(asn_files,\n", + " exclude_dirs=[]):\n", + " \"\"\"\n", + " Update the expname field in an ASN file with its absolute path.\n", + "\n", + " The absolute path is determined by locating matching files in the\n", + " current directory while excluding any files found in the specified\n", + " excluded directories. Absolute paths ensure that the pipeline correctly\n", + " locates the files, regardless of the ASN file's location.\n", + "\n", + " Parameters\n", + " ----------\n", + " asn_files : list of str\n", + " List of ASN files to update.\n", + " exclude_dirs : list of str, optional\n", + " List of directories to exclude in the search.\n", + "\n", + " Returns\n", + " -------\n", + " None.\n", + " \"\"\"\n", + "\n", + " # Ensure asn_files is a list.\n", + " asn_files = [asn_files] if isinstance(asn_files, str) else asn_files\n", + "\n", + " for asn in asn_files:\n", + " with open(asn, 'r') as file:\n", + " data = json.load(file)\n", + "\n", + " update = False\n", + "\n", + " # Loop through each product and its members.\n", + " for product in data['products']:\n", + " for member in product['members']:\n", + " search_pattern = f\"**/{member['expname']}\"\n", + " filtered_files = [f for f in glob.glob(os.path.join(os.getcwd(),\n", + " search_pattern), recursive=True) if not\n", + " any(f.startswith(exc) for exc in exclude_dirs)]\n", + " if len(filtered_files) > 0:\n", + " member['expname'] = filtered_files[0]\n", + " update = True\n", + "\n", + " basn = os.path.basename(asn)\n", + " if update:\n", + " try:\n", + " with open(asn, 'w') as json_file:\n", + " json.dump(data, json_file, indent=4)\n", + " print(f\"{basn} 'expname' paths have been updated to absolute paths!\")\n", + " except Exception as e:\n", + " print(f\"Error saving updated file {asn}: {e}\")\n", + " else:\n", + " print(f\"{basn} 'expname' paths NOT updated! Paths may be absolute.\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "bba3a01a-41c7-4f47-990b-14cb04f38336", + "metadata": {}, + "outputs": [], + "source": [ + "def write_asn(sci_files,\n", + " bg_files=None,\n", + " leakcal_files=None,\n", + " asn_output_dir='.',\n", + " level=2):\n", + " \"\"\"\n", + " Write a JWST Level 2b or Level 3 association (ASN) file for the given science\n", + " files, with optional background or imprint files (applicable for IFU/MOS Spec2 ASN).\n", + "\n", + " Parameters\n", + " ----------\n", + " sci_files : str or list\n", + " Path to the science file or a list of paths to science files.\n", + " bg_files : str or list, optional\n", + " Path to a background file or a list of paths to background files.\n", + " If None, no background files are included.\n", + " leakcal_files : str or list, optional\n", + " Path to an imprint file or a list of paths to imprint files.\n", + " If None, no imprint files are included.\n", + " asn_output_dir : str\n", + " Output path for the ASN JSON file.\n", + " level : int, optional\n", + " ASN level (2 for Spec2, 3 for Spec3).\n", + "\n", + " Returns\n", + " -------\n", + " None.\n", + " \"\"\"\n", + "\n", + " # Ensure inputs are lists.\n", + " sci_files = ([sci_files] if isinstance(sci_files, str) else sci_files)\n", + " bg_files = ([bg_files] if isinstance(bg_files, str) else bg_files or [])\n", + " leakcal_files = ([leakcal_files] if isinstance(leakcal_files, str)\n", + " else leakcal_files or [])\n", + "\n", + " # Header information.\n", + " prgm = fits.getval(sci_files[0], 'PROGRAM')\n", + " obs = fits.getval(sci_files[0], 'OBSERVTN')\n", + " filt = fits.getval(sci_files[0], 'FILTER').lower()\n", + " grating = fits.getval(sci_files[0], 'GRATING').lower()\n", + " name = f'jw{prgm}-o{obs}_{{source_id}}_nirspec_{filt}-{grating}-{{slit_name}}'\n", + "\n", + " if level == 3:\n", + " # Create Level 3 ASN file.\n", + " asn = afl.asn_from_list([os.path.abspath(f) for f in sci_files],\n", + " rule=DMS_Level3_Base,\n", + " product_name=name)\n", + "\n", + " # Add background files to the association.\n", + " for bg_file in bg_files:\n", + " asn['products'][0]['members'].append({\n", + " 'expname': os.path.abspath(bg_file), 'exptype': 'background'\n", + " })\n", + "\n", + " primary_prefix = \"_\".join(os.path.basename(sci_files[0]).split(\"_\")[:1])\n", + " asn_name = primary_prefix+\"_manual_spec3_asn.json\"\n", + "\n", + " # Serialize and write the ASN to a JSON file.\n", + " _, serialized = asn.dump()\n", + " asn_output = os.path.join(asn_output_dir, asn_name)\n", + " with open(asn_output, 'w') as outfile:\n", + " outfile.write(serialized)\n", + " print(f\"Saved Association File: {asn_output}\")\n", + "\n", + " elif level == 2:\n", + " # Create Level 2 ASN files.\n", + " # Loop over each science file, create ASN file for each.\n", + " # If nodded, rotate the primary science exposure for\n", + " # pixel-to-pixel background subtraction.\n", + " for i, primary_sci_file in enumerate(sci_files):\n", + " # Create the basic association for the science files.\n", + " product_name = os.path.basename(primary_sci_file).rsplit('_', 1)[0]\n", + " asn = afl.asn_from_list([os.path.abspath(primary_sci_file)],\n", + " rule=DMSLevel2bBase,\n", + " product_name=product_name)\n", + "\n", + " # Detector NRS1 or NRS2.\n", + " det = fits.getval(primary_sci_file, 'DETECTOR')\n", + "\n", + " # Set up ASN file for pixel-to-pixel background subtraction if nodded.\n", + " try:\n", + " NOD_TYPE = fits.getval(primary_sci_file, 'NOD_TYPE')\n", + " except KeyError:\n", + " NOD_TYPE = None # Return None if the keyword does not exist.\n", + " PATTTYPE = fits.getval(primary_sci_file, 'PATTTYPE')\n", + " if (PATTTYPE and 'NOD' in PATTTYPE) or NOD_TYPE:\n", + " # Add background files for nodded exposures by matching the prefix.\n", + " nod_backgrounds = [\n", + " f for f in sci_files if f != primary_sci_file and det.lower() in f\n", + " ]\n", + " asn['products'][0]['members'].extend(\n", + " {'expname': os.path.abspath(nod_bg), 'exptype': 'background'}\n", + " for nod_bg in nod_backgrounds\n", + " )\n", + "\n", + " # If background files are provided, add those.\n", + " if bg_files:\n", + " asn['products'][0]['members'].extend(\n", + " {'expname': os.path.abspath(bg_file), 'exptype': 'background'}\n", + " for bg_file in bg_files\n", + " )\n", + "\n", + " # If leakcal files are provided, add those.\n", + " if leakcal_files:\n", + " asn['products'][0]['members'].extend(\n", + " {'expname': os.path.abspath(imp_file), 'exptype': 'imprint'}\n", + " for imp_file in leakcal_files\n", + " )\n", + "\n", + " asn_name = f\"jw{prgm}-o{obs}_000{i+1}_{det}_manual_spec2_asn.json\"\n", + "\n", + " # Serialize and write the ASN to a JSON file.\n", + " _, serialized = asn.dump()\n", + " asn_output = os.path.join(asn_output_dir, asn_name)\n", + " with open(asn_output, 'w') as outfile:\n", + " outfile.write(serialized)\n", + " print(f\"Saved Association File: {asn_output}\")\n", + " else:\n", + " raise ValueError(\"Invalid level for ASN function.\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "dfc00e09-1468-4e1f-8e4e-7274ad05364a", + "metadata": {}, + "outputs": [], + "source": [ + "# Start a timer to keep track of runtime.\n", + "time0 = time.perf_counter()" + ] + }, + { + "cell_type": "markdown", + "id": "060a63f8-3183-4408-99e5-461564c93ebc", + "metadata": {}, + "source": [ + "---\n", + "\n", + "## 3. Demo Mode Setup (ignore if not using demo data)\n", + "\n", + "
\n", + "\n", + " The data in this notebook is public and does not require a token. For other data sets, you may need to provide a token. For more infomation visit [astroquery](https://astroquery.readthedocs.io/en/latest/index.html#) documentation.\n", + " \n", + "
\n", + "\n", + "If running in demonstration mode, set up the program information to retrieve the uncalibrated data (`_uncal.fits`) automatically from MAST using `astroquery`. MAST provides flexibility by allowing searches based on proposal ID and observation ID, rather than relying solely on filenames. More information about the JWST file naming conventions can be found [here](https://jwst-pipeline.readthedocs.io/en/latest/jwst/data_products/file_naming.html).\n", + "\n", + "\n", + "The FS demo data in this notebook is from the [NIRSpec calibration program 1128](https://www.stsci.edu/jwst/science-execution/program-information?id=1128) and features observations of HD1808347 (point source) using the G235M grism. The program setup is briefly summarized in the table below.\n", + "\n", + "\n", + "| Demo Target: HD1808347 A3V Standard Star | | | \n", + "|:-----------:|:-------:|:---:|\n", + "| PROGRAM | 01128 | Program number | \n", + "| OBSERVTN | 006 | Observation number | \n", + "| [GRATING/FILTER](https://jwst-docs.stsci.edu/jwst-near-infrared-spectrograph/nirspec-observing-modes/nirspec-fixed-slits-spectroscopy#gsc.tab=0:~:text=%C2%A0Table%202.%20NIRSpec%20A%20fixed%20slit%C2%A0instrument%20configurations%2C%20resolutions%2C%20and%20wavelength%20ranges) | G235M/F170LP | λ: 1.66–3.17 μm (a medium resolution, R ~ 1000) |\n", + "| SUBARRAY | SUBS200A1 | Subarray used | \n", + "| NINTS | 2 | Number of integrations in exposure | \n", + "| NGROUPS | 30 | Number of groups in integration |\n", + "| DURATION | 96.637 [s] | Total duration of one exposure | \n", + "| READPATT | NRSRAPID | Readout pattern | \n", + "| PATTTYPE | 3-POINT-NOD | Primary dither pattern type | \n", + "| NUMDTHPT | 3 | Total number of points in pattern | \n", + "| SRCTYAPT | POINT | Source type selected in APT |\n", + "\n", + "> **Note:** The presence of a physical gap between detectors affects high-resolution FS observations because the spectra are long enough to span both NIRSpec detectors. [More Info ...](https://jwst-docs.stsci.edu/jwst-near-infrared-spectrograph/nirspec-operations/nirspec-fs-operations/nirspec-fs-wavelength-ranges-and-gaps#gsc.tab=0)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "095668db-97a9-4403-a17b-c45f41ad611f", + "metadata": {}, + "outputs": [], + "source": [ + "# Set up the program information and directories to collect\n", + "# the data in demo_mode.\n", + "if demo_mode:\n", + "\n", + " print('Running in demonstration mode. '\n", + " 'Example data will be downloaded from MAST!')\n", + "\n", + " # For non public data sets, you may need to provide a token.\n", + " # However, for security it is not recommended to enter tokens into\n", + " # a terminal or Jupyter notebook.\n", + " # Observations.login(token=\"your-token\")\n", + "\n", + " # --------------Program and observation information--------------\n", + " program = \"01128\"\n", + " sci_observtn = \"006\"\n", + " bg_observtn = None\n", + " filters = [\"F170LP;G235M\"]\n", + "\n", + " # ----------Define the base and observation directories----------\n", + " basedir = os.path.join('.', f'fs_data_{program}')\n", + " sci_dir = os.path.join(basedir, f'Obs{sci_observtn}')\n", + " asn_dir = os.path.join(basedir, 'asn/')\n", + " uncal_dir = os.path.join(sci_dir, 'uncal/')\n", + "\n", + " # If no background observation, leave blank.\n", + " bg_dir = os.path.join(basedir, f'Obs{bg_observtn}') if bg_observtn else ''\n", + " uncal_bgdir = os.path.join(bg_dir, 'uncal/') if bg_observtn else ''\n", + "\n", + " # ------Ensure directories for downloading MAST data exists------\n", + " os.makedirs(uncal_dir, exist_ok=True)\n", + " os.makedirs(asn_dir, exist_ok=True)\n", + " # Makes directory only when a background observation is provided.\n", + " if bg_observtn:\n", + " os.makedirs(uncal_bgdir, exist_ok=True)" + ] + }, + { + "cell_type": "markdown", + "id": "e4aacb35-f53b-4095-882c-bd66e11403cd", + "metadata": {}, + "source": [ + "Identify a list of science (SCI) and background (BG) uncalibrated files associated with visits. Note that this `demo_mode` case doesn't have background observations." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "ba4f6291-943a-425a-aafc-d1232ce0956f", + "metadata": {}, + "outputs": [], + "source": [ + "# Obtain a list of observation IDs for the specified demo program.\n", + "if demo_mode:\n", + "\n", + " # --------------------SCIENCE Observation--------------------\n", + " sci_obs_id_table = Observations.query_criteria(instrument_name=['NIRSPEC/SLIT'],\n", + " provenance_name=[\"CALJWST\"],\n", + " obs_id=[f'*{program}*{sci_observtn}*'])\n", + "\n", + " # ------------------BACKGROUND Observation-------------------\n", + " bg_obs_id_table = Observations.query_criteria(instrument_name=['NIRSPEC/SLIT'],\n", + " provenance_name=[\"CALJWST\"],\n", + " obs_id=[f'*{program}*{bg_observtn}*'])" + ] + }, + { + "cell_type": "markdown", + "id": "9b4ece3c-4200-40cd-928e-1e2e7af20ae8", + "metadata": { + "jp-MarkdownHeadingCollapsed": true + }, + "source": [ + "\n", + "\n", + "The demo dataset consists of six `_uncal.fits` files, each approximately 15 MB in size." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "7e420e58-e282-4345-bbf1-37ad6013f909", + "metadata": {}, + "outputs": [], + "source": [ + "# Convert visits into a list of uncalibrated data and ASN files.\n", + "\n", + "if demo_mode:\n", + "\n", + " file_criteria = {'filters': filters, 'calib_level': [1],\n", + " 'productSubGroupDescription': 'UNCAL'}\n", + " asn_criteria = {'filters': filters, 'calib_level': [2, 3],\n", + " 'productSubGroupDescription': 'ASN'}\n", + "\n", + " # Initialize lists for science, background, and ASN files.\n", + " sci_downloads, bg_downloads, asn_downloads = [], [], []\n", + "\n", + " pfilter = Observations.filter_products # Alias for filter_products method.\n", + "\n", + " # ----------Identify uncalibrated SCIENCE files associated with each visit----------\n", + " for exposure in sci_obs_id_table:\n", + " sci_products = Observations.get_product_list(exposure)\n", + " asn_downloads.extend(pfilter(sci_products, **asn_criteria)['dataURI'])\n", + "\n", + " # Filter for full-size science files (exclude smaller confirmation images).\n", + " avg_sci_size = np.nanmean(sci_products['size'])\n", + " sci_products = sci_products[sci_products['size'] > avg_sci_size]\n", + " sci_downloads.extend(pfilter(sci_products, **file_criteria)['dataURI'])\n", + "\n", + " # --------Identify uncalibrated BACKGROUND files associated with each visit---------\n", + " for exposure in bg_obs_id_table:\n", + " bg_products = Observations.get_product_list(exposure)\n", + "\n", + " # Filter for full-size background files (exclude smaller confirmation images).\n", + " avg_bg_size = np.nanmean(bg_products['size'])\n", + " bg_products = sci_products[bg_products['size'] > avg_bg_size]\n", + " bg_downloads.extend(pfilter(bg_products, **file_criteria)['dataURI'])" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "e37a48cb-d189-4f94-93ba-f92bcfbfcf66", + "metadata": {}, + "outputs": [], + "source": [ + "if demo_mode:\n", + " # Filter out other observations and remove duplicates.\n", + " sci_downloads = {f for f in sci_downloads if f\"jw{program}{sci_observtn}\" in f}\n", + " bg_downloads = {f for f in bg_downloads if f\"jw{program}{bg_observtn}\" in f}\n", + "\n", + " # Selects observation and candidate associations.\n", + " asn_downloads = {\n", + " f for f in asn_downloads\n", + " if any(f\"-{p}{sci_observtn}_\" in f for p in [\"o\", \"c?\"])\n", + " }\n", + "\n", + " print(f\"Science files selected for downloading: {len(sci_downloads)}\")\n", + " print(f\"Background files selected for downloading: {len(bg_downloads)}\")\n", + " print(f\"ASN files selected for downloading: {len(asn_downloads)}\")" + ] + }, + { + "cell_type": "markdown", + "id": "f244c138-ec1f-4abb-a0fd-3b00202a4c34", + "metadata": {}, + "source": [ + "Download the data.\n", + "
\n", + "\n", + "**Warning**: If this notebook is halted during this step, the downloaded file may be incomplete, and cause crashes later on!\n", + "
" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "0f642e1d-71d8-4f10-922e-01a5c6d02ee6", + "metadata": { + "scrolled": true + }, + "outputs": [], + "source": [ + "# Download data and place them into the appropriate directories.\n", + "if demo_mode:\n", + " for file in sci_downloads:\n", + " sci_manifest = Observations.download_file(file, local_path=uncal_dir)\n", + " for file in bg_downloads:\n", + " bg_manifest = Observations.download_file(file, local_path=uncal_bgdir)\n", + " for file in asn_downloads:\n", + " asn_manifest = Observations.download_file(file, local_path=asn_dir)" + ] + }, + { + "cell_type": "markdown", + "id": "23ecf36a-cbe2-433c-a056-87bc72fcd76f", + "metadata": {}, + "source": [ + "---\n", + "\n", + "## 4. Directory Setup\n", + "Set up detailed paths to input/output stages here." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "685ca7d0-fc74-4774-bbf0-a26a2aacefbc", + "metadata": {}, + "outputs": [], + "source": [ + "# Define/create output subdirectories to keep data products organized.\n", + "\n", + "# -----------------------------Science Directories------------------------------\n", + "uncal_dir = os.path.join(sci_dir, 'uncal/') # Uncalibrated pipeline inputs.\n", + "det1_dir = os.path.join(sci_dir, 'stage1/') # calwebb_detector1 pipeline outputs.\n", + "spec2_dir = os.path.join(sci_dir, 'stage2/') # calwebb_spec2 pipeline outputs.\n", + "spec3_dir = os.path.join(sci_dir, 'stage3/') # calwebb_spec3 pipeline outputs.\n", + "\n", + "# Creates the directories if target directory does not exist.\n", + "os.makedirs(det1_dir, exist_ok=True)\n", + "os.makedirs(spec2_dir, exist_ok=True)\n", + "os.makedirs(spec3_dir, exist_ok=True)\n", + "\n", + "# ---------------------------Background Directories-----------------------------\n", + "uncal_bgdir = os.path.join(bg_dir, 'uncal/') # Uncalibrated pipeline inputs.\n", + "det1_bgdir = os.path.join(bg_dir, 'stage1/') # calwebb_detector1 pipeline outputs.\n", + "spec2_bgdir = os.path.join(bg_dir, 'stage2/') # calwebb_spec2 pipeline outputs.\n", + "\n", + "# Creates directories if background observations are provided and do not already exist.\n", + "if bg_dir:\n", + " os.makedirs(det1_bgdir, exist_ok=True)\n", + " os.makedirs(spec2_bgdir, exist_ok=True)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "050f53e4-f225-48b2-bdfb-c6220628a80b", + "metadata": {}, + "outputs": [], + "source": [ + "# Print out the time benchmark.\n", + "time1 = time.perf_counter()\n", + "print(f\"Runtime so far: {round((time1-time0)/60.0, 1):0.4f} min\")" + ] + }, + { + "cell_type": "markdown", + "id": "d41b71ee-f4aa-45ef-865b-9ba83dd82186", + "metadata": {}, + "source": [ + "---\n", + "\n", + "## 5. Stage 1: `Detector1Pipeline` (`calwebb_detector1`)\n", + "\n", + "In this section, we process the data through the `calwebb_detector1` pipeline to create Stage 1 data products.\n", + "\n", + "* **Input**: Raw exposure (`_uncal.fits`) containing original data from all detector readouts (ncols x nrows x ngroups x nintegrations).\n", + "* **Output**: Uncalibrated countrate (slope) image in units of DN/s:\n", + " * `_rate.fits`: A single countrate image averaged over multiple integrations (if available).\n", + " * `_rateints.fits`: Countrate images for each integration, saved in multiple extensions.\n", + "\n", + "The `Detector1Pipeline` applies basic detector-level corrections on a group-by-group basis, followed by ramp fitting for all exposure types, commonly referred to as \"ramps-to-slopes\" processing. \n", + "\n", + "---" + ] + }, + { + "cell_type": "markdown", + "id": "b135e183-3cf6-4e46-9fbe-7af2e4439982", + "metadata": {}, + "source": [ + "### 5.1 Configure `Detector1Pipeline`\n", + "\n", + "The `Detector1Pipeline` has the following steps available for NIRSpec FS:\n", + "\n", + "> * `group_scale` : Rescales pixel values to correct for improper onboard frame averaging.\n", + "> * `dq_init` : Initializes the data quality (DQ) flags for the input data.\n", + "> * `saturation` : Flags pixels at or below the A/D floor or above the saturation threshold.\n", + "> * `superbias` : Subtracts the superbias reference file from the input data.\n", + "> * `refpix` : Use reference pixels to correct bias drifts.\n", + "> * `linearity` : Applies a correction for non-linear detector response. \n", + "> * `dark_current` : Subtracts the dark current reference file from the input data.\n", + "> * `jump` : Performs CR/jump detection on each ramp integration within an exposure.\n", + "> * `clean_flicker_noise`: Removes flicker (1/f) noise from calibrated ramp images (similar to `nsclean` in spec2).\n", + "> * `ramp_fit` : Determines the mean count rate (counts per second) for each pixel by performing a linear fit to the input data.\n", + "> * `gain_scale` : Corrects pixel values for non-standard gain settings, primarily in NIRSpec subarray data.\n", + "\n", + "For more information about each step and a full list of step arguments, please refer to the official documentation: [JDox](https://jwst-docs.stsci.edu/jwst-science-calibration-pipeline-overview/stages-of-jwst-data-processing/calwebb_detector1) and\n", + "[ReadtheDocs](https://jwst-pipeline.readthedocs.io/en/stable/jwst/pipeline/calwebb_detector1.html)\n", + "\n", + "Below, we set up a dictionary that defines how the `Detector1Pipeline` should be configured for FS data. " + ] + }, + { + "cell_type": "markdown", + "id": "a36fab53-21c8-46e8-85a8-db632520c0fb", + "metadata": {}, + "source": [ + "
\n", + " To override specific steps and reference files, use the examples below.\n", + "
\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "449e61a0-2c0d-4469-a719-078ad33b2938", + "metadata": {}, + "outputs": [], + "source": [ + "# Set up a dictionary to define how the Detector1 pipeline should be configured.\n", + "\n", + "# -------------------------Boilerplate dictionary setup-------------------------\n", + "det1dict = {}\n", + "det1dict['group_scale'], det1dict['dq_init'], det1dict['saturation'] = {}, {}, {}\n", + "det1dict['superbias'], det1dict['refpix'] = {}, {}\n", + "det1dict['linearity'], det1dict['dark_current'] = {}, {}\n", + "det1dict['jump'], det1dict['ramp_fit'], det1dict['gain_scale'] = {}, {}, {}\n", + "\n", + "# ---------------------------Override reference files---------------------------\n", + "\n", + "# Overrides for various reference files (example).\n", + "# Files should be in the base local directory or provide full path.\n", + "#det1dict['dq_init']['override_mask'] = 'myfile.fits' # Bad pixel mask\n", + "#det1dict['superbias']['override_superbias'] = 'myfile.fits' # Bias subtraction\n", + "#det1dict['dark_current']['override_dark'] = 'myfile.fits' # Dark current subtraction\n", + "\n", + "# -----------------------------Set step parameters------------------------------\n", + "\n", + "# Overrides for whether or not certain steps should be skipped (example).\n", + "det1dict['linearity']['skip'] = False # This is the default.\n", + "\n", + "# Turn on multi-core processing (off by default).\n", + "# Choose what fraction of cores to use (quarter, half, or all).\n", + "det1dict['jump']['maximum_cores'] = 'half'\n", + "#det1dict['ramp_fit']['maximum_cores'] = 'half'" + ] + }, + { + "cell_type": "markdown", + "id": "c9da565f-1a16-4853-ba6c-1552d361b907", + "metadata": {}, + "source": [ + "
\n", + " \n", + "Many exposures are affected by artifacts known as [snowballs](https://jwst-docs.stsci.edu/known-issues-with-jwst-data/shower-and-snowball-artifacts#gsc.tab=0) caused by large cosmic ray events. These artifacts are particularly significant in deep exposures with long integration times, with an estimated rate of one snowball per detector (FULL FRAME) per 20 seconds. To expand the number of pixels flagged as jumps around large cosmic ray events, set `expand_large_events` to True. An `expand_factor` of 3 works well for NIRSpec observations to cover most snowballs.\n", + "
" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "96bd19ad-c315-404c-93ff-5d7d24ffecc0", + "metadata": {}, + "outputs": [], + "source": [ + "# Turn on detection of cosmic ray snowballs (on by default).\n", + "det1dict['jump']['expand_large_events'] = True\n", + "det1dict['jump']['expand_factor'] = 3 # (default 2)\n", + "\n", + "# Suppress computations for saturated ramps with\n", + "# only one good (unsaturated) sample (default True).\n", + "# The demo data has some saturation.\n", + "det1dict['ramp_fit']['suppress_one_group'] = False" + ] + }, + { + "cell_type": "markdown", + "id": "6f1d283b-d710-45e1-89ca-c4a84b3e2f65", + "metadata": {}, + "source": [ + "
\n", + " \n", + "JWST detector readout electronics (a.k.a. SIDECAR ASICs) generate significant 1/f noise during detector operations and signal digitization. This noise manifests as faint banding along the detector's slow axis and varies from column to column. For NIRSpec data, the primary pipeline algorithm to address 1/f noise is `nsclean` in the `Spec2Pipeline`. (Rauscher 2023) but is off by default.\n", + "\n", + "An additional 1/f noise-cleaning algorithm, [`clean_flicker_noise`](https://jwst-pipeline.readthedocs.io/en/stable/jwst/clean_flicker_noise/main.html#overview), has been implemented at the group stage in the `Detector1Pipeline`. It removes flicker noise from calibrated ramp images after the jump step and prior to performing the ramp_fitting step. This step is also off by default.\n", + "
" + ] + }, + { + "cell_type": "markdown", + "id": "006d3560-e4cb-4c30-bb60-e551d3f9ef26", + "metadata": {}, + "source": [ + "---\n", + "\n", + "### 5.2 Run `Detector1Pipeline`\n", + "\n", + "Run the science files, nods, and, if available, any dedicated background files through the `calwebb_detector1` pipeline using the `.call()` method. \n", + "\n", + "We use `.call()` instead of `.run()` to ensure that the latest default parameters defined via reference files in CRDS, are applied ([ReadtheDocs](https://jwst-pipeline.readthedocs.io/en/latest/jwst/stpipe/call_via_run.html)).\n", + "\n", + "This stage takes approximately 2 minutes to process six `_uncal.fits` files from the demo data." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "87fc7152-75b0-4374-81ce-85b4f3c7e58c", + "metadata": {}, + "outputs": [], + "source": [ + "# Final list of UNCAL files ready for Stage 1 processing.\n", + "uncal_sci = sorted(glob.glob(uncal_dir + '*uncal.fits'))\n", + "uncal_bg = sorted(glob.glob(uncal_bgdir + '*uncal.fits'))\n", + "\n", + "print(f\"Science UNCAL Files:\\n{'-'*20}\\n\" + \"\\n\".join(uncal_sci))\n", + "print(f\"Background UNCAL Files:\\n{'-'*20}\\n\" + \"\\n\".join(uncal_bg))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "a4ea658b-a4ad-4e78-9eea-d8f57596f741", + "metadata": {}, + "outputs": [], + "source": [ + "time_det1 = time.perf_counter() # Tracks runtime for Stage 1." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "0081dc17-5285-4512-8d0c-d6fe6e817048", + "metadata": { + "scrolled": true + }, + "outputs": [], + "source": [ + "# Run the Stage 1 pipeline with the custom det1dict dictionary.\n", + "# Specify the following:\n", + "# - output_dir: Directory to save *_rate.fits files.\n", + "# - save_results: Set to True to ensure rate files are saved.\n", + "\n", + "if dodet1:\n", + " # --------------------------Science UNCAL files--------------------------\n", + " for uncal_file in uncal_sci:\n", + " print(f\"Applying Stage 1 Corrections & Calibrations to: \"\n", + " f\"{os.path.basename(uncal_file)}\")\n", + "\n", + " det1_result = Detector1Pipeline.call(uncal_file,\n", + " save_results=True,\n", + " steps=det1dict,\n", + " output_dir=det1_dir)\n", + " print(\"Stage 1 has been completed!\\n\")\n", + "else:\n", + " print(\"Skipping Stage 1 processing for Science data.\\n\")\n", + "\n", + "\n", + "if (pixel_bg or master_bg) and dodet1:\n", + " # ------------------------Background UNCAL files-------------------------\n", + " for uncal_file in uncal_bg:\n", + " print(f\"Applying Stage 1 Corrections & Calibrations to: \"\n", + " f\"{os.path.basename(uncal_file)}\")\n", + "\n", + " det1bg_result = Detector1Pipeline.call(uncal_file,\n", + " save_results=True,\n", + " steps=det1dict,\n", + " output_dir=det1_bgdir)\n", + " print(\"Stage 1 for background has been completed!\\n\")\n", + "else:\n", + " print(\"Skipping Stage 1 processing for background data.\\n\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "48e779c4-c972-4f30-8930-ef9a299e104f", + "metadata": {}, + "outputs": [], + "source": [ + "# Print out the time benchmark.\n", + "time1 = time.perf_counter()\n", + "print(f\"Runtime so far: {round((time1-time0)/60.0, 1):0.4f} min\")\n", + "print(f\"Runtime for Det1: {round((time1-time_det1)/60.0, 1):0.4f} min\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "abdf1b4f-6081-4cb9-89d8-b0d55ba1831b", + "metadata": {}, + "outputs": [], + "source": [ + "# Print output result details:\n", + "#det1_result.__dict__ # View entire contents.\n", + "#det1_result.meta.filename\n", + "#det1_result.data.shape" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "cff537e6-0519-45ea-b050-f9c2e0436c1b", + "metadata": {}, + "outputs": [], + "source": [ + "# Final list of RATE[INTS] files ready for Stage 2 processing.\n", + "rate_sci = sorted(glob.glob(det1_dir + '*_rate.fits'))\n", + "rate_bg = sorted(glob.glob(det1_bgdir + '*_rate.fits'))\n", + "rateints_sci = sorted(glob.glob(det1_dir + '*_rateints.fits'))\n", + "rateints_bg = sorted(glob.glob(det1_bgdir + '*_rateints.fits'))\n", + "\n", + "print(f\"SCIENCE | RATE[INTS] Files:\\n{'-'*30}\\n\" + \"\\n\".join(rate_sci + rateints_sci))\n", + "print(f\"BACKGROUND | RATE[INTS] Files:\\n{'-'*30}\\n\" + \"\\n\".join(rate_bg + rateints_bg))" + ] + }, + { + "cell_type": "markdown", + "id": "56207fe4-9b04-489e-8e72-8c208d4e6ba9", + "metadata": {}, + "source": [ + "---\n", + "\n", + "## 6. Stage 2: `Spec2Pipeline` (`calwebb_spec2`)\n", + "\n", + "In this section, we process our countrate (slope) image products from Stage 1 (`calwebb_detector1`) through the Spec2 (`calwebb_spec2`) pipeline to create Stage 2 data products.\n", + "\n", + "* **Input**: A single countrate (slope) image (`_rate[ints].fits`) or an association file listing multiple inputs.\n", + "* **Output**: Calibrated products (rectified and unrectified) and 1D spectra.\n", + "\t* `_cal[ints].fits`: Calibrated 2D (unrectified) spectra (ncols x nrows).\n", + "\t* `_s2d.fits`: Resampled (rectified) 2D spectra (ncols x nrows). \n", + "\t* `_x1d[ints].fits`: Extracted 1D spectroscopic data (wavelength vs. flux).\n", + " \n", + "The `Spec2Pipeline` applies additional instrumental corrections and calibrations (e.g., slit loss, path loss, etc.,) to countrate products that result in a fully calibrated individual exposure (per nod/dither position). The `Spec2Pipeline` also converts countrate products from units of DN/s to flux (Jy) for point sources and surface brightness (MJy/sr) for extended sources.\n", + "\n", + "---" + ] + }, + { + "cell_type": "markdown", + "id": "68f1e8f3-874f-48f3-9f20-9125ebda6f5e", + "metadata": {}, + "source": [ + "### 6.1 Configure `Spec2Pipeline`\n", + "\n", + "The `Spec2Pipeline` has the following steps available for NIRSpec FS:\n", + "\n", + "> * `assign_wcs`: Assigns wavelength solution for spectra.\n", + "> * `nsclean`: Cleans 1/f noise.\n", + "> * `bkg_subtract`: Performs image subtraction for background removal.\n", + "> * `extract_2d` : Extracts 2D arrays from spectral images.\n", + "> * `srctype`: Determines whether a spectroscopic source should be classified as a point or extended object.\n", + "> * `wavecorr` : Updates wavelengths for FS and MOS point sources that are offset in the dispersion direction within their slit.\n", + "> * `flat_field`: Applies flat-field corrections to the input science dataset.\n", + "> * `pathloss`: Calculates and applies corrections for signal loss in spectroscopic data.\n", + "> * `photom`: Applies photometric calibrations to convert data from countrate to surface brightness or flux density.\n", + "> * `pixel_replace`: Interpolates and estimates flux values for pixels flagged as DO_NOT_USE in 2D extracted spectra.\n", + "> * `resample_spec`: Resamples each input 2D spectral image using WCS and distortion information.\n", + "> * `extract_1d`: Extracts a 1D signal from 2D or 3D datasets.\n", + "\n", + "For more information about each step and a full list of step arguments, please refer to the official documentation: [JDox](https://jwst-docs.stsci.edu/jwst-science-calibration-pipeline-overview/stages-of-jwst-data-processing/calwebb_spec2) and\n", + "[ReadtheDocs](https://jwst-pipeline.readthedocs.io/en/latest/jwst/pipeline/calwebb_spec2.html)\n", + "\n", + "Below, we set up a dictionary that defines how the `Spec2Pipeline` should be configured for FS data. \n" + ] + }, + { + "cell_type": "markdown", + "id": "14bd43d0-91d2-4fe5-9876-6317ad81b913", + "metadata": {}, + "source": [ + "
\n", + "\n", + "If pixel-to-pixel background subtraction was chosen above, it will be applied during this stage.
\n", + "To override specific steps and reference files, use the examples below. \n", + "\n", + "
" + ] + }, + { + "cell_type": "markdown", + "id": "6fc9dc54-0f58-45e4-96c3-1a86764d2f74", + "metadata": {}, + "source": [ + "
\n", + " \n", + "Resampling 2D spectra can sometimes introduce artificial noise and reduce the signal-to-noise ratio (SNR) in the resulting 1D spectra when using `weight_type='ivm'` ([see known issues](https://jwst-docs.stsci.edu/known-issues-with-jwst-data/nirspec-known-issues/nirspec-fs-known-issues#NIRSpecFSKnownIssues-Resamplingof2-Dspectra:~:text=noise%20workaround%20notebook.-,Resampling%20of%202%2DD%20spectra,-Stages%202%20and)). The default is now set to set weight type to 'exptime'. Consider the following when selecting a `weight_type` parameter:\n", + "* **'ivm'**: Inverse variant scaling based on read noise (VAR_RNOISE), ideal for rejecting outliers and better suited for faint sources.\n", + "* **'exptime'**: Uses exposure time for scaling, improving SNR for bright sources.\n", + "
\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "c44854a3-1cff-43a8-a115-9a1a6414f014", + "metadata": {}, + "outputs": [], + "source": [ + "# Set up a dictionary to define how the Spec2 pipeline should be configured.\n", + "\n", + "# -------------------------Boilerplate dictionary setup-------------------------\n", + "spec2dict = {}\n", + "spec2dict['assign_wcs'], spec2dict['nsclean'] = {}, {}\n", + "spec2dict['extract_2d'], spec2dict['bkg_subtract'] = {}, {}\n", + "spec2dict['srctype'], spec2dict['wavecorr'] = {}, {}\n", + "spec2dict['flat_field'], spec2dict['pathloss'] = {}, {}\n", + "spec2dict['photom'], spec2dict['pixel_replace'] = {}, {}\n", + "spec2dict['resample_spec'], spec2dict['extract_1d'] = {}, {}\n", + "\n", + "# ---------------------------Override reference files---------------------------\n", + "\n", + "# Overrides for various reference files (example).\n", + "# Files should be in the base local directory or provide full path.\n", + "#spec2dict['extract_1d']['override_extract1d'] = 'myfile.json'\n", + "\n", + "# -----------------------------Set step parameters------------------------------\n", + "\n", + "# Overrides for whether or not certain steps should be skipped (example).\n", + "spec2dict['bkg_subtract']['skip'] = not pixel_bg # Runs if pixel-to-pixel bkg selected.\n", + "\n", + "# Run pixel replacement code to extrapolate values for otherwise bad pixels.\n", + "# This can help mitigate 5-10% negative dips in spectra of bright sources.\n", + "# Use the 'fit_profile' algorithm.\n", + "#spec2dict['pixel_replace']['skip'] = False\n", + "#spec2dict['pixel_replace']['n_adjacent_cols'] = 5\n", + "#spec2dict['pixel_replace']['algorithm'] = 'fit_profile'\n", + "\n", + "# Run nsclean for 1/f noise.\n", + "#spec2dict['nsclean']['skip'] = False\n", + "#spec2dict['nsclean']['n_sigma'] = 2\n", + "\n", + "# Resample weight_type.\n", + "#spec2dict['resample_spec']['weight_type'] = 'exptime'" + ] + }, + { + "cell_type": "markdown", + "id": "ce8ff935-bc98-4664-ae55-427954a419cb", + "metadata": {}, + "source": [ + "
\n", + "\n", + "To correct for 1/f noise with `nsclean` in Stage 2, see the **FS_NSClean_example** demo notebook for FS data [here](https://github.com/spacetelescope/jdat_notebooks/tree/main/notebooks/NIRSpec/NIRSpec_NSClean).\n", + "\n", + "
" + ] + }, + { + "cell_type": "markdown", + "id": "1b5cda90-50c3-42f4-8b40-1fe9404da201", + "metadata": {}, + "source": [ + "---\n", + "\n", + "### 6.2 Create `Spec2Pipeline` Association Files\n", + "\n", + "[Association (ASN) files](https://jwst-pipeline.readthedocs.io/en/stable/jwst/associations/overview.html) define the relationships between multiple exposures, allowing them to get processed as a set rather than individually. Processing an ASN file enables the exposures to be calibrated, archived, retrieved, and reprocessed as a set rather than as individual objects.\n", + "\n", + "[Stage 2 ASN files](https://jwst-pipeline.readthedocs.io/en/latest/jwst/associations/level2_asn_technical.html) for FS data can include `science` and `background` exposure types. A Stage 2 ASN file requires at least one `science` file but can contain multiple `background` files that enable pixel-to-pixel background subtraction in `calwebb_spec2`.\n", + "\n", + "Types of Associations:\n", + "* Observation: Includes data from a single observation.\n", + "* Candidate: Includes data from multiple observations (e.g., dedicated background exposures).\n", + "\n", + "In practice, Stage 2 ASN files can be downloaded directly from MAST, as we did with the demo data for convenience. Below, we provide an example of manually creating Stage 2 ASN files using the `write_asn` function defined above. To run the function, set the `spec2_asn_manual` flag below.\n", + "\n", + "
\n", + "\n", + "Background subtraction may not be correctly applied if more than *one* `science` file is included in the association. Additionally, pixel-to-pixel background subtraction will only be performed if the grating wheel has not moved between the target and off-scene associated background exposures. If the grating wheel moved between the target and background exposures (as would be the case if they were in different visits), pipeline processing will follow a more involved \"master background\" subtraction done in Stage 3.\n", + "\n", + "
" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "717ae8ae-a506-4815-8759-00ee6116e8e7", + "metadata": {}, + "outputs": [], + "source": [ + "# Set `spec2_manual` to True to create Stage 2 ASN files.\n", + "spec2_asn_manual = False\n", + "if dospec2 and spec2_asn_manual:\n", + " write_asn(rate_sci, background_files=rate_bg, asn_output_dir=asn_dir, level=2)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "36d229cf-23f2-4944-9e5e-9b2116f2f2d6", + "metadata": {}, + "outputs": [], + "source": [ + "# Get list of all spec2 association (ASN) files and categorize them.\n", + "spec2_asn_all = glob.glob(f\"{asn_dir}*spec2*asn.json\")\n", + "candidate_asn = [asn for asn in spec2_asn_all if \"-c\" in asn]\n", + "obs_asn = [asn for asn in spec2_asn_all if \"-o\" in asn and \"man\" not in asn]\n", + "\n", + "# If ASN files were created manually using the function above:\n", + "man_asn = [asn for asn in spec2_asn_all if \"manual\" in asn]\n", + "\n", + "# Choose ASN files: candidate if available, otherwise observation ASN.\n", + "spec2_asn = man_asn if spec2_asn_manual else candidate_asn or obs_asn\n", + "print(f\"Stage 2 ASN Files:\\n{'-'*20}\\n\" + \"\\n\".join(spec2_asn))" + ] + }, + { + "cell_type": "markdown", + "id": "98ad3e07-8ac6-4df8-8530-a4a1965e6474", + "metadata": {}, + "source": [ + "ASN files downloaded from MAST expect the input files to be in the same directory, which is incompatible with our directory structure. In the cell below, we update the `expname` fields in the ASN files to use absolute paths to ensure the pipeline looks in the correct locations.\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "3fc12ab2-2179-4ec1-bf82-7c3162cf749e", + "metadata": {}, + "outputs": [], + "source": [ + "# Convert 'expname' paths in the ASN file to absolute paths.\n", + "# Ensures the pipeline can locate the files,\n", + "# regardless of the ASN file's location.\n", + "update_asn_paths(spec2_asn, exclude_dirs=[])" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "2900c561-9dee-4d33-8e77-c808ba1e086c", + "metadata": { + "scrolled": true + }, + "outputs": [], + "source": [ + "# Open an ASN file as an example.\n", + "# Check that file paths have been correctly updated.\n", + "with open(spec2_asn[0], 'r') as f_obj:\n", + " asnfile_data = json.load(f_obj)\n", + "\n", + "JSON(asnfile_data, expanded=True)" + ] + }, + { + "cell_type": "markdown", + "id": "25fcc478-5310-4d82-a16a-68abe63c61cc", + "metadata": {}, + "source": [ + "---\n", + "\n", + "### 6.3 Run `Spec2Pipeline`\n", + "\n", + "Run the science files, associated nods and, if available, any background files through the `calwebb_spec2` pipeline using the `.call()` method." + ] + }, + { + "cell_type": "markdown", + "id": "cca1b66b-a7ce-4ebd-b46c-84c0b76e5e0d", + "metadata": {}, + "source": [ + "
\n", + "Perform pixel-to-pixel background subtraction (if desired) here in Stage 2. Otherwise, reduce the backgrounds individually for master background subtraction in Stage 3 (if desired).\n", + "
" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "746d7456-efa2-43fc-a2d3-1c9486efc116", + "metadata": {}, + "outputs": [], + "source": [ + "time_spec2 = time.perf_counter()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "bf1c24c9-1de0-4c0e-8882-ecb712d25067", + "metadata": {}, + "outputs": [], + "source": [ + "# Process dedicated background data to use in the master background step.\n", + "\n", + "if dospec2 and master_bg:\n", + " # ------------------------Background RATE files------------------------\n", + " for rate in rate_bg:\n", + " print(f\"Applying Stage 2 Corrections & Calibrations to: \"\n", + " f\"{os.path.basename(rate)}\")\n", + "\n", + " spec2bg_result = Spec2Pipeline.call(rate,\n", + " save_results=True,\n", + " steps=spec2dict,\n", + " output_dir=spec2_bgdir)\n", + "else:\n", + " print(\"Skipping Spec2 processing for background data.\")" + ] + }, + { + "cell_type": "markdown", + "id": "5cfbb257-dc2a-47e1-bff6-1f2d7690de9a", + "metadata": {}, + "source": [ + "Process the science data." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "147be713-5c4c-46f7-ab84-ca77b877cfa6", + "metadata": {}, + "outputs": [], + "source": [ + "# To save on runtime turns off creation of quicklook 2d/1d spectra for science data.\n", + "spec2dict['resample_spec']['skip'] = False # S2D products.\n", + "spec2dict['extract_1d']['skip'] = False # X1D products." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "9dbaa79f-4a3c-4f11-80ce-d060af0b5f60", + "metadata": { + "scrolled": true + }, + "outputs": [], + "source": [ + "# Run Stage 2 pipeline.\n", + "\n", + "if dospec2:\n", + " # --------------------------Science ASN files--------------------------\n", + " for asn in spec2_asn:\n", + "\n", + " asn_data = json.load(open(asn))\n", + " sci_file = os.path.basename(asn_data['products'][0]['members'][0]['expname'])\n", + " print(f\"Applying Stage 2 Corrections & Calibrations to: {sci_file}\")\n", + "\n", + " spec2sci_result = Spec2Pipeline.call(asn,\n", + " save_results=True,\n", + " steps=spec2dict,\n", + " output_dir=spec2_dir)\n", + "else:\n", + " print(\"Skipping Spec2 processing.\")\n", + "\n", + "print(\" ... Stage 2 has been completed!\\n\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "7ca91e64-2fb9-46d8-8b95-b8e41b3c0016", + "metadata": {}, + "outputs": [], + "source": [ + "# Print output result details:\n", + "#spec2sci_result.__dict__ # View entire contents.\n", + "#spec2sci_result.meta.filename\n", + "#spec2sci_result.data.shape" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "a1b42e27-aa6b-45c4-baec-9c1157329252", + "metadata": {}, + "outputs": [], + "source": [ + "# Print out the time benchmarks.\n", + "time2 = time.perf_counter()\n", + "print(f\"Runtime so far: {round((time2-time0)/60.0, 1):0.4f} min\")\n", + "print(f\"Runtime for Spec2: {round((time2-time_spec2)/60.0, 1):0.4f} min\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "c553a0f5-5ca8-4a9c-905c-c7897e7dde8c", + "metadata": {}, + "outputs": [], + "source": [ + "# List the Stage 2 products.\n", + "\n", + "# -----------------------------Science files-----------------------------\n", + "sci_cal = sorted(glob.glob(spec2_dir + '*_cal.fits'))\n", + "sci_s2d = sorted(glob.glob(spec2_dir + '*_s2d.fits'))\n", + "sci_x1d = sorted(glob.glob(spec2_dir + '*_x1d.fits'))\n", + "\n", + "print(f\"SCIENCE | Stage 2 CAL Products:\\n{'-'*20}\\n\" + \"\\n\".join(sci_cal))\n", + "print(f\"SCIENCE | Stage 2 S2D Products:\\n{'-'*20}\\n\" + \"\\n\".join(sci_s2d))\n", + "print(f\"SCIENCE | Stage 2 X1D Products:\\n{'-'*20}\\n\" + \"\\n\".join(sci_x1d))\n", + "\n", + "# ----------------------------Background files---------------------------\n", + "bg_cal = sorted(glob.glob(spec2_bgdir + '*_cal.fits'))\n", + "bg_s2d = sorted(glob.glob(spec2_bgdir + '*_s2d.fits'))\n", + "bg_x1d = sorted(glob.glob(spec2_bgdir + '*_x1d.fits'))\n", + "\n", + "print(f\"BACKGROUND | Stage 2 CAL Products:\\n{'-'*20}\\n\" + \"\\n\".join(bg_cal))\n", + "print(f\"BACKGROUND | Stage 2 S2D Products:\\n{'-'*20}\\n\" + \"\\n\".join(bg_s2d))\n", + "print(f\"BACKGROUND | Stage 2 X1D Products:\\n{'-'*20}\\n\" + \"\\n\".join(bg_x1d))" + ] + }, + { + "cell_type": "markdown", + "id": "73ea2288-f1fb-412a-a0c4-20c56a66cc66", + "metadata": {}, + "source": [ + "
\n", + "\n", + "**Flux Calibration Updates (Build 11.1)**: Differences in flux calibration are expected due to recent improvements, including updates to spectral resampling, new F-flat reference files, and improved pathloss corrections for FS data. These changes will be seen in both Stage 2 and Stage 3 products. \n", + "\n", + "> * **Improved Spectral Resampling (`resample_spec`)**: The spectral resampling for NIRSpec was updated to use a grid uniformly sampled in angular units rather than normalized slit units. This update ensures flux is conserved across different slitlet sizes and 2D extraction box dimensions, resulting in improved consistency for MOS slitlets with varying numbers of shutters.\n", + "> * **New MOS/FS F-flat Reference Files**: New F-flat reference files for MOS and FS were delivered to CRDS to work with the updated spectral resampling. These files correct an issue that would have caused an erroneous ~5% change in flux calibration with the new resampling method.\n", + "> * **S1600A1 Calibration Update**: The F-flat for S1600A1 was additionally adjusted to accurately calibrate un-dithered observations (e.g., BOTS) by accounting for over-subtraction caused by neighboring nods. \n", + "> * **Improved Hybrid (Data+Model) FS Pathlosses for Cycle 1/2+ Data**: New pathloss reference files for FS data were delivered to CRDS. These files use data-derived pathlosses at the fixed dither positions and model-based pathlosses elsewhere for slits S200A1/A2. Pathloss corrections are turned off for slits S400A1 and S1600A1, as they are consistent within 1 sigma of observed data at fixed dither positions (minimal pathloss).\n", + "\n", + "\n", + "
" + ] + }, + { + "cell_type": "markdown", + "id": "46615dc5-6efd-4b1f-8d15-bb6e3e2d9eeb", + "metadata": {}, + "source": [ + "---\n", + "\n", + "## 7. Stage 3: `Spec3Pipeline` (`calwebb_spec3`)\n", + "\n", + "In this section, we process our calibrated spectra from Stage 2 (`calwebb_spec2`) through the Spec3 (`calwebb_spec3`) pipeline to create Stage 3 data products.\n", + "\n", + "* **Input**: An ASN file that lists multiple calibrated exposures (`_cal.fits`) in addition to any background exposures (`_x1d.fits`).\n", + "* **Output**: A single calibrated product (rectified and unrectified) and 1D spectrum. These data products have units of MJy/sr (or Jy for extracted point-source spectra).\n", + "\t* `_cal.fits`: Calibrated 2D (unrectified) spectra (ncols x nrows).\n", + " * `_crf.fits`: Calibrated 2D (unrectified) spectra whose DQ array has been updated to flag pixels detected as outliers (ncols x nrows).\n", + " * `_s2d.fits`: Resampled (rectified) 2D spectra (ncols x nrows). \n", + "\t* `_x1d.fits`: Extracted 1D spectroscopic data. \n", + "\n", + "The `Spec3Pipeline` performs additional corrections (e.g., outlier detection, background subtraction) and combines calibrated data from multiple exposures (e.g. a dither/nod pattern) into a single 2D spectral product, as well as a combined 1D spectrum." + ] + }, + { + "cell_type": "markdown", + "id": "9262258e-2fc3-4c83-92e9-0f3a74cd4f86", + "metadata": {}, + "source": [ + "---\n", + "\n", + "### 7.1 Configure `Spec3Pipeline`\n", + "\n", + "The `Spec3Pipeline` has the following steps available for NIRSpec FS:\n", + "\n", + "> * `assign_mtwcs`: Modifies the WCS output frame in each exposure of a Moving Target (MT) observation association.\n", + "> * `master_background`: Master background subtraction.\n", + "> * `outlier_detection` : Identification of bad pixels or cosmic-rays that remain in each of the input images.\n", + "> * `pixel_replace`: Interpolates and estimates flux values for pixels flagged as DO_NOT_USE in 2D extracted spectra.\n", + "> * `resample_spec`: Resamples each input 2D spectral image using WCS and distortion information.\n", + "> * `extract_1d`: Extracts a 1D signal from 2D or 3D datasets\n", + "\n", + "For more information about each step and a full list of step arguments, please refer to the official documentation: [JDox](https://jwst-docs.stsci.edu/jwst-science-calibration-pipeline-overview/stages-of-jwst-data-processing/calwebb_spec3) and\n", + "[ReadtheDocs](https://jwst-pipeline.readthedocs.io/en/latest/jwst/pipeline/calwebb_spec3.html)\n", + "\n", + "Below, we set up a dictionary that defines how the `Spec3Pipeline` should be configured for FS data. \n", + "\n", + "\n", + "
\n", + "\n", + "If master background subtraction was chosen above, it will be applied during this stage.
\n", + "To override specific steps and reference files, use the examples below. \n", + "\n", + "
" + ] + }, + { + "cell_type": "markdown", + "id": "5102e7a7-876f-4b60-84c9-0e603489ae47", + "metadata": {}, + "source": [ + "
\n", + " \n", + "The outlier detection step can be too aggressive for FS data, and better results may be obtained by turning it off.\n", + "
" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "a734be27-787f-4273-abf7-827c7b552e77", + "metadata": {}, + "outputs": [], + "source": [ + "# Set up a dictionary to define how the Spec3 pipeline should be configured.\n", + "\n", + "# -------------------------Boilerplate dictionary setup-------------------------\n", + "spec3dict = {}\n", + "spec3dict['assign_mtwcs'], spec3dict['master_background'] = {}, {}\n", + "spec3dict['outlier_detection'], spec3dict['pixel_replace'] = {}, {}\n", + "spec3dict['resample_spec'], spec3dict['extract_1d'] = {}, {}\n", + "\n", + "# ---------------------------Override reference files---------------------------\n", + "\n", + "# Overrides for various reference files.\n", + "# Files should be in the base local directory or provide full path.\n", + "# spec3dict['extract_1d']['override_extract1d'] = 'myfile.json'\n", + "\n", + "# -----------------------------Set step parameters------------------------------\n", + "\n", + "# Overrides for whether or not certain steps should be skipped (example).\n", + "spec3dict['outlier_detection']['skip'] = True\n", + "\n", + "# Master background usage was set up above, propagate that here.\n", + "spec3dict['master_background']['skip'] = not master_bg\n", + "\n", + "# Run pixel replacement code to extrapolate values for otherwise bad pixels.\n", + "# This can help mitigate 5-10% negative dips in spectra of bright sources.\n", + "# Use the 'fit_profile' algorithm.\n", + "# spec3dict['pixel_replace']['skip'] = False\n", + "# spec3dict['pixel_replace']['n_adjacent_cols'] = 5\n", + "# spec3dict['pixel_replace']['algorithm'] = 'fit_profile'\n", + "\n", + "# Resample weight_type.\n", + "# spec3dict['resample_spec']['weight_type'] = 'exptime'" + ] + }, + { + "cell_type": "markdown", + "id": "81f331fc-f869-4fcb-a796-548545253de1", + "metadata": {}, + "source": [ + "---\n", + "### 7.2 Create `Spec3Pipeline` Association Files\n", + "\n", + "[Stage 3 ASN files](https://jwst-pipeline.readthedocs.io/en/latest/jwst/associations/level3_asn_technical.html) for FS data can include `science` and `background` exposure types. A Stage 3 ASN file requires at least one `science` file (there is usually more than one) but can contain multiple `background` files that enable master background subtraction in `calwebb_spec3`. **Note that the science exposures should be in the `_cal.fits` format, while the background exposures must be in the `_x1d.fits` format.**\n", + "\n", + "In practice, Stage 3 ASN files can be downloaded directly from MAST, as we did with the demo data for convenience. Below, we provide an example of manually creating Stage 3 ASN files using the `write_asn` function defined above. To run the function, set the `spec3_manual` flag below.\n", + "\n", + "\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "c8fc6718-8e06-463a-ac96-4a76f0354da8", + "metadata": {}, + "outputs": [], + "source": [ + "# Set `spec3_manual` to True to create Stage 3 ASN files.\n", + "spec3_manual = False\n", + "if dospec3 and spec3_manual:\n", + " write_asn(sci_cal, background_files=bg_x1d, asn_output_dir=asn_dir, level=3)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "b58bb6b9-56f0-40f2-bb6d-d9325ced4102", + "metadata": {}, + "outputs": [], + "source": [ + "# Get list of all spec3 association (ASN) files and categorize them.\n", + "spec3_asn_all = glob.glob(f\"{asn_dir}*spec3*asn.json\")\n", + "spec3_obs_asn = [asn for asn in spec3_asn_all if \"-o\" in asn and \"manual\" not in asn]\n", + "\n", + "# If ASN files were created manually using the function above:\n", + "spec3_asn_man = [asn for asn in spec3_asn_all if \"manual\" in asn]\n", + "\n", + "# Choose ASN files: manually created or downloaded from MAST.\n", + "spec3_asn = spec3_asn_man if spec3_manual else spec3_obs_asn\n", + "print(f\"Stage 3 ASN Files:\\n{'-'*20}\\n\" + \"\\n\".join(spec3_asn))" + ] + }, + { + "cell_type": "markdown", + "id": "9cbbb928-1f82-4817-967b-62b2420b46a4", + "metadata": {}, + "source": [ + "ASN files downloaded from MAST expect the input files to be in the same directory, which is incompatible with our directory structure. In the cell below, we update the `expname` fields in the ASN files to use absolute paths to ensure the pipeline looks in the correct locations." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "b19744a0-50a7-41c7-a2ad-b7988696daaa", + "metadata": {}, + "outputs": [], + "source": [ + "# Convert 'expname' paths in the ASN file to absolute paths.\n", + "# Ensures the pipeline can locate the files,# regardless of the ASN file's location.\n", + "update_asn_paths(spec3_asn, exclude_dirs=[])" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "b697a75b-3f6c-4509-b5a0-6c5cb0dacda1", + "metadata": { + "scrolled": true + }, + "outputs": [], + "source": [ + "# Open an ASN file as an example.\n", + "# Check that file paths have been correctly updated.\n", + "with open(spec3_asn[0], 'r') as f_obj:\n", + " asnfile_data = json.load(f_obj)\n", + "\n", + "JSON(asnfile_data, expanded=True)" + ] + }, + { + "cell_type": "markdown", + "id": "b6d3a15f-8ade-44e5-9256-9bc0103fcef4", + "metadata": {}, + "source": [ + "---\n", + "\n", + "### 7.3 Run `Spec3Pipeline`\n", + "\n", + "Run the science files and, if available, any background files through the `calwebb_spec3` pipeline using the `.call()` method." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "232579bd-5948-42c1-a00c-03b924c30d03", + "metadata": {}, + "outputs": [], + "source": [ + "time_spec3 = time.perf_counter()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "e08f8275-3c16-4a40-bf08-b99607ebce7d", + "metadata": { + "scrolled": true + }, + "outputs": [], + "source": [ + "# Run Stage 3 pipeline using the custom spec3dict dictionary.\n", + "\n", + "start = time.time()\n", + "\n", + "if dospec3:\n", + " # --------------------------Spec3 ASN files--------------------------\n", + " for s3_asn in spec3_asn:\n", + " print(f\"Applying Stage 3 Corrections & Calibrations to: \"\n", + " f\"{os.path.basename(s3_asn)}\")\n", + " spec3_result = Spec3Pipeline.call(s3_asn,\n", + " save_results=True,\n", + " steps=spec3dict,\n", + " output_dir=spec3_dir)\n", + "else:\n", + " print(\"Skipping Spec3 processing.\")\n", + "\n", + "print(\" ... Stage 3 has been completed!\\n\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "91e6e577-6b6c-4390-b634-4532667b7056", + "metadata": {}, + "outputs": [], + "source": [ + "# Print output result details:\n", + "#spec3_result.__dict__ # View entire contents.\n", + "#spec3_result.meta.filename\n", + "#spec3_result.data.shape" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "4812af22-170f-4c3b-b768-fa62af1b925d", + "metadata": {}, + "outputs": [], + "source": [ + "# Print out the time benchmarks.\n", + "time3 = time.perf_counter()\n", + "print(f\"Runtime so far: {round((time3-time0)/60.0, 1):0.4f} min\")\n", + "print(f\"Runtime for Spec2: {round((time3-time_spec3)/60.0, 1):0.4f} min\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "01f03261-b84e-4e5e-b49a-be80f2eb5272", + "metadata": {}, + "outputs": [], + "source": [ + "# List the Stage 3 products.\n", + "\n", + "stage3_cal = sorted(glob.glob(spec3_dir + '*_cal.fits'))\n", + "stage3_s2d = sorted(glob.glob(spec3_dir + '*_s2d.fits'))\n", + "stage3_x1d = sorted(glob.glob(spec3_dir + '*_x1d.fits'))\n", + "\n", + "print(f\"Stage 3 CAL Products:\\n{'-'*20}\\n\" + \"\\n\".join(stage3_cal))\n", + "print(f\"Stage 3 S3D Products:\\n{'-'*20}\\n\" + \"\\n\".join(stage3_s2d))\n", + "print(f\"Stage 3 X1D Products:\\n{'-'*20}\\n\" + \"\\n\".join(stage3_x1d))" + ] + }, + { + "cell_type": "markdown", + "id": "14a14495-79f7-4f57-a250-cc1e0a51e1f0", + "metadata": {}, + "source": [ + "---\n", + "\n", + "## 8. Visualize the Data\n", + "\n", + "Here, we'll plot the products for different stages.\n", + "\n", + "First, let's define a few convenience functions." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "cdb5e1b7-493c-48b7-a278-136da6a5a6af", + "metadata": {}, + "outputs": [], + "source": [ + "# Function to display the rate files produced by Stage 1 in the\n", + "# JWST Calibration Pipeline.\n", + "def display_rate(rates,\n", + " slits_models=[],\n", + " integration=0,\n", + " extname='data',\n", + " cmap='viridis',\n", + " bad_color=(1, 0.7, 0.7),\n", + " vmin=None,\n", + " vmax=None,\n", + " scale='asinh',\n", + " aspect='auto',\n", + " title_prefix=None,\n", + " title_path=False,\n", + " save_plot=False):\n", + " \"\"\"\n", + " Display countrate images.\n", + "\n", + " Parameters\n", + " ----------\n", + " rates : list of str\n", + " A list of RATE[INTS] files to be displayed.\n", + " slits_models : list of str, optional\n", + " A list of CAL[INTS] or S2D files containing the slit models.\n", + " If provided, slit cutouts will be overlaid on the countrate images.\n", + " integration : {None, 'min', int}, optional\n", + " Specifies the integration to use for multi-integration data.\n", + " If 'min', the minimum value across all integrations is used.\n", + " If an integer, the specific integration index is used (default 0).\n", + " extname : str, optional\n", + " The name of the data extension to extract from ('data', 'dq', etc.).\n", + " cmap : str, optional\n", + " Colormap to use for displaying the image. Default is 'viridis'.\n", + " bad_color : tuple of float, optional\n", + " Color to use for NaN pixels. Default is light red (1, 0.7, 0.7).\n", + " vmin : float, optional\n", + " Minimum value for color scaling. If None, determined from the data.\n", + " vmax : float, optional\n", + " Maximum value for color scaling. If None, determined from the data.\n", + " scale : {'linear', 'log', 'asinh'}, optional\n", + " Scale to use for the image normalization. Default is 'asinh'.\n", + " aspect : str, optional\n", + " Aspect ratio of the plot. Default is 'auto'.\n", + " title_prefix : str, optional\n", + " Optional prefix for the plot title.\n", + " title_path : bool, optional\n", + " If True, uses the full file path for the title;\n", + " otherwise, uses the basename. Default is False.\n", + " save_plot : bool, optional\n", + " If True, saves the plot as a PNG file. Default is False.\n", + "\n", + " Returns\n", + " -------\n", + " None.\n", + " \"\"\"\n", + "\n", + " # -------------------------------Check Inputs-------------------------------\n", + " rates = [rates] if isinstance(rates, str) else rates\n", + " slits_models = [slits_models] if isinstance(slits_models, str) else slits_models\n", + " nrates = len(rates)\n", + "\n", + " # ------------------------------Set up figures------------------------------\n", + " fig, axes = plt.subplots(nrates, 1, figsize=(12, 12 * nrates),\n", + " sharex=True, height_ratios=[1] * nrates)\n", + " fig.subplots_adjust(hspace=0.2, wspace=0.2)\n", + " axes = [axes] if nrates == 1 else axes\n", + "\n", + " cmap = plt.get_cmap(cmap) # Set up colormap and bad pixel color.\n", + " cmap.set_bad(bad_color, 1.0)\n", + "\n", + " # ---------------------------Plot countrate image---------------------------\n", + " for i, (rate, cal) in enumerate(itertools.zip_longest(rates,\n", + " slits_models,\n", + " fillvalue=None)):\n", + "\n", + " # -------------------Open files as JWST datamodels-------------------\n", + " model = datamodels.open(rate)\n", + " slits_model = datamodels.open(cal) if cal else None\n", + "\n", + " # -----------------------Extract the 2D/3D data----------------------\n", + " data_2d = getattr(model, extname)\n", + " if data_2d.ndim == 3: # Handle multi-integration data.\n", + " if integration == 'min':\n", + " data_2d = np.nanmin(data_2d, axis=0)\n", + " elif isinstance(integration, int) and 0 <= integration < data_2d.shape[0]:\n", + " data_2d = data_2d[integration]\n", + " else:\n", + " raise ValueError(f\"Invalid integration '{integration}' for 3D data.\")\n", + "\n", + " # ---------------------------Scale the data-------------------------\n", + " sigma_clipped_data = sigma_clip(data_2d, sigma=5, maxiters=3)\n", + " vmin = np.nanmin(sigma_clipped_data) if vmin is None else vmin\n", + " vmax = np.nanmax(sigma_clipped_data) if vmax is None else vmax\n", + " stretch_map = {'log': LogStretch(), 'linear': LinearStretch(),\n", + " 'asinh': AsinhStretch()}\n", + " if scale in stretch_map:\n", + " norm = ImageNormalize(sigma_clipped_data,\n", + " interval=ManualInterval(vmin=vmin, vmax=vmax),\n", + " stretch=stretch_map[scale])\n", + " else:\n", + " norm = simple_norm(sigma_clipped_data, vmin=vmin, vmax=vmax)\n", + "\n", + " # ----------------Plot the countrate image & colorbar---------------\n", + " plt.subplots_adjust(left=0.05, right=0.85)\n", + " im = axes[i].imshow(data_2d, origin='lower', cmap=cmap,\n", + " norm=norm, aspect=aspect, interpolation='nearest')\n", + " units = model.meta.bunit_data\n", + " cbar_ax = fig.add_axes([axes[i].get_position().x1 + 0.02,\n", + " axes[i].get_position().y0, 0.02,\n", + " axes[i].get_position().height])\n", + " cbar = fig.colorbar(im, cax=cbar_ax)\n", + " cbar.set_label(units, fontsize=12)\n", + "\n", + " # -----------------Draw slits and label source ids------------------\n", + " # slits_model can be s2d/cal from spec2 - contains slit models for all sources.\n", + " if slits_model:\n", + " slit_patches = []\n", + " for slit in slits_model.slits:\n", + " slit_patch = Rectangle((slit.xstart, slit.ystart),\n", + " slit.xsize, slit.ysize)\n", + " slit_patches.append(slit_patch)\n", + " y = slit.ystart + slit.ysize / 2\n", + " x = slit.xstart if 'nrs1' in rate else slit.xstart + slit.xsize\n", + " ha = 'right' if 'nrs1' in rate else 'left'\n", + " plt.text(x, y, slit.source_id, color='w', ha=ha, va='center',\n", + " fontsize=7, path_effects=[], weight='bold')\n", + " axes[i].add_collection(PatchCollection(slit_patches, ec='r', fc='None'))\n", + "\n", + " # -----------------Construct title and axis labels------------------\n", + " filename = model.meta.filename\n", + " title = (f\"{title_prefix + ' ' if title_prefix else ''}\"\n", + " f\"{filename if title_path else os.path.basename(filename)}\")\n", + " if integration is not None:\n", + " title = title.replace('rateints', f'rateints[{integration}]')\n", + " axes[i].set_title(title, fontsize=14)\n", + " axes[i].set_xlabel(\"Pixel Column\", fontsize=12)\n", + " axes[i].set_ylabel(\"Pixel Row\", fontsize=12)\n", + "\n", + " # -------------------------Save the figure?-------------------------\n", + " if save_plot:\n", + " save_plot = rate.replace('fits', 'png')\n", + " if integration:\n", + " save_plot = save_plot.replace('.png', '%s.png' % integration)\n", + " fig.savefig(save_plot, dpi=200)\n", + "\n", + " fig.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "ecb9caf9-9ceb-478c-9854-fa777d1e8e21", + "metadata": {}, + "outputs": [], + "source": [ + "# Function to display the spectra generated in stage 2/3 of the\n", + "# JWST Calibration Pipeline.\n", + "def display_spectra(spectra,\n", + " compare_x1d=None,\n", + " compare_mast=None,\n", + " integration=None,\n", + " extname='data',\n", + " source_id=1,\n", + " source_type=None,\n", + " expand_wavelength_gap=True,\n", + " plot_resample=True,\n", + " plot_errors=False,\n", + " cmap='viridis',\n", + " bad_color=(1, 0.7, 0.7),\n", + " aspect='auto',\n", + " vmin=None,\n", + " vmax=None,\n", + " scale='asinh',\n", + " title_prefix=None,\n", + " title_path=False,\n", + " y_limits=None,\n", + " is_stage3=False):\n", + "\n", + " \"\"\"\n", + " Display 2D and 1D spectra (Stage 2/3).\n", + "\n", + " Parameters\n", + " ----------\n", + " spectra : list of str\n", + " A list of data products (e.g., CAL, S2D, X1D files).\n", + " compare_x1d : list of str, optional\n", + " A list of 1D spectra for comparison (X1D files).\n", + " compare_mast : list of str, optional\n", + " A list of 1D spectra from MAST for comparison (X1D files).\n", + " integration : {None, 'min', int}, optional\n", + " Specifies the integration to use for multi-integration data.\n", + " If 'min', the minimum value across all integrations is used.\n", + " If an integer, the specific integration index is used (default 0).\n", + " extname : str, optional\n", + " The name of the data extension to extract ('data', 'dq', etc.).\n", + " source_id : int or str, optional\n", + " Identifier for the source/slit to be displayed. Default is 1.\n", + " source_type : str, optional\n", + " Override data source type ('POINT' or 'EXTENDED').\n", + " expand_wavelength_gap : bool, optional\n", + " If True, expands gaps in the wavelength data for better visualization.\n", + " plot_resample : bool, optional\n", + " If True, plots resampled (S2D) data products;\n", + " otherwise, plots calibrated (CAL) data. Default is True.\n", + " plot_errors : bool, optional\n", + " If True, plots the error bands for the 1D spectra. Default is False.\n", + " cmap : str, optional\n", + " Colormap to use for displaying the images. Default is 'viridis'.\n", + " bad_color : tuple of float, optional\n", + " Color to use for bad pixels. Default is light red (1, 0.7, 0.7).\n", + " aspect : str, optional\n", + " Aspect ratio of the plot. Default is 'auto'.\n", + " vmin : float, optional\n", + " Minimum value for color scaling. If None, determined from the data.\n", + " vmax : float, optional\n", + " Maximum value for color scaling. If None, determined from the data.\n", + " scale : {'linear', 'log', 'asinh'}, optional\n", + " Scale to use for the image normalization. Default is 'asinh'.\n", + " title_prefix : str, optional\n", + " Optional prefix for the plot title.\n", + " title_path : bool, optional\n", + " If True, uses the full file path for the title;\n", + " otherwise, uses the basename. Default is False.\n", + " y_limits : tuple of float, optional\n", + " Limits for the y-axis of the 1D spectrum plot.\n", + " If None, limits are determined from the data.\n", + " is_stage3 : bool, optional\n", + " Plot stage 3 products? Default is False.\n", + "\n", + " Returns\n", + " -------\n", + " None.\n", + " \"\"\"\n", + "\n", + " # ---------------------------------Check Inputs---------------------------------\n", + " spectra = [spectra] if isinstance(spectra, str) else spectra\n", + " compare_x1d = [compare_x1d] if isinstance(compare_x1d, str) else compare_x1d\n", + " compare_mast = [compare_mast] if isinstance(compare_mast, str) else compare_mast\n", + "\n", + " # Plot stage 3 products?\n", + " if is_stage3:\n", + "\n", + " # Stage 3 products should include the source_id in the filename.\n", + " # Sort based on filename rather than open all.\n", + " def filter_prod(products, source_id):\n", + " \"\"\"Filter products based on the source_id.\"\"\"\n", + " return [f for f in products if source_id.lower() in f and\n", + " ('FXD_SLIT' not in fits.getheader(f, ext=0) or\n", + " fits.getheader(f, ext=0)['FXD_SLIT'].lower() == source_id.lower())]\n", + "\n", + " spectra = filter_prod(spectra, source_id)\n", + " compare_x1d = filter_prod(compare_x1d, source_id) if compare_x1d else None\n", + " compare_mast = filter_prod(compare_mast, source_id) if compare_mast else None\n", + "\n", + " ftypes = {ftype: [f for f in spectra\n", + " if ftype in f] for ftype in [\"cal\", \"s2d\", \"x1d\"]}\n", + " products = sorted(ftypes['s2d']) if plot_resample else sorted(ftypes['cal'])\n", + " if not products:\n", + " raise ValueError(\"No valid data products found for plotting.\")\n", + "\n", + " # --------------------------------Set up figures-------------------------------\n", + " total_plots = len(products) + bool(ftypes['x1d'])\n", + " height_ratios = [1] * len(products) + ([3] if bool(ftypes['x1d']) else [])\n", + " fig, axes = plt.subplots(total_plots, 1, figsize=(15, 5*total_plots),\n", + " sharex=False, height_ratios=height_ratios)\n", + " fig.subplots_adjust(hspace=0.2, wspace=0.2)\n", + " ax2d, ax1d = (axes[:-1], axes[-1]) if bool(ftypes['x1d']) else (axes, None)\n", + "\n", + " cmap = plt.get_cmap(cmap) # Set up colormap and bad pixel color.\n", + " cmap.set_bad(bad_color, 1.0)\n", + " colors = plt.get_cmap('tab10').colors\n", + " color_cycle = itertools.cycle(colors)\n", + "\n", + " # ---------------------------------Plot spectra--------------------------------\n", + " for i, product in enumerate(products):\n", + " model = datamodels.open(product) # Open files as JWST datamodels.\n", + "\n", + " # Extract the correct 2D source spectrum if there are multiple.\n", + " slit_m = model\n", + " if 'slits' in model:\n", + " slits = model.slits\n", + " slit_m = next((s for s in slits\n", + " if getattr(s, 'name', None) == source_id), None)\n", + " slit_m = slit_m or next((s for s in model.slits\n", + " if s.source_id == source_id), None)\n", + " if not slit_m:\n", + " print(f\"'{source_id}' not found/invalid.\")\n", + " print(f\"Available source_ids: {[s.source_id for s in slits][:5]}\")\n", + " break\n", + "\n", + " # Check if 'fixed_slit' exists, otherwise fall back to 'slitlet_id'\n", + " slit_name = (f\"SLIT: {getattr(slit_m, 'name', None) or slit_m.slitlet_id}, \"\n", + " f\"SOURCE: {getattr(slit_m, 'source_id', '')}\")\n", + "\n", + " # -----------------------Extract the 2D/3D data----------------------\n", + " data_2d = getattr(slit_m, extname)\n", + " if data_2d.ndim == 3: # Handle multi-integration data.\n", + " if integration == 'min':\n", + " data_2d = np.nanmin(data_2d, axis=0)\n", + " elif isinstance(integration, int) and 0 <= integration < data_2d.shape[0]:\n", + " data_2d = data_2d[integration]\n", + " else:\n", + " raise ValueError(f\"Invalid integration '{integration}' for 3D data.\")\n", + "\n", + " # -----------Convert from pixels to wavelength (x-axis)--------------\n", + " wcsobj = slit_m.meta.wcs # Obtaining the WCS object from the meta data.\n", + " y, x = np.mgrid[:slit_m.data.shape[0], :slit_m.data.shape[1]]\n", + " # Coordinate transform from detector space (pixels) to sky (RA, DEC).\n", + " det2sky = wcsobj.get_transform('detector', 'world')\n", + " ra, dec, s2dwave = det2sky(x, y) # RA/Dec, wavelength (microns) for each pixel.\n", + " s2dwaves = s2dwave[0, :] # Single row since this is the rectified spectrum.\n", + " x_arr = np.arange(0, slit_m.data.shape[1], int(len(slit_m.data[1]) / 4))\n", + " wav = np.round(s2dwaves[x_arr], 2) # Populating the wavelength array.\n", + " ax2d[i].set_xticks(x_arr, wav)\n", + "\n", + " # xticks = np.arange(np.ceil(wave_1d[0]), wave_1d[-1], 0.2)\n", + " # xtick_pos = np.interp(xticks, wave_1d, np.arange(num_waves))\n", + " # ax1d.set_xticks(xtick_pos)\n", + " # ax1d.set_xticklabels([f'{xtick:.1f}' for xtick in xticks])\n", + "\n", + " # ---------------------------Scale the data-------------------------\n", + " sigma_clipped_data = sigma_clip(data_2d, sigma=5, maxiters=3)\n", + " vmin = np.nanmin(sigma_clipped_data) if vmin is None else vmin\n", + " vmax = np.nanmax(sigma_clipped_data) if vmax is None else vmax\n", + " stretch_map = {'log': LogStretch(), 'linear': LinearStretch(),\n", + " 'asinh': AsinhStretch()}\n", + " if scale in stretch_map:\n", + " norm = ImageNormalize(sigma_clipped_data,\n", + " interval=ManualInterval(vmin=vmin, vmax=vmax),\n", + " stretch=stretch_map[scale])\n", + " else:\n", + " norm = simple_norm(sigma_clipped_data, vmin=vmin, vmax=vmax)\n", + "\n", + " # -------------------------Plot 1D Spectra-------------------------\n", + " for prods_1d, prefix in [(sorted(ftypes['x1d']), f'{title_prefix} '),\n", + " (compare_x1d, 'RE-EXTRACTION '),\n", + " (compare_mast, 'MAST ')]:\n", + " if prods_1d:\n", + "\n", + " model_1d = datamodels.open(prods_1d[i])\n", + " specs = model_1d.spec\n", + " spec = next((s for s in specs if\n", + " getattr(s, 'name', None) == source_id), None)\n", + " spec = spec or next((s for s in specs\n", + " if s.source_id == source_id), None)\n", + "\n", + " if spec:\n", + " tab = spec.spec_table\n", + " source_type = source_type if source_type else slit_m.source_type\n", + " wave = tab.WAVELENGTH\n", + " flux = tab.FLUX if source_type == 'POINT' else tab.SURF_BRIGHT\n", + " errs = tab.FLUX_ERROR if source_type == 'POINT' else tab.SB_ERROR\n", + "\n", + " # Expand the array to visualize the wavelength gap.\n", + " if expand_wavelength_gap:\n", + " dx1d_wave = wave[1:] - wave[:-1]\n", + " igap = np.argmax(dx1d_wave)\n", + " dx_replace = (dx1d_wave[igap-1] + dx1d_wave[igap+1]) / 2.\n", + " nfill = int(np.round(np.nanmax(dx1d_wave) / dx_replace))\n", + "\n", + " if nfill > 1:\n", + " print(f\"Expanding wavelength gap {wave[igap]:.2f} \"\n", + " f\"-- {wave[igap+1]:.2f} μm\")\n", + "\n", + " wave_fill = np.mgrid[wave[igap]:wave[igap+1]:(nfill+1)*1j]\n", + " wave = np.concatenate([wave[:igap+1],\n", + " wave_fill[1:-1],\n", + " wave[igap+1:]])\n", + "\n", + " if prefix != 'RE-EXTRACTION ':\n", + " num_rows, num_waves = data_2d.shape\n", + " fill_2d = np.zeros(shape=(num_rows, nfill-1))*np.nan\n", + " data_2d = np.concatenate([data_2d[:, :igap+1],\n", + " fill_2d, data_2d[:, igap+1:]],\n", + " axis=1)\n", + "\n", + " fill = np.zeros(shape=(nfill-1)) * np.nan\n", + " flux = np.concatenate([flux[:igap+1], fill, flux[igap+1:]])\n", + " errs = np.concatenate([errs[:igap+1], fill, errs[igap+1:]])\n", + " else:\n", + " nfill = 0\n", + "\n", + " # ----------------Construct legends and annotations-----------------\n", + " detector = slit_m.meta.instrument.detector\n", + " ffilter = slit_m.meta.instrument.filter\n", + " grating = slit_m.meta.instrument.grating\n", + " dither = model.meta.dither.position_number\n", + " label_2d = f'{grating}/{ffilter}'\n", + " label_1d = f'{detector} ({grating}/{ffilter})'\n", + " if not is_stage3:\n", + " label_2d = f'Dither/Nod {dither} ({label_2d})'\n", + " label_1d = (f'{prefix} Dither/Nod {dither} {label_1d}')\n", + " else:\n", + " label_1d = f'{prefix}{label_1d}'\n", + " ax2d[i].annotate(label_2d, xy=(1, 1), xycoords='axes fraction',\n", + " xytext=(-10, -10), textcoords='offset points',\n", + " bbox=dict(boxstyle=\"round,pad=0.3\",\n", + " edgecolor='white',\n", + " facecolor='white', alpha=0.8),\n", + " fontsize=12, ha='right', va='top')\n", + "\n", + " title_2d = (f\"{title_prefix + ' ' if title_prefix else ''}\"\n", + " f\"{model.meta.filename} | {slit_name}\")\n", + " if integration:\n", + " title_2d = title_2d.replace('.fits', f'[{integration}].fits')\n", + " ax2d[i].set_title(title_2d, fontsize=14)\n", + " if not bool(ftypes['x1d']):\n", + " ax2d[i].set_xlabel(\"Wavelength (μm)\", fontsize=12)\n", + " ax2d[i].set_ylabel(\"Pixel Row\", fontsize=12)\n", + " ax2d[i].legend(fontsize=12)\n", + "\n", + " # ------------------------------------------------------------------\n", + "\n", + " num_waves = len(wave)\n", + " color = next(color_cycle)\n", + " ax1d.step(wave, flux, lw=1, label=label_1d, color=color)\n", + " if plot_errors:\n", + " ax1d.fill_between(np.arange(num_waves), flux - errs,\n", + " flux + errs, color='grey', alpha=0.3)\n", + " ax1d.legend(fontsize=12)\n", + " ax1d.set_title(f\"{title_prefix + ' ' if title_prefix else ''}\"\n", + " f\"Extracted 1D Spectra | {slit_name}\", fontsize=14)\n", + " ax1d.set_ylabel(\"Flux (Jy)\" if source_type == 'POINT'\n", + " else \"Surface Brightness (MJy/sr)\", fontsize=12)\n", + " ax1d.set_xlabel(\"Wavelength (μm)\", fontsize=12)\n", + "\n", + " ax1d.set_ylim(y_limits or (np.nanpercentile(flux, 1),\n", + " np.nanpercentile(flux, 99.5)))\n", + "\n", + " # --------------------Plot the 2D spectra & colorbar---------------\n", + " plt.subplots_adjust(left=0.05, right=0.85)\n", + " im = ax2d[i].imshow(data_2d, origin='lower', cmap=cmap, norm=norm,\n", + " aspect=aspect, interpolation='nearest')\n", + " units = slit_m.meta.bunit_data\n", + " cbar_ax = fig.add_axes([ax2d[i].get_position().x1 + 0.02,\n", + " ax2d[i].get_position().y0, 0.02,\n", + " ax2d[i].get_position().height])\n", + " cbar = fig.colorbar(im, cax=cbar_ax)\n", + " cbar.set_label(units, fontsize=12)\n", + "\n", + " # ----------------------Add extraction region---------------------\n", + " ystart, ystop, xstart, xstop = (spec.extraction_ystart - 1,\n", + " spec.extraction_ystop - 1,\n", + " spec.extraction_xstart - 1,\n", + " spec.extraction_xstop - 1)\n", + " extract_width = ystop - ystart + 1\n", + " box = Rectangle((xstart, ystart), xstop - xstart+nfill,\n", + " extract_width, fc='None', ec=color,\n", + " lw=2, label=prefix)\n", + " ax2d[i].add_patch(box)\n", + " ax2d[i].legend()" + ] + }, + { + "cell_type": "markdown", + "id": "360c10c8-6594-4308-ab60-578c4b2614d1", + "metadata": {}, + "source": [ + "---\n", + "\n", + "### 8.1 Display `Detector1Pipeline` Products\n", + "\n", + "Inspect the Stage 1 slope products. " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "11e52fa4-eade-46f0-817c-cd69deb73334", + "metadata": {}, + "outputs": [], + "source": [ + "if doviz:\n", + " rate_file = rate_sci[-1] # Show the last rate file, as an example.\n", + " display_rate(rate_file, vmin=-0.1, vmax=1, scale='asinh',\n", + " aspect=10, title_prefix='REPROCESSED') # , extname='dq')" + ] + }, + { + "cell_type": "markdown", + "id": "3a8f322b-8f8f-46b2-906b-d5d1011d747d", + "metadata": {}, + "source": [ + "---\n", + "\n", + "### 8.2 Display `Spec2Pipeline` Products\n", + "\n", + "Inspect the Stage 2 calibrated spectra. " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "cbf1fbbf-0967-4b81-8fe6-7beb6a81c1e6", + "metadata": { + "scrolled": true + }, + "outputs": [], + "source": [ + "if doviz:\n", + " display_spectra(sci_s2d+sci_x1d, source_id='S200A1', scale='log',\n", + " vmin=-0.1e-9, vmax=3e-8, title_prefix='REPROCESSED')" + ] + }, + { + "cell_type": "markdown", + "id": "080c3014-5c7e-4b5a-a042-55679ee63c64", + "metadata": {}, + "source": [ + "---\n", + "\n", + "### 8.3 Display `Spec3Pipeline` Products\n", + "\n", + "Inspect the Stage 3 combined calibrated spectra. \n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "72b0e78a-4ff8-46fd-b4d4-cefe6a90ae2b", + "metadata": {}, + "outputs": [], + "source": [ + "# Display stage 3 products.\n", + "if doviz:\n", + " display_spectra(stage3_s2d+stage3_x1d, source_id='S200A1', scale='log',\n", + " vmin=-0.1e-9, vmax=3e-8, title_prefix='REPROCESSED', is_stage3=True)" + ] + }, + { + "cell_type": "markdown", + "id": "4d96aa47-dfda-4803-bbb0-ec21c384187a", + "metadata": {}, + "source": [ + "---\n", + "\n", + "## 9. Modifying the EXTRACT1D Reference File (as needed)\n", + "\n", + "The `extract_1d` step's `use_source_pos` parameter in Stage 2 generally centers the 1D extraction box on the actual source location effectively and thus doesn't usually require manual adjustment. However, in some cases, adjusting the position of the extraction box by modifying the EXTRACT1D reference file may be useful. The following section demonstrates how to do this.\n", + "\n", + "The EXTRACT1D reference file, along with several other parameter files, can be found in the `CRDS_PATH` directory. While some files, like `.json` files, can be manually edited, we modify them using Python.\n", + "\n", + "
\n", + " \n", + "**Warning**: Currently, there is no aperture correction in place for NIRSpec, so the `extract_width` parameter **MUST** remain unchanged (6 pixels wide; 5 for S1600A1) to ensure proper flux calibration! The extraction box limits (`ystart` and `ystop`) can be modified; however, if `ystart` and `ystop` do not match the `extract_width`, the `extract_width` takes precedence and is applied symmetrically around the midpoint between `ystart` and `ystop`.\n", + "\n", + "
\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "9181d5d6-e5f0-4a23-82ed-18ccaf9b2b23", + "metadata": {}, + "outputs": [], + "source": [ + "# Modify the EXTRACT1D reference file.\n", + "\n", + "# If you don't know the reference file name this should work.\n", + "# extract_1d_ref = Spec3Pipeline().get_reference_file(stage3_s2d, 'extract1d')\n", + "\n", + "refs = api.dump_references(crds_client.get_context_used('jwst'),\n", + " ['jwst_nirspec_extract1d_0008.json'])\n", + "extract_1d_ref = refs['jwst_nirspec_extract1d_0008.json']\n", + "\n", + "# Open EXTRACT1D reference file in read-mode.\n", + "with open(extract_1d_ref, \"r\") as ref_file:\n", + " params = json.load(ref_file)\n", + "\n", + " yshift = -4 # Applied shift in pixels as example.\n", + "\n", + " # S200A1\n", + " params[\"apertures\"][0][\"extract_width\"] = 6\n", + " params[\"apertures\"][0][\"ystart\"] += yshift\n", + " params[\"apertures\"][0][\"ystop\"] += yshift\n", + "\n", + " # S200B1\n", + " params[\"apertures\"][1][\"extract_width\"] = 6\n", + " params[\"apertures\"][1][\"ystart\"] = 26.5\n", + " params[\"apertures\"][1][\"ystop\"] = 31.5\n", + "\n", + " # S200A2\n", + " params[\"apertures\"][2][\"extract_width\"] = 6\n", + " params[\"apertures\"][2][\"ystart\"] = 26.5\n", + " params[\"apertures\"][2][\"ystop\"] = 31.5\n", + "\n", + " # S400A1\n", + " params[\"apertures\"][3][\"extract_width\"] = 6\n", + " params[\"apertures\"][3][\"ystart\"] = 31\n", + " params[\"apertures\"][3][\"ystop\"] = 36\n", + "\n", + " # S1600A1\n", + " params[\"apertures\"][4][\"extract_width\"] = 5\n", + " params[\"apertures\"][4][\"ystart\"] = 14\n", + " params[\"apertures\"][4][\"ystop\"] = 18\n", + "\n", + "# Write changes to a new file.\n", + "newData = json.dumps(params, indent=4)\n", + "# Add the suffix '_fs' to distinguish the file from the default version.\n", + "basename = os.path.basename(extract_1d_ref)[:-5]\n", + "extract_1d_ref_mod = os.path.join(basedir, basename + \"_fs.json\")\n", + "with open(extract_1d_ref_mod, \"w\") as file:\n", + " file.write(newData)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "86d136e3-5dbd-4a20-9d50-0925aaee61d5", + "metadata": { + "scrolled": true + }, + "outputs": [], + "source": [ + "# Inspect the EXTRACT1D reference file.\n", + "with open(extract_1d_ref_mod, 'r') as f_obj:\n", + " extract_1d_ref_mod_data = json.load(f_obj)\n", + "\n", + "JSON(extract_1d_ref_mod_data, expanded=True)" + ] + }, + { + "cell_type": "markdown", + "id": "c1cab3dc-25f7-41b7-9d76-734811231625", + "metadata": {}, + "source": [ + "Now, we re-extract the 1D spectrum by running the `Extract1dStep` and overriding the reference file.\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "8928c904-9a7b-4686-9f19-2ef6937b9888", + "metadata": {}, + "outputs": [], + "source": [ + "Extract1dStep.call(stage3_s2d,\n", + " save_results=True,\n", + " output_dir=spec3_dir,\n", + " output_use_model=True,\n", + " suffix='x1d_mod', # Change suffix to easily find modified file.\n", + " use_source_posn=False,\n", + " override_extract1d=extract_1d_ref_mod)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "0c86652f-3d0d-44dc-a804-28a4aa7e90a4", + "metadata": {}, + "outputs": [], + "source": [ + "stage3_x1ds_mod = sorted(glob.glob(spec3_dir + '*_x1d_mod.fits'))\n", + "display_spectra(stage3_s2d+stage3_x1d, compare_x1d=stage3_x1ds_mod, source_id='S200A1',\n", + " scale='log', vmin=-0.1e-9, vmax=3e-8,\n", + " title_prefix='REPROCESSED', is_stage3=True)" + ] + }, + { + "cell_type": "markdown", + "id": "ddf937b0-34ce-404a-86bf-0d5e55968ef0", + "metadata": {}, + "source": [ + "As expected, the spectrum extracted in the shifted location has lower flux that the spectrum extracted in the center of the 2D spectral trace. \n", + "\n", + "---\n", + "\n", + "
\n", + " \"Space\n", + "
\n", + " \n", + "[Top of Page](#NIRSpec-FS-Pipeline-Notebook)" + ] + } + ], + "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.11.10" + }, + "widgets": { + "application/vnd.jupyter.widget-state+json": { + "state": {}, + "version_major": 2, + "version_minor": 0 + } + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/notebooks/NIRSPEC/FSlit/requirements.txt b/notebooks/NIRSPEC/FSlit/requirements.txt new file mode 100644 index 0000000..2b5b1b0 --- /dev/null +++ b/notebooks/NIRSPEC/FSlit/requirements.txt @@ -0,0 +1,7 @@ +astropy>=6.1.4 +astroquery>=0.4.8.dev9474 +jwst==1.16.0 +matplotlib>=3.9.2 +numpy==1.26.4 +requests>=2.32.3 + diff --git a/notebooks/NIRSPEC/requirements.txt b/notebooks/NIRSPEC/requirements.txt deleted file mode 100644 index 6dd201d..0000000 --- a/notebooks/NIRSPEC/requirements.txt +++ /dev/null @@ -1,3 +0,0 @@ -numpy -crds -jwst>=1.13.3