From a90fc66b20f408ac9e3ca24771dbe8e4d6843cc2 Mon Sep 17 00:00:00 2001 From: Hatice Karatay Date: Mon, 15 Apr 2024 16:43:38 -0400 Subject: [PATCH] Consolidate convenience functions into a script --- .../NIRSpec_NSClean/FS_NSClean_example.ipynb | 830 ++++++------------ .../NIRSpec_NSClean/IFU_NSClean_example.ipynb | 610 +------------ .../NIRSpec_NSClean/MOS_NSClean_example.ipynb | 610 +------------ notebooks/NIRSpec_NSClean/utils.py | 525 +++++++++++ 4 files changed, 822 insertions(+), 1753 deletions(-) create mode 100644 notebooks/NIRSpec_NSClean/utils.py diff --git a/notebooks/NIRSpec_NSClean/FS_NSClean_example.ipynb b/notebooks/NIRSpec_NSClean/FS_NSClean_example.ipynb index cbe8b0756..97a2e0e41 100644 --- a/notebooks/NIRSpec_NSClean/FS_NSClean_example.ipynb +++ b/notebooks/NIRSpec_NSClean/FS_NSClean_example.ipynb @@ -19,19 +19,18 @@ "\n", "* 1. [Introduction](#introduction)\n", "* 2. [Import Library](#imports)\n", - "* 3. [Convenience Functions](#functions)\n", - "* 4. [Download the Data](#data)\n", - "* 5. [Running `Spec2Pipeline` without NSClean (Original Data)](#nsclean_skipped)\n", - "* 6. [Clean up 1/f Noise with NSClean (Default Pipeline Mask)](#nsclean_default)\n", - " * 6.1 [Verify the Mask (Default Pipeline Mask)](#verify_default_mask)\n", - " * 6.2 [Comparing Original vs. Cleaned Data (Default Pipeline Mask)](#nsclean_default_compare)\n", - "* 7. [Clean up 1/f Noise with NSClean (Alternate Mask)](#nsclean_alternate)\n", - " * 7.1 [Verify the Mask (Alternate Mask)](#verify_alternate_mask)\n", - " * 7.2 [Comparing Original vs. Cleaned Data (Alternate Mask)](#nsclean_alternate_compare)\n", - "* 8. [Clean up 1/f Noise with NSClean (Hand-Modified Mask)](#nsclean_modified)\n", - " * 8.1 [Verify the Mask (Hand-Modified Mask)](#verify_modified_mask)\n", - " * 8.2 [Comparing Original vs. Cleaned Data (Hand-Modified Mask)](#nsclean_modified_compare)\n", - "* 9. [Conclusion](#conclusion)\n", + "* 3. [Download the Data](#data)\n", + "* 4. [Running `Spec2Pipeline` without NSClean (Original Data)](#nsclean_skipped)\n", + "* 5. [Clean up 1/f Noise with NSClean (Default Pipeline Mask)](#nsclean_default)\n", + " * 5.1 [Verify the Mask (Default Pipeline Mask)](#verify_default_mask)\n", + " * 5.2 [Comparing Original vs. Cleaned Data (Default Pipeline Mask)](#nsclean_default_compare)\n", + "* 6. [Clean up 1/f Noise with NSClean (Alternate Mask)](#nsclean_alternate)\n", + " * 6.1 [Verify the Mask (Alternate Mask)](#verify_alternate_mask)\n", + " * 6.2 [Comparing Original vs. Cleaned Data (Alternate Mask)](#nsclean_alternate_compare)\n", + "* 7. [Clean up 1/f Noise with NSClean (Hand-Modified Mask)](#nsclean_modified)\n", + " * 7.1 [Verify the Mask (Hand-Modified Mask)](#verify_modified_mask)\n", + " * 7.2 [Comparing Original vs. Cleaned Data (Hand-Modified Mask)](#nsclean_modified_compare)\n", + "* 8. [Conclusion](#conclusion)\n", "\n", "\n", "\n", @@ -89,584 +88,27 @@ "import time as tt\n", "import logging\n", "import warnings\n", - "import requests\n", "\n", "# ------ JWST Calibration Pipeline Imports ------\n", "from jwst.pipeline.calwebb_spec2 import Spec2Pipeline\n", "\n", - "# ------ Plotting/Stats Imports ------\n", + "# # ------ Plotting/Stats Imports ------\n", "from matplotlib import pyplot as plt\n", - "from astropy.stats import sigma_clip\n", "from astropy.io import fits\n", "\n", + "from utils import get_jwst_file, plot_dark_data, plot_cleaned_data, plot_spectra\n", + "\n", "# Hide all log and warning messages.\n", "logging.disable(logging.ERROR)\n", "warnings.simplefilter(\"ignore\", RuntimeWarning)" ] }, - { - "cell_type": "markdown", - "id": "574fab53-aa67-434c-bcbe-17360421e3b4", - "metadata": {}, - "source": [ - "## 3. Convenience Functions \n", - "
" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "71fb474f-e349-43ed-b3e5-b0ac46fc923e", - "metadata": {}, - "outputs": [], - "source": [ - "def get_jwst_file(name, mast_api_token=None, save_directory=\".\", redownload=False):\n", - " \"\"\"\n", - " Retrieve a JWST data file from MAST archive and save it to a specified directory.\n", - "\n", - " Parameters:\n", - " ----------\n", - " name: str\n", - " File name.\n", - " mast_api_token: str\n", - " MAST authorization token. Get your MAST Token Here: https://auth.mast.stsci.edu/token.\n", - " save_directory: str\n", - " Save directory path.\n", - " redownload: bool\n", - " Redownload the data even if it exsits already?\n", - " \"\"\"\n", - " mast_url = \"https://mast.stsci.edu/api/v0.1/Download/file\"\n", - " params = dict(uri=f\"mast:JWST/product/{name}\")\n", - "\n", - " if mast_api_token:\n", - " headers = dict(Authorization=f\"token {mast_api_token}\")\n", - " else:\n", - " headers = {}\n", - "\n", - " file_path = os.path.join(save_directory, name)\n", - "\n", - " # Check if the file already exists in the save directory.\n", - " if os.path.exists(file_path) and not redownload:\n", - " print(f\"The file {name} already exists in the directory. Skipping download.\")\n", - " return\n", - "\n", - " r = requests.get(mast_url, params=params, headers=headers, stream=True)\n", - " r.raise_for_status()\n", - "\n", - " with open(file_path, \"wb\") as fobj:\n", - " for chunk in r.iter_content(chunk_size=1024000):\n", - " fobj.write(chunk)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "2c996d1b-9d98-4f40-b254-f1daa252926c", - "metadata": {}, - "outputs": [], - "source": [ - "# Define a helper function for plotting the dark areas (the mask).\n", - "\n", - "def plot_dark_data(\n", - " rate_file,\n", - " mask_file,\n", - " cmap=\"viridis\",\n", - " scale=0.2,\n", - " slice_index=1,\n", - " axis=0,\n", - " aspect=\"auto\",\n", - " layout=\"rows\",\n", - "):\n", - " \"\"\"\n", - " Plot the dark areas on the detector, masking out the illuminated areas.\n", - " This function can take 2D (_rate.fits) and 3D (_rateints.fits) files as input.\n", - "\n", - " Parameters:\n", - " ----------\n", - " rate_file: str\n", - " File path to the FITS file containing countrate data.\n", - " mask_file: str\n", - " File path to the FITS file containing mask data.\n", - " slice_index: int\n", - " Index of the slice to be plotted, 2D data default is 1.\n", - " axis: int (0-2)\n", - " Axis along which the slice will be taken\n", - " (0 for x-axis, 1 for y-axis, 2 for z-axis).\n", - " aspect: int or 'auto'\n", - " Plot's aspect ratio.\n", - " cmap: str\n", - " Color map.\n", - " layout: str\n", - " Layout subplots in rows or columns?\n", - " scale: float\n", - " Scaling factor applied to determine the intensity range for visualization.\n", - " \"\"\"\n", - "\n", - " # Make a matplotlib figure.\n", - " if layout == \"rows\":\n", - " fig, ax = plt.subplots(3, 1, figsize=(10, 5))\n", - " else:\n", - " fig, ax = plt.subplots(1, 3, figsize=(20, 5))\n", - " plt.suptitle(f\"Dark Areas for {os.path.basename(rate_file)}\", fontsize=15)\n", - "\n", - " # Open the mask and science data.\n", - " with fits.open(mask_file) as hdulist:\n", - " mask = hdulist[1].data\n", - "\n", - " with fits.open(rate_file) as hdulist:\n", - " sci = hdulist[\"SCI\"].data\n", - "\n", - " # Determine if countrate data is 2D or 3D.\n", - " if len(sci.shape) == 3:\n", - " dim = 3\n", - " else:\n", - " dim = 2\n", - "\n", - " # Get data limits from the dark data.\n", - " masked_sci = sci.copy()\n", - " masked_sci[mask == 0] = 0\n", - " vmin = np.nanpercentile(masked_sci, scale)\n", - " vmax = np.nanpercentile(masked_sci, 100 - scale)\n", - "\n", - " # Plot the science image with limits from the dark data.\n", - " sci[np.isnan(sci)] = 0\n", - "\n", - " # If the countrate data is 3D, determine the integration (slice) to plot.\n", - " if dim == 3:\n", - " sci = np.take(sci, slice_index, axis=axis)\n", - " mask = np.take(mask, slice_index, axis=axis)\n", - " masked_sci = np.take(masked_sci, slice_index, axis=axis)\n", - "\n", - " ax[0].set_title(\"Original Rate Data: Integration [{},:,:]\".format(slice_index))\n", - " ax[1].set_title(\n", - " \"Illuminated Pixel Mask: Integration [{},:,:]\".format(slice_index)\n", - " )\n", - " ax[2].set_title(\n", - " \"Dark Data (Illuminated Pixels Masked): Integration [{},:,:]\".format(\n", - " slice_index\n", - " )\n", - " )\n", - "\n", - " else:\n", - " sci = sci\n", - " ax[0].set_title(\"Original Rate Data\")\n", - " ax[1].set_title(\"Illuminated Pixel Mask\")\n", - " ax[2].set_title(\"Dark Data (Illuminated Pixels Masked)\")\n", - "\n", - " # Plot the science image with limits from the dark data.\n", - " ax[0].set_ylabel(\"Pixel Row\", fontsize=12)\n", - " ax[0].set_xlabel(\"Pixel Column\", fontsize=12)\n", - " im1 = ax[0].imshow(\n", - " sci, cmap=cmap, origin=\"lower\", aspect=aspect, vmin=vmin, vmax=vmax\n", - " )\n", - " cbar1 = fig.colorbar(im1, ax=ax[0], pad=0.05, shrink=0.7, label=\"DN/s\")\n", - "\n", - " # Plot the mask: values are 1 or 0.\n", - " ax[1].set_ylabel(\"Pixel Row\", fontsize=12)\n", - " ax[1].set_xlabel(\"Pixel Column\", fontsize=12)\n", - " im2 = ax[1].imshow(mask, cmap=cmap, aspect=aspect, origin=\"lower\", vmin=0, vmax=1)\n", - " cbar2 = fig.colorbar(im2, ax=ax[1], pad=0.05, shrink=0.7, ticks=[0, 1])\n", - " cbar2.set_ticklabels([\"Illuminated Pixel Masked\", \"Un-Illuminated Pixel\"])\n", - "\n", - " # Plot the dark data with the same limits as the science data.\n", - " ax[2].set_ylabel(\"Pixel Row\", fontsize=12)\n", - " ax[2].set_xlabel(\"Pixel Column\", fontsize=12)\n", - " im3 = ax[2].imshow(\n", - " masked_sci, cmap=cmap, aspect=aspect, origin=\"lower\", vmin=vmin, vmax=vmax\n", - " )\n", - " cbar3 = fig.colorbar(im3, ax=ax[2], pad=0.05, shrink=0.7, label=\"DN/s\")\n", - "\n", - " if layout == \"rows\": # Adjust the multiplier as needed.\n", - " cbar1.ax.figure.set_size_inches(cbar1.ax.figure.get_size_inches() * 1.2)\n", - " cbar2.ax.figure.set_size_inches(cbar2.ax.figure.get_size_inches() * 1.2)\n", - " cbar3.ax.figure.set_size_inches(cbar3.ax.figure.get_size_inches() * 1.2)\n", - "\n", - " fig.tight_layout() # Adjusted tight_layout for better suptitle spacing." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "3a209311-2eba-45e4-b450-0c89d67b5f70", - "metadata": {}, - "outputs": [], - "source": [ - "# Define a helper function for plotting the cleaned data.\n", - "\n", - "def plot_cleaned_data(\n", - " rate_file,\n", - " cleaned_file,\n", - " slice_indx=1,\n", - " aspect=\"auto\",\n", - " cmap=\"viridis\",\n", - " scale=0.2,\n", - " layout=\"rows\",\n", - "):\n", - " \"\"\"\n", - " Plot the 2D cleaned rate data (or 2D slice from a 3D cleaned data cube)\n", - " and compare against the original data (not applying NSClean).\n", - "\n", - " Parameters:\n", - " ----------\n", - " rate_file: str\n", - " File path to the FITS file containing original countrate data.\n", - " cleaned_file: str\n", - " File path to the FITS file containing cleaned countrate data.\n", - " slice_indx: int\n", - " Index of the slice to be plotted, 2D data default is 1.\n", - " aspect: int or 'auto'\n", - " Plot's aspect ratio.\n", - " cmap: str\n", - " Color map.\n", - " layout: str\n", - " Layout subplots in rows or columns?\n", - " scale: float\n", - " Scaling factor applied to determine the intensity range for visualization.\n", - " \"\"\"\n", - "\n", - " # Make a matplotlib figure.\n", - " if layout == \"rows\":\n", - " fig, ax = plt.subplots(4, 1, figsize=(10, 12))\n", - " else:\n", - " fig, ax = plt.subplots(2, 2, figsize=(12, 10))\n", - " ax = ax.flatten()\n", - " plt.suptitle(\n", - " f\"Cleaned Data (1/f noise removed) for {os.path.basename(rate_file)}\",\n", - " fontsize=15,\n", - " y=1,\n", - " )\n", - "\n", - " # Open the original and cleaned data.\n", - " with fits.open(rate_file) as hdulist:\n", - " original = hdulist[\"SCI\"].data\n", - "\n", - " with fits.open(cleaned_file) as hdulist:\n", - " cleaned = hdulist[\"SCI\"].data\n", - " # Determine if rate data is 2D or 3D.\n", - " if len(cleaned.shape) == 3:\n", - " dim = 3\n", - " else:\n", - " dim = 2\n", - "\n", - " # Define image limits from the original data.\n", - " vmin = np.nanpercentile(original, scale)\n", - " vmax = np.nanpercentile(original, 100 - scale)\n", - "\n", - " original[np.isnan(original)] = 0\n", - " cleaned[np.isnan(cleaned)] = 0\n", - "\n", - " # If the rate data is 3D, determine the integration (slice) to plot.\n", - " if dim == 3:\n", - " original = original[slice_indx, :, :]\n", - " cleaned = cleaned[slice_indx, :, :]\n", - " ax[0].set_title(\n", - " \"Original Rate Data: Integration [{},:,:]\".format(slice_indx), fontsize=15\n", - " )\n", - " ax[1].set_title(\n", - " \"Cleaned Rate Data: Integration [{},:,:]\".format(slice_indx), fontsize=15\n", - " )\n", - "\n", - " else:\n", - " ax[0].set_title(\"Original Rate Data\", fontsize=12)\n", - " ax[1].set_title(\"Cleaned Rate Data\", fontsize=12)\n", - "\n", - " # Plot the original rate data.\n", - " ax[0].set_xlabel(\"Pixel Column\", fontsize=10)\n", - " ax[0].set_ylabel(\"Pixel Row\", fontsize=10)\n", - " fig.colorbar(\n", - " ax[0].imshow(\n", - " original, cmap=cmap, origin=\"lower\", aspect=aspect, vmin=vmin, vmax=vmax\n", - " ),\n", - " ax=ax[0],\n", - " pad=0.05,\n", - " shrink=0.7,\n", - " label=\"DN/s\",\n", - " )\n", - "\n", - " # Plot the cleaned data with the same image limits.\n", - " ax[1].set_xlabel(\"Pixel Column\", fontsize=10)\n", - " ax[1].set_ylabel(\"Pixel Row\", fontsize=10)\n", - " fig.colorbar(\n", - " ax[1].imshow(\n", - " cleaned, cmap=cmap, origin=\"lower\", aspect=aspect, vmin=vmin, vmax=vmax\n", - " ),\n", - " ax=ax[1],\n", - " pad=0.05,\n", - " shrink=0.7,\n", - " label=\"DN/s\",\n", - " )\n", - "\n", - " # Plot the relative difference between the original and cleaned data.\n", - " diff = (original - cleaned) / original\n", - " diff[~np.isfinite(diff)] = 0\n", - " vmin_diff = np.nanpercentile(diff, scale)\n", - " vmax_diff = np.nanpercentile(diff, 100 - scale)\n", - " ax[2].set_title(\"Relative Difference (Original - Cleaned Rate Data)\", fontsize=12)\n", - " ax[2].set_xlabel(\"Pixel Column\", fontsize=10)\n", - " ax[2].set_ylabel(\"Pixel Row\", fontsize=10)\n", - " fig.colorbar(\n", - " ax[2].imshow(\n", - " diff,\n", - " cmap=cmap,\n", - " origin=\"lower\",\n", - " aspect=aspect,\n", - " vmin=vmin_diff,\n", - " vmax=vmax_diff,\n", - " ),\n", - " ax=ax[2],\n", - " pad=0.05,\n", - " shrink=0.7,\n", - " label=\"DN/s\",\n", - " )\n", - "\n", - " # Plot the relative difference again,\n", - " # this time hiding the outliers so that low-level\n", - " # background changes can be seen.\n", - " hide_outliers = np.ma.filled(sigma_clip(diff, masked=True), fill_value=0)\n", - " vmin_outliers = np.nanpercentile(hide_outliers, scale)\n", - " vmax_outliers = np.nanpercentile(hide_outliers, 100 - scale)\n", - " ax[3].set_title(\n", - " \"Relative Difference (Original - Cleaned Rate Data) \\n with Outliers Hidden\",\n", - " fontsize=12,\n", - " )\n", - " ax[3].set_xlabel(\"Pixel Column\", fontsize=10)\n", - " ax[3].set_ylabel(\"Pixel Row\", fontsize=10)\n", - " fig.colorbar(\n", - " ax[3].imshow(\n", - " hide_outliers,\n", - " cmap=cmap,\n", - " origin=\"lower\",\n", - " aspect=aspect,\n", - " vmin=vmin_outliers,\n", - " vmax=vmax_outliers,\n", - " ),\n", - " ax=ax[3],\n", - " pad=0.05,\n", - " shrink=0.7,\n", - " label=\"DN/s\",\n", - " )\n", - "\n", - " fig.tight_layout()" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "6662a39d-6a4b-49c1-b010-25908fa53fa3", - "metadata": {}, - "outputs": [], - "source": [ - "# Define a helper function for plotting the 1D spectra.\n", - "\n", - "def plot_spectra(\n", - " spec_list,\n", - " ext_num=1,\n", - " scale_percent=2.0,\n", - " wavelength_range=None,\n", - " flux_range=None,\n", - " xlim_low=None,\n", - "):\n", - " \"\"\"\n", - " Plot the cleaned extracted 1D spectrum and compare against\n", - " the original (not applying NSClean) extracted 1D spectrum.\n", - "\n", - " Parameters:\n", - " ----------\n", - " spec_list: list\n", - " List of paths to the FITS files containing original/cleaned 1D extracted data.\n", - " scale_percent: float\n", - " Scaling factor applied to determine the intensity range for visualization.\n", - " ext_num: int\n", - " Index/extension of the slice to be plotted.\n", - " The EXTVER header value. The default is 1 for 2D data.\n", - " wavelength_range: dict\n", - " Wavelength range (x-axis) {'nrs1': [3.6, 3.65], 'nrs2': [1.65, 1.75]}.\n", - " flux_range: dict\n", - " Flux range (y-axis) {'nrs1': [1, 2], 'nrs2': [1, 2]}.\n", - " xlim_low: int\n", - " Define a lower wavelength end for the 1D spectrum. Helpful for BOTS data.\n", - " \"\"\"\n", - "\n", - " if wavelength_range is None:\n", - " wavelength_range = {}\n", - "\n", - " # Open the FITS files.\n", - " original_hdul = fits.open(spec_list[0])\n", - " cleaned_hdul = fits.open(spec_list[1])\n", - " if len(spec_list) == 3:\n", - " alternate_hdul = fits.open(spec_list[2])\n", - " elif len(spec_list) == 4:\n", - " alternate_hdul = fits.open(spec_list[2])\n", - " handmod_hdul = fits.open(spec_list[3])\n", - " else:\n", - " alternate_hdul = None\n", - " handmod_hdul = None\n", - "\n", - " # Find the spectral extension (EXTRACT1D).\n", - " for extnum in range(len(original_hdul)):\n", - " hdu = original_hdul[extnum]\n", - " if hdu.name != \"EXTRACT1D\":\n", - " continue\n", - " slit_name = hdu.header[\"SLTNAME\"]\n", - "\n", - " if hdu.header[\"EXTVER\"] == ext_num:\n", - "\n", - " # Plot the original and cleaned spectra together.\n", - " fig, ax = plt.subplots(1, 2, figsize=(10, 5))\n", - " plt.suptitle(\n", - " f\"Compare 1D Spectra for {os.path.basename(spec_list[0])};\"\n", - " f\"EXP_TYPE/Slit = {slit_name}\",\n", - " fontsize=15,\n", - " )\n", - "\n", - " if \"nrs1\" in spec_list[0] and xlim_low is not None:\n", - " xlim_low = xlim_low\n", - " else:\n", - " xlim_low = 0\n", - "\n", - " ax[0].step(\n", - " hdu.data[\"WAVELENGTH\"][xlim_low:],\n", - " hdu.data[\"FLUX\"][xlim_low:],\n", - " linewidth=1,\n", - " label=\"Original\",\n", - " )\n", - " ax[0].step(\n", - " cleaned_hdul[extnum].data[\"WAVELENGTH\"][xlim_low:],\n", - " cleaned_hdul[extnum].data[\"FLUX\"][xlim_low:],\n", - " linewidth=1,\n", - " label=\"Cleaned\",\n", - " )\n", - " ax[0].set_xlabel(f\"Wavelength ({hdu.header['TUNIT1']})\", fontsize=12)\n", - " ax[0].set_ylabel(f\"Flux ({hdu.header['TUNIT2']})\", fontsize=12)\n", - " ax[0].set_title(\"1D Extracted Spectra\", fontsize=12)\n", - "\n", - " # Plot the difference between the spectra as a ratio.\n", - " diff = (\n", - " cleaned_hdul[extnum].data[\"FLUX\"][xlim_low:]\n", - " / hdu.data[\"FLUX\"][xlim_low:]\n", - " )\n", - " ax[1].step(\n", - " hdu.data[\"WAVELENGTH\"][xlim_low:],\n", - " diff,\n", - " linewidth=1,\n", - " label=\"Cleaned/Original\",\n", - " )\n", - " ax[1].set_xlabel(f\"Wavelength ({hdu.header['TUNIT1']})\", fontsize=12)\n", - " ax[1].set_ylabel(\"Cleaned/Original\", fontsize=12)\n", - " ax[1].set_title(\"Difference After Cleaning\", fontsize=12)\n", - " # print(\"Average Difference = {}\".format(np.nanmean(diff)))\n", - "\n", - " # If available, also plot the alternate spectra.\n", - " if alternate_hdul is not None:\n", - " ax[0].step(\n", - " alternate_hdul[extnum].data[\"WAVELENGTH\"][xlim_low:],\n", - " alternate_hdul[extnum].data[\"FLUX\"][xlim_low:],\n", - " linewidth=1,\n", - " label=\"Cleaned (alternate mask)\",\n", - " )\n", - " diff2 = (\n", - " alternate_hdul[extnum].data[\"FLUX\"][xlim_low:]\n", - " / hdu.data[\"FLUX\"][xlim_low:]\n", - " )\n", - " ax[1].step(\n", - " hdu.data[\"WAVELENGTH\"][xlim_low:],\n", - " diff2,\n", - " linewidth=1,\n", - " label=\"Cleaned (alternate mask)/Original\",\n", - " )\n", - " if handmod_hdul is not None:\n", - "\n", - " ax[0].step(\n", - " handmod_hdul[extnum].data[\"WAVELENGTH\"][xlim_low:],\n", - " handmod_hdul[extnum].data[\"FLUX\"][xlim_low:],\n", - " linewidth=1,\n", - " label=\"Cleaned (hand-modified mask)\",\n", - " color=\"Purple\"\n", - " )\n", - " diff3 = (\n", - " handmod_hdul[extnum].data[\"FLUX\"][xlim_low:]\n", - " / hdu.data[\"FLUX\"][xlim_low:]\n", - " )\n", - " ax[1].step(\n", - " hdu.data[\"WAVELENGTH\"][xlim_low:],\n", - " diff3,\n", - " linewidth=1,\n", - " label=\"Cleaned (hand-modified mask)/Original\",\n", - " )\n", - "\n", - " # Set the y-range of the plot if needed.\n", - " if flux_range is not None:\n", - " for key, y_range in flux_range.items():\n", - " if key.lower() in spec_list[0].lower() and y_range is not None:\n", - " if len(y_range) == 2:\n", - " ax[0].set_ylim(y_range)\n", - " ax[1].set_ylim(\n", - " [\n", - " np.nanpercentile(diff, scale_percent),\n", - " np.nanpercentile(diff, 100 - scale_percent),\n", - " ]\n", - " )\n", - " # ax[1].set_ylim(0.5, 1.5)\n", - "\n", - " else:\n", - " all_flux = [\n", - " hdu.data[\"FLUX\"],\n", - " cleaned_hdul[extnum].data[\"FLUX\"],\n", - " ]\n", - " y_range_ax0 = [\n", - " np.nanpercentile(all_flux[0], scale_percent),\n", - " np.nanpercentile(all_flux[0], 100 - scale_percent),\n", - " ]\n", - " ax[0].set_ylim(y_range_ax0)\n", - " # ax[1].set_ylim(0.5, 1.5)\n", - "\n", - " ax[1].set_ylim(\n", - " [\n", - " np.nanpercentile(diff, scale_percent),\n", - " np.nanpercentile(diff, 100 - scale_percent),\n", - " ]\n", - " )\n", - "\n", - " else:\n", - " all_flux = [hdu.data[\"FLUX\"], cleaned_hdul[extnum].data[\"FLUX\"]]\n", - " y_range_ax0 = [\n", - " np.nanpercentile(all_flux[0], scale_percent),\n", - " np.nanpercentile(all_flux[0], 100 - scale_percent),\n", - " ]\n", - " ax[0].set_ylim(y_range_ax0)\n", - " # ax[1].set_ylim(0.5, 1.5)\n", - "\n", - " ax[1].set_ylim(\n", - " [\n", - " np.nanpercentile(diff, scale_percent),\n", - " np.nanpercentile(diff, 100 - scale_percent),\n", - " ]\n", - " )\n", - "\n", - " # Set the x-range of the plot if needed.\n", - " for key in wavelength_range:\n", - " if key in spec_list[0] and wavelength_range[key] is not None:\n", - " ax[0].set_xlim(wavelength_range[key])\n", - " ax[1].set_xlim(wavelength_range[key])\n", - "\n", - " fig.tight_layout()\n", - "\n", - " ax[0].legend()\n", - " ax[1].legend()\n", - "\n", - " original_hdul.close()\n", - " cleaned_hdul.close()\n", - " if alternate_hdul is not None:\n", - " alternate_hdul.close()\n", - " handmod_hdul.close()" - ] - }, { "cell_type": "markdown", "id": "de1354d4-d360-4485-a466-4d459510a128", "metadata": {}, "source": [ - "## 4. Download the Data \n", + "## 3. Download the Data \n", "
\n", "\n", "The input data for this notebook features fixed slit (FS) observations with the S200A1 subarray (SUBS200A1) and full-frame observations with S200A1 as the primary slit. The FS observation of F dwarf GSC 03162-00665 with subarray SUBS200A1 and grating/filter G395M/F290LP is part of the GTO program [2757](https://www.stsci.edu/jwst/science-execution/program-information?id=2757), specifically observation 2. It consists of 15 integrations (5 integration/exposure; 3-POINT-NOD) with 4 groups each. The FS observation of brown dwarf J03480772-6022270 with full-frame readout and grating/filter G140M/F100LP is part of the GTO program [1189](https://www.stsci.edu/jwst/science-execution/program-information?id=1189), specifically observation 1. It consists of 3 integrations (3-POINT-NOD) with 11 groups each. \n", @@ -727,7 +169,7 @@ "id": "803df294-7a40-4cf6-bb2c-b1796d875979", "metadata": {}, "source": [ - "## 5. Running `Spec2Pipeline` without NSClean (Original Data) \n", + "## 4. Running `Spec2Pipeline` without NSClean (Original Data) \n", "
\n", "\n", "The cell below executes the `Spec2Pipeline`, explicitly skipping the NSClean step during processing. The level 2 products generated will serve as a reference point to illustrate how the countrate images and final extracted spectra appear without the removal of 1/f noise." @@ -784,7 +226,7 @@ "id": "476e2326-4895-4851-b1bd-7e9f6bc443cc", "metadata": {}, "source": [ - "## 6. Clean up 1/f Noise with NSClean (Default Pipeline Mask) \n", + "## 5. Clean up 1/f Noise with NSClean (Default Pipeline Mask) \n", "
\n", "\n", "If a user-supplied mask file is not provided to the [NSClean step](https://jwst-pipeline.readthedocs.io/en/latest/jwst/nsclean/index.html) in the `Spec2Pipeline`, the pipeline will generate a mask based on default parameters. This mask will identify any pixel that is unilluminated. That is, the mask must contain True and False values, where True indicates that the pixel is dark, and False indicates that the pixel is illuminated (not dark).\n", @@ -870,7 +312,7 @@ "id": "0562b0de-fff0-4420-982f-50898ad61f58", "metadata": {}, "source": [ - "### 6.1 Verify the Mask (Default Pipeline Mask) \n", + "### 5.1 Verify the Mask (Default Pipeline Mask) \n", "
\n", "\n", "Check the mask against the rate data to make sure it keeps only dark areas of the detector.\n", @@ -911,7 +353,7 @@ "id": "4fe9880c-db9f-4d19-9a10-5a81d9f29d28", "metadata": {}, "source": [ - "### 6.2 Comparing Original vs. Cleaned Data (Default Pipeline Mask) \n", + "### 5.2 Comparing Original vs. Cleaned Data (Default Pipeline Mask) \n", "
\n", "\n", "We can now compare the cleaned data (with the default pipeline mask) to the original rate file and verify that the 1/f noise has been reduced.\n", @@ -1054,7 +496,7 @@ "id": "79f8833f-fabf-482c-8edb-7b37413938e7", "metadata": {}, "source": [ - "## 7. Clean up 1/f Noise with NSClean (Alternate Mask) \n", + "## 6. Clean up 1/f Noise with NSClean (Alternate Mask) \n", "
\n", "\n", "For some data sets, masking the entire science region may excessively mask dark areas of the detector that could be used to improve the background fit. Excessive masking can introduce some high frequency noise in the cleaning process that appears as vertical striping over the spectral traces. For example, see the cleaned rate data for the subarray (PID 2757) above in section 6.2. Since this region is completely masked, some residual artifacts are introduced. \n", @@ -1128,7 +570,7 @@ "id": "3e23569f-f004-4bc6-9588-04769bc1744f", "metadata": {}, "source": [ - "### 7.1 Verify the Mask (Alternate Mask) \n", + "### 6.1 Verify the Mask (Alternate Mask) \n", "
\n", "\n", "Check the mask against the rate data to make sure it keeps only dark areas of the detector.\n" @@ -1165,7 +607,7 @@ "id": "f57c16a1-5561-4b96-97e7-9c9dc3591531", "metadata": {}, "source": [ - "### 7.2 Comparing Original vs. Cleaned Data (Alternate Mask) \n", + "### 6.2 Comparing Original vs. Cleaned Data (Alternate Mask) \n", "
" ] }, @@ -1292,7 +734,7 @@ "id": "fc01eb3d-4eeb-48d2-bb87-ffcee182d074", "metadata": {}, "source": [ - "## 8. Clean up 1/f Noise with NSClean (Hand-Modified Mask) \n", + "## 7. Clean up 1/f Noise with NSClean (Hand-Modified Mask) \n", "
\n", "\n", "In certain scenarios, manual generation of a mask may be required. Here, we present **one** approach to manually modify the mask for both datasets (shrinking some of the masked FS regions to include more background), starting with the default mask output from the pipeline. It is worth noting that the mask modified using this method may not necessarily outperform the two previous options.\n", @@ -1359,7 +801,7 @@ "id": "775ae18f-eec4-4789-a6c3-2316825ae67a", "metadata": {}, "source": [ - "### 8.1 Verify the Mask (Hand-Modified Mask) \n", + "### 7.1 Verify the Mask (Hand-Modified Mask) \n", "
\n", "\n", "Check the mask against the rate data to make sure it keeps only dark areas of the detector.\n" @@ -1455,7 +897,7 @@ "id": "1c71d4cc-eccc-493a-8bf6-66d06562b94e", "metadata": {}, "source": [ - "### 8.2 Comparing Original vs. Cleaned Data (Hand-Modified Mask) \n", + "### 7.2 Comparing Original vs. Cleaned Data (Hand-Modified Mask) \n", "
" ] }, @@ -1583,7 +1025,7 @@ "id": "3ff38306-8096-4334-b9b8-b55387d406ca", "metadata": {}, "source": [ - "## 9. Conclusion \n", + "## 8. Conclusion \n", "
\n", "\n", "The final plots below show the countrate images and the resulting 1D extracted spectra side-by-side to compare the different cleaning methods: the original (no NSClean applied), the cleaned countrate image (with the default pipeline mask), the cleaned countrate image (with an alternate pipeline mask), and finally, the cleaned countrate image (with the hand-modified mask).\n", @@ -1666,6 +1108,216 @@ "plt.show()" ] }, + { + "cell_type": "code", + "execution_count": null, + "id": "f3e9f088-8dcb-4a80-88ec-d74b51275a63", + "metadata": {}, + "outputs": [], + "source": [ + "def plot_spectra(\n", + " spec_list,\n", + " ext_num=1,\n", + " scale_percent=2.0,\n", + " wavelength_range=None,\n", + " flux_range=None,\n", + " xlim_low=None,\n", + "):\n", + " \"\"\"\n", + " Plot the cleaned extracted 1D spectrum and compare against\n", + " the original (not applying NSClean) extracted 1D spectrum.\n", + "\n", + " Parameters:\n", + " ----------\n", + " spec_list: list\n", + " List of paths to the FITS files containing original/cleaned 1D extracted data.\n", + " scale_percent: float\n", + " Scaling factor applied to determine the intensity range for visualization.\n", + " ext_num: int\n", + " Index/extension of the slice to be plotted.\n", + " The EXTVER header value. The default is 1 for 2D data.\n", + " wavelength_range: dict\n", + " Wavelength range (x-axis) {'nrs1': [3.6, 3.65], 'nrs2': [1.65, 1.75]}.\n", + " flux_range: dict\n", + " Flux range (y-axis) {'nrs1': [1, 2], 'nrs2': [1, 2]}.\n", + " xlim_low: int\n", + " Define a lower wavelength end for the 1D spectrum. Helpful for BOTS data.\n", + " \"\"\"\n", + "\n", + " if wavelength_range is None:\n", + " wavelength_range = {}\n", + "\n", + " # Open the FITS files.\n", + " original_hdul = fits.open(spec_list[0])\n", + " cleaned_hdul = fits.open(spec_list[1])\n", + " if len(spec_list) == 3:\n", + " alternate_hdul = fits.open(spec_list[2])\n", + " elif len(spec_list) == 4:\n", + " alternate_hdul = fits.open(spec_list[2])\n", + " handmod_hdul = fits.open(spec_list[3])\n", + " else:\n", + " alternate_hdul = None\n", + " handmod_hdul = None\n", + "\n", + " # Find the spectral extension (EXTRACT1D).\n", + " for extnum in range(len(original_hdul)):\n", + " hdu = original_hdul[extnum]\n", + " if hdu.name != \"EXTRACT1D\":\n", + " continue\n", + " slit_name = hdu.header[\"SLTNAME\"]\n", + "\n", + " if hdu.header[\"EXTVER\"] == ext_num:\n", + "\n", + " # Plot the original and cleaned spectra together.\n", + " fig, ax = plt.subplots(1, 2, figsize=(10, 5))\n", + " plt.suptitle(\n", + " f\"Compare 1D Spectra for {os.path.basename(spec_list[0])};\"\n", + " f\"EXP_TYPE/Slit = {slit_name}\",\n", + " fontsize=15,\n", + " )\n", + "\n", + " if \"nrs1\" in spec_list[0] and xlim_low is not None:\n", + " xlim_low = xlim_low\n", + " else:\n", + " xlim_low = 0\n", + "\n", + " ax[0].step(\n", + " hdu.data[\"WAVELENGTH\"][xlim_low:],\n", + " hdu.data[\"FLUX\"][xlim_low:],\n", + " linewidth=1,\n", + " label=\"Original\",\n", + " )\n", + " ax[0].step(\n", + " cleaned_hdul[extnum].data[\"WAVELENGTH\"][xlim_low:],\n", + " cleaned_hdul[extnum].data[\"FLUX\"][xlim_low:],\n", + " linewidth=1,\n", + " label=\"Cleaned\",\n", + " )\n", + " ax[0].set_xlabel(f\"Wavelength ({hdu.header['TUNIT1']})\", fontsize=12)\n", + " ax[0].set_ylabel(f\"Flux ({hdu.header['TUNIT2']})\", fontsize=12)\n", + " ax[0].set_title(\"1D Extracted Spectra\", fontsize=12)\n", + "\n", + " # Plot the difference between the spectra as a ratio.\n", + " diff = (\n", + " cleaned_hdul[extnum].data[\"FLUX\"][xlim_low:]\n", + " / hdu.data[\"FLUX\"][xlim_low:]\n", + " )\n", + " ax[1].step(\n", + " hdu.data[\"WAVELENGTH\"][xlim_low:],\n", + " diff,\n", + " linewidth=1,\n", + " label=\"Cleaned/Original\",\n", + " )\n", + " ax[1].set_xlabel(f\"Wavelength ({hdu.header['TUNIT1']})\", fontsize=12)\n", + " ax[1].set_ylabel(\"Cleaned/Original\", fontsize=12)\n", + " ax[1].set_title(\"Difference After Cleaning\", fontsize=12)\n", + " # print(\"Average Difference = {}\".format(np.nanmean(diff)))\n", + "\n", + " # If available, also plot the alternate spectra.\n", + " if alternate_hdul is not None:\n", + " ax[0].step(\n", + " alternate_hdul[extnum].data[\"WAVELENGTH\"][xlim_low:],\n", + " alternate_hdul[extnum].data[\"FLUX\"][xlim_low:],\n", + " linewidth=1,\n", + " label=\"Cleaned (alternate mask)\",\n", + " )\n", + " diff2 = (\n", + " alternate_hdul[extnum].data[\"FLUX\"][xlim_low:]\n", + " / hdu.data[\"FLUX\"][xlim_low:]\n", + " )\n", + " ax[1].step(\n", + " hdu.data[\"WAVELENGTH\"][xlim_low:],\n", + " diff2,\n", + " linewidth=1,\n", + " label=\"Cleaned (alternate mask)/Original\",\n", + " )\n", + " if handmod_hdul is not None:\n", + "\n", + " ax[0].step(\n", + " handmod_hdul[extnum].data[\"WAVELENGTH\"][xlim_low:],\n", + " handmod_hdul[extnum].data[\"FLUX\"][xlim_low:],\n", + " linewidth=1,\n", + " label=\"Cleaned (hand-modified mask)\",\n", + " color=\"Purple\"\n", + " )\n", + " diff3 = (\n", + " handmod_hdul[extnum].data[\"FLUX\"][xlim_low:]\n", + " / hdu.data[\"FLUX\"][xlim_low:]\n", + " )\n", + " ax[1].step(\n", + " hdu.data[\"WAVELENGTH\"][xlim_low:],\n", + " diff3,\n", + " linewidth=1,\n", + " label=\"Cleaned (hand-modified mask)/Original\",\n", + " )\n", + "\n", + " # Set the y-range of the plot if needed.\n", + " if flux_range is not None:\n", + " for key, y_range in flux_range.items():\n", + " if key.lower() in spec_list[0].lower() and y_range is not None:\n", + " if len(y_range) == 2:\n", + " ax[0].set_ylim(y_range)\n", + " ax[1].set_ylim(\n", + " [\n", + " np.nanpercentile(diff, scale_percent),\n", + " np.nanpercentile(diff, 100 - scale_percent),\n", + " ]\n", + " )\n", + " # ax[1].set_ylim(0.5, 1.5)\n", + "\n", + " else:\n", + " all_flux = [\n", + " hdu.data[\"FLUX\"],\n", + " cleaned_hdul[extnum].data[\"FLUX\"],\n", + " ]\n", + " y_range_ax0 = [\n", + " np.nanpercentile(all_flux[0], scale_percent),\n", + " np.nanpercentile(all_flux[0], 100 - scale_percent),\n", + " ]\n", + " ax[0].set_ylim(y_range_ax0)\n", + " # ax[1].set_ylim(0.5, 1.5)\n", + "\n", + " ax[1].set_ylim(\n", + " [\n", + " np.nanpercentile(diff, scale_percent),\n", + " np.nanpercentile(diff, 100 - scale_percent),\n", + " ]\n", + " )\n", + "\n", + " else:\n", + " all_flux = [hdu.data[\"FLUX\"], cleaned_hdul[extnum].data[\"FLUX\"]]\n", + " y_range_ax0 = [\n", + " np.nanpercentile(all_flux[0], scale_percent),\n", + " np.nanpercentile(all_flux[0], 100 - scale_percent),\n", + " ]\n", + " ax[0].set_ylim(y_range_ax0)\n", + " # ax[1].set_ylim(0.5, 1.5)\n", + "\n", + " ax[1].set_ylim(\n", + " [\n", + " np.nanpercentile(diff, scale_percent),\n", + " np.nanpercentile(diff, 100 - scale_percent),\n", + " ]\n", + " )\n", + "\n", + " # Set the x-range of the plot if needed.\n", + " for key in wavelength_range:\n", + " if key in spec_list[0] and wavelength_range[key] is not None:\n", + " ax[0].set_xlim(wavelength_range[key])\n", + " ax[1].set_xlim(wavelength_range[key])\n", + "\n", + " fig.tight_layout()\n", + "\n", + " ax[0].legend()\n", + " ax[1].legend()\n", + "\n", + " original_hdul.close()\n", + " cleaned_hdul.close()\n", + " if alternate_hdul is not None:\n", + " alternate_hdul.close()\n", + " handmod_hdul.close()" + ] + }, { "cell_type": "code", "execution_count": null, @@ -1739,6 +1391,14 @@ "\n", "" ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "4518f00a-c911-4f26-8f38-07a03a076886", + "metadata": {}, + "outputs": [], + "source": [] } ], "metadata": { diff --git a/notebooks/NIRSpec_NSClean/IFU_NSClean_example.ipynb b/notebooks/NIRSpec_NSClean/IFU_NSClean_example.ipynb index 678f73bd3..92c303850 100644 --- a/notebooks/NIRSpec_NSClean/IFU_NSClean_example.ipynb +++ b/notebooks/NIRSpec_NSClean/IFU_NSClean_example.ipynb @@ -18,19 +18,18 @@ "\n", "* 1. [Introduction](#introduction)\n", "* 2. [Import Library](#imports)\n", - "* 3. [Convenience Functions](#functions)\n", - "* 4. [Download the Data](#data)\n", - "* 5. [Running `Spec2Pipeline` without NSClean (Original Data)](#nsclean_skipped)\n", - "* 6. [Clean up 1/f Noise with NSClean (Default Pipeline Mask)](#nsclean_default)\n", - " * 6.1 [Verify the Mask (Default Pipeline Mask)](#verify_default_mask)\n", - " * 6.2 [Comparing Original vs. Cleaned Data (Default Pipeline Mask)](#nsclean_default_compare)\n", - "* 7. [Clean up 1/f Noise with NSClean (Alternate Mask)](#nsclean_alternate)\n", - " * 7.1 [Verify the Mask (Alternate Mask)](#verify_alternate_mask)\n", - " * 7.2 [Comparing Original vs. Cleaned Data (Alternate Mask)](#nsclean_alternate_compare)\n", - "* 8. [Clean up 1/f Noise with NSClean (Hand-Modified Mask)](#nsclean_modified)\n", - " * 8.1 [Verify the Mask (Hand-Modified Mask)](#verify_modified_mask)\n", - " * 8.2 [Comparing Original vs. Cleaned Data (Hand-Modified Mask)](#nsclean_modified_compare)\n", - "* 9. [Conclusion](#conclusion)\n", + "* 3. [Download the Data](#data)\n", + "* 4. [Running `Spec2Pipeline` without NSClean (Original Data)](#nsclean_skipped)\n", + "* 5. [Clean up 1/f Noise with NSClean (Default Pipeline Mask)](#nsclean_default)\n", + " * 5.1 [Verify the Mask (Default Pipeline Mask)](#verify_default_mask)\n", + " * 5.2 [Comparing Original vs. Cleaned Data (Default Pipeline Mask)](#nsclean_default_compare)\n", + "* 6. [Clean up 1/f Noise with NSClean (Alternate Mask)](#nsclean_alternate)\n", + " * 6.1 [Verify the Mask (Alternate Mask)](#verify_alternate_mask)\n", + " * 6.2 [Comparing Original vs. Cleaned Data (Alternate Mask)](#nsclean_alternate_compare)\n", + "* 7. [Clean up 1/f Noise with NSClean (Hand-Modified Mask)](#nsclean_modified)\n", + " * 7.1 [Verify the Mask (Hand-Modified Mask)](#verify_modified_mask)\n", + " * 7.2 [Comparing Original vs. Cleaned Data (Hand-Modified Mask)](#nsclean_modified_compare)\n", + "* 8. [Conclusion](#conclusion)\n", "\n", "## 1. Introduction \n", "
\n", @@ -84,584 +83,27 @@ "import time as tt\n", "import logging\n", "import warnings\n", - "import requests\n", "\n", "# ------ JWST Calibration Pipeline Imports ------\n", "from jwst.pipeline.calwebb_spec2 import Spec2Pipeline\n", "\n", "# ------ Plotting/Stats Imports ------\n", "from matplotlib import pyplot as plt\n", - "from astropy.stats import sigma_clip\n", "from astropy.io import fits\n", "\n", + "from utils import get_jwst_file, plot_dark_data, plot_cleaned_data, plot_spectra\n", + "\n", "# Hide all log and warning messages.\n", "logging.disable(logging.ERROR)\n", "warnings.simplefilter(\"ignore\", RuntimeWarning)" ] }, - { - "cell_type": "markdown", - "id": "b52a3be7-5a74-4420-8221-c4b3f8b2324c", - "metadata": {}, - "source": [ - "## 3. Convenience Functions \n", - "
" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "6d4f3ba0-f15f-4c24-bdbc-aaf77837f21e", - "metadata": {}, - "outputs": [], - "source": [ - "def get_jwst_file(name, mast_api_token=None, save_directory=\".\", redownload=False):\n", - " \"\"\"\n", - " Retrieve a JWST data file from MAST archive and save it to a specified directory.\n", - "\n", - " Parameters:\n", - " ----------\n", - " name: str\n", - " File name.\n", - " mast_api_token: str\n", - " MAST authorization token. Get your MAST Token Here: https://auth.mast.stsci.edu/token.\n", - " save_directory: str\n", - " Save directory path.\n", - " redownload: bool\n", - " Redownload the data even if it exsits already?\n", - " \"\"\"\n", - " mast_url = \"https://mast.stsci.edu/api/v0.1/Download/file\"\n", - " params = dict(uri=f\"mast:JWST/product/{name}\")\n", - "\n", - " if mast_api_token:\n", - " headers = dict(Authorization=f\"token {mast_api_token}\")\n", - " else:\n", - " headers = {}\n", - "\n", - " file_path = os.path.join(save_directory, name)\n", - "\n", - " # Check if the file already exists in the save directory.\n", - " if os.path.exists(file_path) and not redownload:\n", - " print(f\"The file {name} already exists in the directory. Skipping download.\")\n", - " return\n", - "\n", - " r = requests.get(mast_url, params=params, headers=headers, stream=True)\n", - " r.raise_for_status()\n", - "\n", - " with open(file_path, \"wb\") as fobj:\n", - " for chunk in r.iter_content(chunk_size=1024000):\n", - " fobj.write(chunk)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "3015c6d5-07a3-4206-9495-3e9c9aae71a0", - "metadata": {}, - "outputs": [], - "source": [ - "# Define a helper function for plotting the dark areas (the mask).\n", - "\n", - "def plot_dark_data(\n", - " rate_file,\n", - " mask_file,\n", - " cmap=\"viridis\",\n", - " scale=0.2,\n", - " slice_index=1,\n", - " axis=0,\n", - " aspect=\"auto\",\n", - " layout=\"rows\",\n", - "):\n", - " \"\"\"\n", - " Plot the dark areas on the detector, masking out the illuminated areas.\n", - " This function can take 2D (_rate.fits) and 3D (_rateints.fits) files as input.\n", - "\n", - " Parameters:\n", - " ----------\n", - " rate_file: str\n", - " File path to the FITS file containing countrate data.\n", - " mask_file: str\n", - " File path to the FITS file containing mask data.\n", - " slice_index: int\n", - " Index of the slice to be plotted, 2D data default is 1.\n", - " axis: int (0-2)\n", - " Axis along which the slice will be taken\n", - " (0 for x-axis, 1 for y-axis, 2 for z-axis).\n", - " aspect: int or 'auto'\n", - " Plot's aspect ratio.\n", - " cmap: str\n", - " Color map.\n", - " layout: str\n", - " Layout subplots in rows or columns?\n", - " scale: float\n", - " Scaling factor applied to determine the intensity range for visualization.\n", - " \"\"\"\n", - "\n", - " # Make a matplotlib figure.\n", - " if layout == \"rows\":\n", - " fig, ax = plt.subplots(3, 1, figsize=(10, 5))\n", - " else:\n", - " fig, ax = plt.subplots(1, 3, figsize=(20, 5))\n", - " plt.suptitle(f\"Dark Areas for {os.path.basename(rate_file)}\", fontsize=15)\n", - "\n", - " # Open the mask and science data.\n", - " with fits.open(mask_file) as hdulist:\n", - " mask = hdulist[1].data\n", - "\n", - " with fits.open(rate_file) as hdulist:\n", - " sci = hdulist[\"SCI\"].data\n", - "\n", - " # Determine if countrate data is 2D or 3D.\n", - " if len(sci.shape) == 3:\n", - " dim = 3\n", - " else:\n", - " dim = 2\n", - "\n", - " # Get data limits from the dark data.\n", - " masked_sci = sci.copy()\n", - " masked_sci[mask == 0] = 0\n", - " vmin = np.nanpercentile(masked_sci, scale)\n", - " vmax = np.nanpercentile(masked_sci, 100 - scale)\n", - "\n", - " # Plot the science image with limits from the dark data.\n", - " sci[np.isnan(sci)] = 0\n", - "\n", - " # If the countrate data is 3D, determine the integration (slice) to plot.\n", - " if dim == 3:\n", - " sci = np.take(sci, slice_index, axis=axis)\n", - " mask = np.take(mask, slice_index, axis=axis)\n", - " masked_sci = np.take(masked_sci, slice_index, axis=axis)\n", - "\n", - " ax[0].set_title(\"Original Rate Data: Integration [{},:,:]\".format(slice_index))\n", - " ax[1].set_title(\n", - " \"Illuminated Pixel Mask: Integration [{},:,:]\".format(slice_index)\n", - " )\n", - " ax[2].set_title(\n", - " \"Dark Data (Illuminated Pixels Masked): Integration [{},:,:]\".format(\n", - " slice_index\n", - " )\n", - " )\n", - "\n", - " else:\n", - " sci = sci\n", - " ax[0].set_title(\"Original Rate Data\")\n", - " ax[1].set_title(\"Illuminated Pixel Mask\")\n", - " ax[2].set_title(\"Dark Data (Illuminated Pixels Masked)\")\n", - "\n", - " # Plot the science image with limits from the dark data.\n", - " ax[0].set_ylabel(\"Pixel Row\", fontsize=12)\n", - " ax[0].set_xlabel(\"Pixel Column\", fontsize=12)\n", - " im1 = ax[0].imshow(\n", - " sci, cmap=cmap, origin=\"lower\", aspect=aspect, vmin=vmin, vmax=vmax\n", - " )\n", - " cbar1 = fig.colorbar(im1, ax=ax[0], pad=0.05, shrink=0.7, label=\"DN/s\")\n", - "\n", - " # Plot the mask: values are 1 or 0.\n", - " ax[1].set_ylabel(\"Pixel Row\", fontsize=12)\n", - " ax[1].set_xlabel(\"Pixel Column\", fontsize=12)\n", - " im2 = ax[1].imshow(mask, cmap=cmap, aspect=aspect, origin=\"lower\", vmin=0, vmax=1)\n", - " cbar2 = fig.colorbar(im2, ax=ax[1], pad=0.05, shrink=0.7, ticks=[0, 1])\n", - " cbar2.set_ticklabels([\"Illuminated Pixel Masked\", \"Un-Illuminated Pixel\"])\n", - "\n", - " # Plot the dark data with the same limits as the science data.\n", - " ax[2].set_ylabel(\"Pixel Row\", fontsize=12)\n", - " ax[2].set_xlabel(\"Pixel Column\", fontsize=12)\n", - " im3 = ax[2].imshow(\n", - " masked_sci, cmap=cmap, aspect=aspect, origin=\"lower\", vmin=vmin, vmax=vmax\n", - " )\n", - " cbar3 = fig.colorbar(im3, ax=ax[2], pad=0.05, shrink=0.7, label=\"DN/s\")\n", - "\n", - " if layout == \"rows\": # Adjust the multiplier as needed.\n", - " cbar1.ax.figure.set_size_inches(cbar1.ax.figure.get_size_inches() * 1.2)\n", - " cbar2.ax.figure.set_size_inches(cbar2.ax.figure.get_size_inches() * 1.2)\n", - " cbar3.ax.figure.set_size_inches(cbar3.ax.figure.get_size_inches() * 1.2)\n", - "\n", - " fig.tight_layout() # Adjusted tight_layout for better suptitle spacing." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "5b195d53-9fed-4bba-b2f4-36e07a2e2392", - "metadata": {}, - "outputs": [], - "source": [ - "# Define a helper function for plotting the cleaned data.\n", - "\n", - "def plot_cleaned_data(\n", - " rate_file,\n", - " cleaned_file,\n", - " slice_indx=1,\n", - " aspect=\"auto\",\n", - " cmap=\"viridis\",\n", - " scale=0.2,\n", - " layout=\"rows\",\n", - "):\n", - " \"\"\"\n", - " Plot the 2D cleaned rate data (or 2D slice from a 3D cleaned data cube)\n", - " and compare against the original data (not applying NSClean).\n", - "\n", - " Parameters:\n", - " ----------\n", - " rate_file: str\n", - " File path to the FITS file containing original countrate data.\n", - " cleaned_file: str\n", - " File path to the FITS file containing cleaned countrate data.\n", - " slice_indx: int\n", - " Index of the slice to be plotted, 2D data default is 1.\n", - " aspect: int or 'auto'\n", - " Plot's aspect ratio.\n", - " cmap: str\n", - " Color map.\n", - " layout: str\n", - " Layout subplots in rows or columns?\n", - " scale: float\n", - " Scaling factor applied to determine the intensity range for visualization.\n", - " \"\"\"\n", - "\n", - " # Make a matplotlib figure.\n", - " if layout == \"rows\":\n", - " fig, ax = plt.subplots(4, 1, figsize=(10, 12))\n", - " else:\n", - " fig, ax = plt.subplots(2, 2, figsize=(12, 10))\n", - " ax = ax.flatten()\n", - " plt.suptitle(\n", - " f\"Cleaned Data (1/f noise removed) for {os.path.basename(rate_file)}\",\n", - " fontsize=15,\n", - " y=1,\n", - " )\n", - "\n", - " # Open the original and cleaned data.\n", - " with fits.open(rate_file) as hdulist:\n", - " original = hdulist[\"SCI\"].data\n", - "\n", - " with fits.open(cleaned_file) as hdulist:\n", - " cleaned = hdulist[\"SCI\"].data\n", - " # Determine if rate data is 2D or 3D.\n", - " if len(cleaned.shape) == 3:\n", - " dim = 3\n", - " else:\n", - " dim = 2\n", - "\n", - " # Define image limits from the original data.\n", - " vmin = np.nanpercentile(original, scale)\n", - " vmax = np.nanpercentile(original, 100 - scale)\n", - "\n", - " original[np.isnan(original)] = 0\n", - " cleaned[np.isnan(cleaned)] = 0\n", - "\n", - " # If the rate data is 3D, determine the integration (slice) to plot.\n", - " if dim == 3:\n", - " original = original[slice_indx, :, :]\n", - " cleaned = cleaned[slice_indx, :, :]\n", - " ax[0].set_title(\n", - " \"Original Rate Data: Integration [{},:,:]\".format(slice_indx), fontsize=15\n", - " )\n", - " ax[1].set_title(\n", - " \"Cleaned Rate Data: Integration [{},:,:]\".format(slice_indx), fontsize=15\n", - " )\n", - "\n", - " else:\n", - " ax[0].set_title(\"Original Rate Data\", fontsize=12)\n", - " ax[1].set_title(\"Cleaned Rate Data\", fontsize=12)\n", - "\n", - " # Plot the original rate data.\n", - " ax[0].set_xlabel(\"Pixel Column\", fontsize=10)\n", - " ax[0].set_ylabel(\"Pixel Row\", fontsize=10)\n", - " fig.colorbar(\n", - " ax[0].imshow(\n", - " original, cmap=cmap, origin=\"lower\", aspect=aspect, vmin=vmin, vmax=vmax\n", - " ),\n", - " ax=ax[0],\n", - " pad=0.05,\n", - " shrink=0.7,\n", - " label=\"DN/s\",\n", - " )\n", - "\n", - " # Plot the cleaned data with the same image limits.\n", - " ax[1].set_xlabel(\"Pixel Column\", fontsize=10)\n", - " ax[1].set_ylabel(\"Pixel Row\", fontsize=10)\n", - " fig.colorbar(\n", - " ax[1].imshow(\n", - " cleaned, cmap=cmap, origin=\"lower\", aspect=aspect, vmin=vmin, vmax=vmax\n", - " ),\n", - " ax=ax[1],\n", - " pad=0.05,\n", - " shrink=0.7,\n", - " label=\"DN/s\",\n", - " )\n", - "\n", - " # Plot the relative difference between the original and cleaned data.\n", - " diff = (original - cleaned) / original\n", - " diff[~np.isfinite(diff)] = 0\n", - " vmin_diff = np.nanpercentile(diff, scale)\n", - " vmax_diff = np.nanpercentile(diff, 100 - scale)\n", - " ax[2].set_title(\"Relative Difference (Original - Cleaned Rate Data)\", fontsize=12)\n", - " ax[2].set_xlabel(\"Pixel Column\", fontsize=10)\n", - " ax[2].set_ylabel(\"Pixel Row\", fontsize=10)\n", - " fig.colorbar(\n", - " ax[2].imshow(\n", - " diff,\n", - " cmap=cmap,\n", - " origin=\"lower\",\n", - " aspect=aspect,\n", - " vmin=vmin_diff,\n", - " vmax=vmax_diff,\n", - " ),\n", - " ax=ax[2],\n", - " pad=0.05,\n", - " shrink=0.7,\n", - " label=\"DN/s\",\n", - " )\n", - "\n", - " # Plot the relative difference again,\n", - " # this time hiding the outliers so that low-level\n", - " # background changes can be seen.\n", - " hide_outliers = np.ma.filled(sigma_clip(diff, masked=True), fill_value=0)\n", - " vmin_outliers = np.nanpercentile(hide_outliers, scale)\n", - " vmax_outliers = np.nanpercentile(hide_outliers, 100 - scale)\n", - " ax[3].set_title(\n", - " \"Relative Difference (Original - Cleaned Rate Data) \\n with Outliers Hidden\",\n", - " fontsize=12,\n", - " )\n", - " ax[3].set_xlabel(\"Pixel Column\", fontsize=10)\n", - " ax[3].set_ylabel(\"Pixel Row\", fontsize=10)\n", - " fig.colorbar(\n", - " ax[3].imshow(\n", - " hide_outliers,\n", - " cmap=cmap,\n", - " origin=\"lower\",\n", - " aspect=aspect,\n", - " vmin=vmin_outliers,\n", - " vmax=vmax_outliers,\n", - " ),\n", - " ax=ax[3],\n", - " pad=0.05,\n", - " shrink=0.7,\n", - " label=\"DN/s\",\n", - " )\n", - "\n", - " fig.tight_layout()" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "ad6fae03-79fc-479c-9802-d3484b23283a", - "metadata": {}, - "outputs": [], - "source": [ - "# Define a helper function for plotting the 1D spectra.\n", - "\n", - "def plot_spectra(\n", - " spec_list,\n", - " ext_num=1,\n", - " scale_percent=2.0,\n", - " wavelength_range=None,\n", - " flux_range=None,\n", - " xlim_low=None,\n", - "):\n", - " \"\"\"\n", - " Plot the cleaned extracted 1D spectrum and compare against\n", - " the original (not applying NSClean) extracted 1D spectrum.\n", - "\n", - " Parameters:\n", - " ----------\n", - " spec_list: list\n", - " List of paths to the FITS files containing original/cleaned 1D extracted data.\n", - " scale_percent: float\n", - " Scaling factor applied to determine the intensity range for visualization.\n", - " ext_num: int\n", - " Index/extension of the slice to be plotted.\n", - " The EXTVER header value. The default is 1 for 2D data.\n", - " wavelength_range: dict\n", - " Wavelength range (x-axis) {'nrs1': [3.6, 3.65], 'nrs2': [1.65, 1.75]}.\n", - " flux_range: dict\n", - " Flux range (y-axis) {'nrs1': [1, 2], 'nrs2': [1, 2]}.\n", - " xlim_low: int\n", - " Define a lower wavelength end for the 1D spectrum. Helpful for BOTS data.\n", - " \"\"\"\n", - "\n", - " if wavelength_range is None:\n", - " wavelength_range = {}\n", - "\n", - " # Open the FITS files.\n", - " original_hdul = fits.open(spec_list[0])\n", - " cleaned_hdul = fits.open(spec_list[1])\n", - " if len(spec_list) == 3:\n", - " alternate_hdul = fits.open(spec_list[2])\n", - " elif len(spec_list) == 4:\n", - " alternate_hdul = fits.open(spec_list[2])\n", - " handmod_hdul = fits.open(spec_list[3])\n", - " else:\n", - " alternate_hdul = None\n", - " handmod_hdul = None\n", - "\n", - " # Find the spectral extension (EXTRACT1D).\n", - " for extnum in range(len(original_hdul)):\n", - " hdu = original_hdul[extnum]\n", - " if hdu.name != \"EXTRACT1D\":\n", - " continue\n", - " slit_name = hdu.header[\"SLTNAME\"]\n", - "\n", - " if hdu.header[\"EXTVER\"] == ext_num:\n", - "\n", - " # Plot the original and cleaned spectra together.\n", - " fig, ax = plt.subplots(1, 2, figsize=(10, 5))\n", - " plt.suptitle(\n", - " f\"Compare 1D Spectra for {os.path.basename(spec_list[0])};\"\n", - " f\"EXP_TYPE/Slit = {slit_name}\",\n", - " fontsize=15,\n", - " )\n", - "\n", - " if \"nrs1\" in spec_list[0] and xlim_low is not None:\n", - " xlim_low = xlim_low\n", - " else:\n", - " xlim_low = 0\n", - "\n", - " ax[0].step(\n", - " hdu.data[\"WAVELENGTH\"][xlim_low:],\n", - " hdu.data[\"FLUX\"][xlim_low:],\n", - " linewidth=1,\n", - " label=\"Original\",\n", - " )\n", - " ax[0].step(\n", - " cleaned_hdul[extnum].data[\"WAVELENGTH\"][xlim_low:],\n", - " cleaned_hdul[extnum].data[\"FLUX\"][xlim_low:],\n", - " linewidth=1,\n", - " label=\"Cleaned\",\n", - " )\n", - " ax[0].set_xlabel(f\"Wavelength ({hdu.header['TUNIT1']})\", fontsize=12)\n", - " ax[0].set_ylabel(f\"Flux ({hdu.header['TUNIT2']})\", fontsize=12)\n", - " ax[0].set_title(\"1D Extracted Spectra\", fontsize=12)\n", - "\n", - " # Plot the difference between the spectra as a ratio.\n", - " diff = (\n", - " cleaned_hdul[extnum].data[\"FLUX\"][xlim_low:]\n", - " / hdu.data[\"FLUX\"][xlim_low:]\n", - " )\n", - " ax[1].step(\n", - " hdu.data[\"WAVELENGTH\"][xlim_low:],\n", - " diff,\n", - " linewidth=1,\n", - " label=\"Cleaned/Original\",\n", - " )\n", - " ax[1].set_xlabel(f\"Wavelength ({hdu.header['TUNIT1']})\", fontsize=12)\n", - " ax[1].set_ylabel(\"Cleaned/Original\", fontsize=12)\n", - " ax[1].set_title(\"Difference After Cleaning\", fontsize=12)\n", - " # print(\"Average Difference = {}\".format(np.nanmean(diff)))\n", - "\n", - " # If available, also plot the alternate spectra.\n", - " if alternate_hdul is not None:\n", - " ax[0].step(\n", - " alternate_hdul[extnum].data[\"WAVELENGTH\"][xlim_low:],\n", - " alternate_hdul[extnum].data[\"FLUX\"][xlim_low:],\n", - " linewidth=1,\n", - " label=\"Cleaned (alternate mask)\",\n", - " )\n", - " diff2 = (\n", - " alternate_hdul[extnum].data[\"FLUX\"][xlim_low:]\n", - " / hdu.data[\"FLUX\"][xlim_low:]\n", - " )\n", - " ax[1].step(\n", - " hdu.data[\"WAVELENGTH\"][xlim_low:],\n", - " diff2,\n", - " linewidth=1,\n", - " label=\"Cleaned (alternate mask)/Original\",\n", - " )\n", - " if handmod_hdul is not None:\n", - "\n", - " ax[0].step(\n", - " handmod_hdul[extnum].data[\"WAVELENGTH\"][xlim_low:],\n", - " handmod_hdul[extnum].data[\"FLUX\"][xlim_low:],\n", - " linewidth=1,\n", - " label=\"Cleaned (hand-modified mask)\",\n", - " color=\"Purple\"\n", - " )\n", - " diff3 = (\n", - " handmod_hdul[extnum].data[\"FLUX\"][xlim_low:]\n", - " / hdu.data[\"FLUX\"][xlim_low:]\n", - " )\n", - " ax[1].step(\n", - " hdu.data[\"WAVELENGTH\"][xlim_low:],\n", - " diff3,\n", - " linewidth=1,\n", - " label=\"Cleaned (hand-modified mask)/Original\",\n", - " )\n", - "\n", - " # Set the y-range of the plot if needed.\n", - " if flux_range is not None:\n", - " for key, y_range in flux_range.items():\n", - " if key.lower() in spec_list[0].lower() and y_range is not None:\n", - " if len(y_range) == 2:\n", - " ax[0].set_ylim(y_range)\n", - " ax[1].set_ylim(\n", - " [\n", - " np.nanpercentile(diff, scale_percent),\n", - " np.nanpercentile(diff, 100 - scale_percent),\n", - " ]\n", - " )\n", - " # ax[1].set_ylim(0.5, 1.5)\n", - "\n", - " else:\n", - " all_flux = [\n", - " hdu.data[\"FLUX\"],\n", - " cleaned_hdul[extnum].data[\"FLUX\"],\n", - " ]\n", - " y_range_ax0 = [\n", - " np.nanpercentile(all_flux[0], scale_percent),\n", - " np.nanpercentile(all_flux[0], 100 - scale_percent),\n", - " ]\n", - " ax[0].set_ylim(y_range_ax0)\n", - " # ax[1].set_ylim(0.5, 1.5)\n", - "\n", - " ax[1].set_ylim(\n", - " [\n", - " np.nanpercentile(diff, scale_percent),\n", - " np.nanpercentile(diff, 100 - scale_percent),\n", - " ]\n", - " )\n", - "\n", - " else:\n", - " all_flux = [hdu.data[\"FLUX\"], cleaned_hdul[extnum].data[\"FLUX\"]]\n", - " y_range_ax0 = [\n", - " np.nanpercentile(all_flux[0], scale_percent),\n", - " np.nanpercentile(all_flux[0], 100 - scale_percent),\n", - " ]\n", - " ax[0].set_ylim(y_range_ax0)\n", - " # ax[1].set_ylim(0.5, 1.5)\n", - "\n", - " ax[1].set_ylim(\n", - " [\n", - " np.nanpercentile(diff, scale_percent),\n", - " np.nanpercentile(diff, 100 - scale_percent),\n", - " ]\n", - " )\n", - "\n", - " # Set the x-range of the plot if needed.\n", - " for key in wavelength_range:\n", - " if key in spec_list[0] and wavelength_range[key] is not None:\n", - " ax[0].set_xlim(wavelength_range[key])\n", - " ax[1].set_xlim(wavelength_range[key])\n", - "\n", - " fig.tight_layout()\n", - "\n", - " ax[0].legend()\n", - " ax[1].legend()\n", - "\n", - " original_hdul.close()\n", - " cleaned_hdul.close()\n", - " if alternate_hdul is not None:\n", - " alternate_hdul.close()\n", - " handmod_hdul.close()" - ] - }, { "cell_type": "markdown", "id": "58bafbdb-d55b-4e31-bc4e-97e0fb8f67b1", "metadata": {}, "source": [ - "## 4. Download the Data \n", + "## 3. Download the Data \n", "
\n", " \n", "The input data for this notebook features an IFU observation of quasar XID2028 with grating/filter G140H/F100LP. The dataset is part of the [JWST Early Release Science program ERS-1335](https://www.stsci.edu/jwst/science-execution/approved-programs/dd-ers/program-1335), specifically observation 4. It consists of 9 integrations (9 dither points) with 16 groups each. This notebook focuses on the second dithered exposure (00002) as an example. However, it's important to note that before proceeding to the `Spec3Pipeline`, all exposures must first be processed through the `Spec2Pipeline`." @@ -711,7 +153,7 @@ "id": "302c90c6-3c42-4710-becc-49533b69de70", "metadata": {}, "source": [ - "## 5. Running `Spec2Pipeline` without NSClean (Original Data) \n", + "## 4. Running `Spec2Pipeline` without NSClean (Original Data) \n", "
\n", "\n", "The cell below executes the `Spec2Pipeline`, explicitly skipping the NSClean step during processing. The level 2 products generated will serve as a reference point to illustrate how the countrate images and final extracted spectra appear without the removal of 1/f noise." @@ -768,7 +210,7 @@ "id": "2b99dc71-32ca-4599-b598-e64ea206f9d2", "metadata": {}, "source": [ - "## 6. Clean up 1/f Noise with NSClean (Default Pipeline Mask) \n", + "## 5. Clean up 1/f Noise with NSClean (Default Pipeline Mask) \n", "
\n", "\n", "If a user-supplied mask file is not provided to the [NSClean step](https://jwst-pipeline.readthedocs.io/en/latest/jwst/nsclean/index.html) in the `Spec2Pipeline`, the pipeline will generate a mask based on default parameters. This mask will identify any pixel that is unilluminated. That is, the mask must contain True and False values, where True indicates that the pixel is dark, and False indicates that the pixel is illuminated (not dark).\n", @@ -852,7 +294,7 @@ "id": "50c5f2cf-d7a3-4606-bc92-46b3b73cfb29", "metadata": {}, "source": [ - "### 6.1 Verify the Mask (Default Pipeline Mask) \n", + "### 5.1 Verify the Mask (Default Pipeline Mask) \n", "
\n", "\n", "Check the mask against the rate data to make sure it keeps only dark areas of the detector.\n", @@ -885,7 +327,7 @@ "id": "9b4c5b7b-4dd6-4d6c-9ae3-5590ab9131b9", "metadata": {}, "source": [ - "### 6.2 Comparing Original vs. Cleaned Data (Default Pipeline Mask) \n", + "### 5.2 Comparing Original vs. Cleaned Data (Default Pipeline Mask) \n", "
\n", "\n", "We can now compare the cleaned data (with the default pipeline mask) to the original rate file and verify that the 1/f noise has been reduced.\n", @@ -974,7 +416,7 @@ "id": "a1f85a17-cb1c-4320-a3bd-1800fa02bab9", "metadata": {}, "source": [ - "## 7. Clean up 1/f Noise with NSClean (Alternate Mask) \n", + "## 6. Clean up 1/f Noise with NSClean (Alternate Mask) \n", "
\n", "\n", "For some data sets, masking the entire science region may excessively mask dark areas of the detector that could be used to improve the background fit. Excessive masking can introduce some high frequency noise in the cleaning process that appears as vertical striping over the spectral traces. Also, for some data sets, there may be several illuminated regions of the detector that are not masked by the IFU slice bounding boxes. This may include transient artifacts like cosmic rays or glow from saturated sources.\n", @@ -1044,7 +486,7 @@ "id": "ef33421e-bd16-41d0-883d-4a68202897b6", "metadata": {}, "source": [ - "### 7.1 Verify the Mask (Alternate Mask) \n", + "### 6.1 Verify the Mask (Alternate Mask) \n", "
\n", "\n", "Check the mask against the rate data to make sure it keeps only dark areas of the detector.\n" @@ -1075,7 +517,7 @@ "id": "b84cca59-2781-437d-a169-fa5c9e012b3a", "metadata": {}, "source": [ - "### 7.2 Comparing Original vs. Cleaned Data (Alternate Mask) \n", + "### 6.2 Comparing Original vs. Cleaned Data (Alternate Mask) \n", "
" ] }, @@ -1150,7 +592,7 @@ "id": "fec75731", "metadata": {}, "source": [ - "## 8. Clean up 1/f Noise with NSClean (Hand-Modified Mask) \n", + "## 7. Clean up 1/f Noise with NSClean (Hand-Modified Mask) \n", "
\n", "\n", "In certain scenarios, manual generation of a mask may be required. Here, we present **one** approach to manually modify the mask (excluding some large snowballs in NRS1/NRS2) starting with the default mask output from the pipeline. It is worth noting that the mask modified using this method may not necessarily outperform the two previous options." @@ -1222,7 +664,7 @@ "id": "463ae320-2f8e-4c5a-9a47-d1698c33704e", "metadata": {}, "source": [ - "### 8.1 Verify the Mask (Hand-Modified Mask) \n", + "### 7.1 Verify the Mask (Hand-Modified Mask) \n", "
\n", "\n", "Check the mask against the rate data to make sure it keeps only dark areas of the detector.\n" @@ -1299,7 +741,7 @@ "id": "36753340-5509-4fec-91a4-3b1dcfe63075", "metadata": {}, "source": [ - "### 8.2 Comparing Original vs. Cleaned Data (Hand-Modified Mask) \n", + "### 7.2 Comparing Original vs. Cleaned Data (Hand-Modified Mask) \n", "
" ] }, @@ -1374,7 +816,7 @@ "id": "f86daa01-29bd-4ce4-9a63-0091e0dbcf2b", "metadata": {}, "source": [ - "## 9. Conclusion \n", + "## 8. Conclusion \n", "
\n", "\n", "The final plots below show the countrate images and the resulting 1D extracted spectra side-by-side to compare the different cleaning methods: the original (no NSClean applied), the cleaned countrate image (with the default pipeline mask), the cleaned countrate image (with an alternate pipeline mask), and finally, the cleaned countrate image (with the hand-modified mask).\n", diff --git a/notebooks/NIRSpec_NSClean/MOS_NSClean_example.ipynb b/notebooks/NIRSpec_NSClean/MOS_NSClean_example.ipynb index c0b4da5e6..fbc4067af 100644 --- a/notebooks/NIRSpec_NSClean/MOS_NSClean_example.ipynb +++ b/notebooks/NIRSpec_NSClean/MOS_NSClean_example.ipynb @@ -18,19 +18,18 @@ "\n", "* 1. [Introduction](#introduction)\n", "* 2. [Import Library](#imports)\n", - "* 3. [Convenience Functions](#functions)\n", - "* 4. [Download the Data](#data)\n", - "* 5. [Running `Spec2Pipeline` without NSClean (Original Data)](#nsclean_skipped)\n", - "* 6. [Clean up 1/f Noise with NSClean (Default Pipeline Mask)](#nsclean_default)\n", - " * 6.1 [Verify the Mask (Default Pipeline Mask)](#verify_default_mask)\n", - " * 6.2 [Comparing Original vs. Cleaned Data (Default Pipeline Mask)](#nsclean_default_compare)\n", - "* 7. [Clean up 1/f Noise with NSClean (Alternate Mask)](#nsclean_alternate)\n", - " * 7.1 [Verify the Mask (Alternate Mask)](#verify_alternate_mask)\n", - " * 7.2 [Comparing Original vs. Cleaned Data (Alternate Mask)](#nsclean_alternate_compare)\n", - "* 8. [Clean up 1/f Noise with NSClean (Hand-Modified Mask)](#nsclean_modified)\n", - " * 8.1 [Verify the Mask (Hand-Modified Mask)](#verify_modified_mask)\n", - " * 8.2 [Comparing Original vs. Cleaned Data (Hand-Modified Mask)](#nsclean_modified_compare)\n", - "* 9. [Conclusion](#conclusion)\n", + "* 3. [Download the Data](#data)\n", + "* 4. [Running `Spec2Pipeline` without NSClean (Original Data)](#nsclean_skipped)\n", + "* 5. [Clean up 1/f Noise with NSClean (Default Pipeline Mask)](#nsclean_default)\n", + " * 5.1 [Verify the Mask (Default Pipeline Mask)](#verify_default_mask)\n", + " * 5.2 [Comparing Original vs. Cleaned Data (Default Pipeline Mask)](#nsclean_default_compare)\n", + "* 6. [Clean up 1/f Noise with NSClean (Alternate Mask)](#nsclean_alternate)\n", + " * 6.1 [Verify the Mask (Alternate Mask)](#verify_alternate_mask)\n", + " * 6.2 [Comparing Original vs. Cleaned Data (Alternate Mask)](#nsclean_alternate_compare)\n", + "* 7. [Clean up 1/f Noise with NSClean (Hand-Modified Mask)](#nsclean_modified)\n", + " * 7.1 [Verify the Mask (Hand-Modified Mask)](#verify_modified_mask)\n", + " * 7.2 [Comparing Original vs. Cleaned Data (Hand-Modified Mask)](#nsclean_modified_compare)\n", + "* 8. [Conclusion](#conclusion)\n", "\n", "\n", "## 1. Introduction \n", @@ -85,584 +84,27 @@ "import time as tt\n", "import logging\n", "import warnings\n", - "import requests\n", "\n", "# ------ JWST Calibration Pipeline Imports ------\n", "from jwst.pipeline.calwebb_spec2 import Spec2Pipeline\n", "\n", "# ------ Plotting/Stats Imports ------\n", "from matplotlib import pyplot as plt\n", - "from astropy.stats import sigma_clip\n", "from astropy.io import fits\n", "\n", + "from utils import get_jwst_file, plot_dark_data, plot_cleaned_data, plot_spectra\n", + "\n", "# Hide all log and warning messages.\n", "logging.disable(logging.ERROR)\n", "warnings.simplefilter(\"ignore\", RuntimeWarning)" ] }, - { - "cell_type": "markdown", - "id": "b1baac4d-ea9b-4d7f-9ea8-5d367cd52ae1", - "metadata": {}, - "source": [ - "## 3. Convenience Functions \n", - "
" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "e3685b99-9dad-45a0-8f76-3eb79621610c", - "metadata": {}, - "outputs": [], - "source": [ - "def get_jwst_file(name, mast_api_token=None, save_directory=\".\", redownload=False):\n", - " \"\"\"\n", - " Retrieve a JWST data file from MAST archive and save it to a specified directory.\n", - "\n", - " Parameters:\n", - " ----------\n", - " name: str\n", - " File name.\n", - " mast_api_token: str\n", - " MAST authorization token. Get your MAST Token Here: https://auth.mast.stsci.edu/token.\n", - " save_directory: str\n", - " Save directory path.\n", - " redownload: bool\n", - " Redownload the data even if it exsits already?\n", - " \"\"\"\n", - " mast_url = \"https://mast.stsci.edu/api/v0.1/Download/file\"\n", - " params = dict(uri=f\"mast:JWST/product/{name}\")\n", - "\n", - " if mast_api_token:\n", - " headers = dict(Authorization=f\"token {mast_api_token}\")\n", - " else:\n", - " headers = {}\n", - "\n", - " file_path = os.path.join(save_directory, name)\n", - "\n", - " # Check if the file already exists in the save directory.\n", - " if os.path.exists(file_path) and not redownload:\n", - " print(f\"The file {name} already exists in the directory. Skipping download.\")\n", - " return\n", - "\n", - " r = requests.get(mast_url, params=params, headers=headers, stream=True)\n", - " r.raise_for_status()\n", - "\n", - " with open(file_path, \"wb\") as fobj:\n", - " for chunk in r.iter_content(chunk_size=1024000):\n", - " fobj.write(chunk)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "dd9a41c1-5fa1-4bd4-bcf3-1f23ec166266", - "metadata": {}, - "outputs": [], - "source": [ - "# Define a helper function for plotting the dark areas (the mask).\n", - "\n", - "def plot_dark_data(\n", - " rate_file,\n", - " mask_file,\n", - " cmap=\"viridis\",\n", - " scale=0.2,\n", - " slice_index=1,\n", - " axis=0,\n", - " aspect=\"auto\",\n", - " layout=\"rows\",\n", - "):\n", - " \"\"\"\n", - " Plot the dark areas on the detector, masking out the illuminated areas.\n", - " This function can take 2D (_rate.fits) and 3D (_rateints.fits) files as input.\n", - "\n", - " Parameters:\n", - " ----------\n", - " rate_file: str\n", - " File path to the FITS file containing countrate data.\n", - " mask_file: str\n", - " File path to the FITS file containing mask data.\n", - " slice_index: int\n", - " Index of the slice to be plotted, 2D data default is 1.\n", - " axis: int (0-2)\n", - " Axis along which the slice will be taken\n", - " (0 for x-axis, 1 for y-axis, 2 for z-axis).\n", - " aspect: int or 'auto'\n", - " Plot's aspect ratio.\n", - " cmap: str\n", - " Color map.\n", - " layout: str\n", - " Layout subplots in rows or columns?\n", - " scale: float\n", - " Scaling factor applied to determine the intensity range for visualization.\n", - " \"\"\"\n", - "\n", - " # Make a matplotlib figure.\n", - " if layout == \"rows\":\n", - " fig, ax = plt.subplots(3, 1, figsize=(10, 5))\n", - " else:\n", - " fig, ax = plt.subplots(1, 3, figsize=(20, 5))\n", - " plt.suptitle(f\"Dark Areas for {os.path.basename(rate_file)}\", fontsize=15)\n", - "\n", - " # Open the mask and science data.\n", - " with fits.open(mask_file) as hdulist:\n", - " mask = hdulist[1].data\n", - "\n", - " with fits.open(rate_file) as hdulist:\n", - " sci = hdulist[\"SCI\"].data\n", - "\n", - " # Determine if countrate data is 2D or 3D.\n", - " if len(sci.shape) == 3:\n", - " dim = 3\n", - " else:\n", - " dim = 2\n", - "\n", - " # Get data limits from the dark data.\n", - " masked_sci = sci.copy()\n", - " masked_sci[mask == 0] = 0\n", - " vmin = np.nanpercentile(masked_sci, scale)\n", - " vmax = np.nanpercentile(masked_sci, 100 - scale)\n", - "\n", - " # Plot the science image with limits from the dark data.\n", - " sci[np.isnan(sci)] = 0\n", - "\n", - " # If the countrate data is 3D, determine the integration (slice) to plot.\n", - " if dim == 3:\n", - " sci = np.take(sci, slice_index, axis=axis)\n", - " mask = np.take(mask, slice_index, axis=axis)\n", - " masked_sci = np.take(masked_sci, slice_index, axis=axis)\n", - "\n", - " ax[0].set_title(\"Original Rate Data: Integration [{},:,:]\".format(slice_index))\n", - " ax[1].set_title(\n", - " \"Illuminated Pixel Mask: Integration [{},:,:]\".format(slice_index)\n", - " )\n", - " ax[2].set_title(\n", - " \"Dark Data (Illuminated Pixels Masked): Integration [{},:,:]\".format(\n", - " slice_index\n", - " )\n", - " )\n", - "\n", - " else:\n", - " sci = sci\n", - " ax[0].set_title(\"Original Rate Data\")\n", - " ax[1].set_title(\"Illuminated Pixel Mask\")\n", - " ax[2].set_title(\"Dark Data (Illuminated Pixels Masked)\")\n", - "\n", - " # Plot the science image with limits from the dark data.\n", - " ax[0].set_ylabel(\"Pixel Row\", fontsize=12)\n", - " ax[0].set_xlabel(\"Pixel Column\", fontsize=12)\n", - " im1 = ax[0].imshow(\n", - " sci, cmap=cmap, origin=\"lower\", aspect=aspect, vmin=vmin, vmax=vmax\n", - " )\n", - " cbar1 = fig.colorbar(im1, ax=ax[0], pad=0.05, shrink=0.7, label=\"DN/s\")\n", - "\n", - " # Plot the mask: values are 1 or 0.\n", - " ax[1].set_ylabel(\"Pixel Row\", fontsize=12)\n", - " ax[1].set_xlabel(\"Pixel Column\", fontsize=12)\n", - " im2 = ax[1].imshow(mask, cmap=cmap, aspect=aspect, origin=\"lower\", vmin=0, vmax=1)\n", - " cbar2 = fig.colorbar(im2, ax=ax[1], pad=0.05, shrink=0.7, ticks=[0, 1])\n", - " cbar2.set_ticklabels([\"Illuminated Pixel Masked\", \"Un-Illuminated Pixel\"])\n", - "\n", - " # Plot the dark data with the same limits as the science data.\n", - " ax[2].set_ylabel(\"Pixel Row\", fontsize=12)\n", - " ax[2].set_xlabel(\"Pixel Column\", fontsize=12)\n", - " im3 = ax[2].imshow(\n", - " masked_sci, cmap=cmap, aspect=aspect, origin=\"lower\", vmin=vmin, vmax=vmax\n", - " )\n", - " cbar3 = fig.colorbar(im3, ax=ax[2], pad=0.05, shrink=0.7, label=\"DN/s\")\n", - "\n", - " if layout == \"rows\": # Adjust the multiplier as needed.\n", - " cbar1.ax.figure.set_size_inches(cbar1.ax.figure.get_size_inches() * 1.2)\n", - " cbar2.ax.figure.set_size_inches(cbar2.ax.figure.get_size_inches() * 1.2)\n", - " cbar3.ax.figure.set_size_inches(cbar3.ax.figure.get_size_inches() * 1.2)\n", - "\n", - " fig.tight_layout() # Adjusted tight_layout for better suptitle spacing." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "acfb32d8-fa5d-438b-9d04-58fd2c78dde9", - "metadata": {}, - "outputs": [], - "source": [ - "# Define a helper function for plotting the cleaned data.\n", - "\n", - "def plot_cleaned_data(\n", - " rate_file,\n", - " cleaned_file,\n", - " slice_indx=1,\n", - " aspect=\"auto\",\n", - " cmap=\"viridis\",\n", - " scale=0.2,\n", - " layout=\"rows\",\n", - "):\n", - " \"\"\"\n", - " Plot the 2D cleaned rate data (or 2D slice from a 3D cleaned data cube)\n", - " and compare against the original data (not applying NSClean).\n", - "\n", - " Parameters:\n", - " ----------\n", - " rate_file: str\n", - " File path to the FITS file containing original countrate data.\n", - " cleaned_file: str\n", - " File path to the FITS file containing cleaned countrate data.\n", - " slice_indx: int\n", - " Index of the slice to be plotted, 2D data default is 1.\n", - " aspect: int or 'auto'\n", - " Plot's aspect ratio.\n", - " cmap: str\n", - " Color map.\n", - " layout: str\n", - " Layout subplots in rows or columns?\n", - " scale: float\n", - " Scaling factor applied to determine the intensity range for visualization.\n", - " \"\"\"\n", - "\n", - " # Make a matplotlib figure.\n", - " if layout == \"rows\":\n", - " fig, ax = plt.subplots(4, 1, figsize=(10, 12))\n", - " else:\n", - " fig, ax = plt.subplots(2, 2, figsize=(12, 10))\n", - " ax = ax.flatten()\n", - " plt.suptitle(\n", - " f\"Cleaned Data (1/f noise removed) for {os.path.basename(rate_file)}\",\n", - " fontsize=15,\n", - " y=1,\n", - " )\n", - "\n", - " # Open the original and cleaned data.\n", - " with fits.open(rate_file) as hdulist:\n", - " original = hdulist[\"SCI\"].data\n", - "\n", - " with fits.open(cleaned_file) as hdulist:\n", - " cleaned = hdulist[\"SCI\"].data\n", - " # Determine if rate data is 2D or 3D.\n", - " if len(cleaned.shape) == 3:\n", - " dim = 3\n", - " else:\n", - " dim = 2\n", - "\n", - " # Define image limits from the original data.\n", - " vmin = np.nanpercentile(original, scale)\n", - " vmax = np.nanpercentile(original, 100 - scale)\n", - "\n", - " original[np.isnan(original)] = 0\n", - " cleaned[np.isnan(cleaned)] = 0\n", - "\n", - " # If the rate data is 3D, determine the integration (slice) to plot.\n", - " if dim == 3:\n", - " original = original[slice_indx, :, :]\n", - " cleaned = cleaned[slice_indx, :, :]\n", - " ax[0].set_title(\n", - " \"Original Rate Data: Integration [{},:,:]\".format(slice_indx), fontsize=15\n", - " )\n", - " ax[1].set_title(\n", - " \"Cleaned Rate Data: Integration [{},:,:]\".format(slice_indx), fontsize=15\n", - " )\n", - "\n", - " else:\n", - " ax[0].set_title(\"Original Rate Data\", fontsize=12)\n", - " ax[1].set_title(\"Cleaned Rate Data\", fontsize=12)\n", - "\n", - " # Plot the original rate data.\n", - " ax[0].set_xlabel(\"Pixel Column\", fontsize=10)\n", - " ax[0].set_ylabel(\"Pixel Row\", fontsize=10)\n", - " fig.colorbar(\n", - " ax[0].imshow(\n", - " original, cmap=cmap, origin=\"lower\", aspect=aspect, vmin=vmin, vmax=vmax\n", - " ),\n", - " ax=ax[0],\n", - " pad=0.05,\n", - " shrink=0.7,\n", - " label=\"DN/s\",\n", - " )\n", - "\n", - " # Plot the cleaned data with the same image limits.\n", - " ax[1].set_xlabel(\"Pixel Column\", fontsize=10)\n", - " ax[1].set_ylabel(\"Pixel Row\", fontsize=10)\n", - " fig.colorbar(\n", - " ax[1].imshow(\n", - " cleaned, cmap=cmap, origin=\"lower\", aspect=aspect, vmin=vmin, vmax=vmax\n", - " ),\n", - " ax=ax[1],\n", - " pad=0.05,\n", - " shrink=0.7,\n", - " label=\"DN/s\",\n", - " )\n", - "\n", - " # Plot the relative difference between the original and cleaned data.\n", - " diff = (original - cleaned) / original\n", - " diff[~np.isfinite(diff)] = 0\n", - " vmin_diff = np.nanpercentile(diff, scale)\n", - " vmax_diff = np.nanpercentile(diff, 100 - scale)\n", - " ax[2].set_title(\"Relative Difference (Original - Cleaned Rate Data)\", fontsize=12)\n", - " ax[2].set_xlabel(\"Pixel Column\", fontsize=10)\n", - " ax[2].set_ylabel(\"Pixel Row\", fontsize=10)\n", - " fig.colorbar(\n", - " ax[2].imshow(\n", - " diff,\n", - " cmap=cmap,\n", - " origin=\"lower\",\n", - " aspect=aspect,\n", - " vmin=vmin_diff,\n", - " vmax=vmax_diff,\n", - " ),\n", - " ax=ax[2],\n", - " pad=0.05,\n", - " shrink=0.7,\n", - " label=\"DN/s\",\n", - " )\n", - "\n", - " # Plot the relative difference again,\n", - " # this time hiding the outliers so that low-level\n", - " # background changes can be seen.\n", - " hide_outliers = np.ma.filled(sigma_clip(diff, masked=True), fill_value=0)\n", - " vmin_outliers = np.nanpercentile(hide_outliers, scale)\n", - " vmax_outliers = np.nanpercentile(hide_outliers, 100 - scale)\n", - " ax[3].set_title(\n", - " \"Relative Difference (Original - Cleaned Rate Data) \\n with Outliers Hidden\",\n", - " fontsize=12,\n", - " )\n", - " ax[3].set_xlabel(\"Pixel Column\", fontsize=10)\n", - " ax[3].set_ylabel(\"Pixel Row\", fontsize=10)\n", - " fig.colorbar(\n", - " ax[3].imshow(\n", - " hide_outliers,\n", - " cmap=cmap,\n", - " origin=\"lower\",\n", - " aspect=aspect,\n", - " vmin=vmin_outliers,\n", - " vmax=vmax_outliers,\n", - " ),\n", - " ax=ax[3],\n", - " pad=0.05,\n", - " shrink=0.7,\n", - " label=\"DN/s\",\n", - " )\n", - "\n", - " fig.tight_layout()" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "92e9da63-3814-4c05-ae66-34e72b9b5e7b", - "metadata": {}, - "outputs": [], - "source": [ - "# Define a helper function for plotting the 1D spectra.\n", - "\n", - "def plot_spectra(\n", - " spec_list,\n", - " ext_num=1,\n", - " scale_percent=2.0,\n", - " wavelength_range=None,\n", - " flux_range=None,\n", - " xlim_low=None,\n", - "):\n", - " \"\"\"\n", - " Plot the cleaned extracted 1D spectrum and compare against\n", - " the original (not applying NSClean) extracted 1D spectrum.\n", - "\n", - " Parameters:\n", - " ----------\n", - " spec_list: list\n", - " List of paths to the FITS files containing original/cleaned 1D extracted data.\n", - " scale_percent: float\n", - " Scaling factor applied to determine the intensity range for visualization.\n", - " ext_num: int\n", - " Index/extension of the slice to be plotted.\n", - " The EXTVER header value. The default is 1 for 2D data.\n", - " wavelength_range: dict\n", - " Wavelength range (x-axis) {'nrs1': [3.6, 3.65], 'nrs2': [1.65, 1.75]}.\n", - " flux_range: dict\n", - " Flux range (y-axis) {'nrs1': [1, 2], 'nrs2': [1, 2]}.\n", - " xlim_low: int\n", - " Define a lower wavelength end for the 1D spectrum. Helpful for BOTS data.\n", - " \"\"\"\n", - "\n", - " if wavelength_range is None:\n", - " wavelength_range = {}\n", - "\n", - " # Open the FITS files.\n", - " original_hdul = fits.open(spec_list[0])\n", - " cleaned_hdul = fits.open(spec_list[1])\n", - " if len(spec_list) == 3:\n", - " alternate_hdul = fits.open(spec_list[2])\n", - " elif len(spec_list) == 4:\n", - " alternate_hdul = fits.open(spec_list[2])\n", - " handmod_hdul = fits.open(spec_list[3])\n", - " else:\n", - " alternate_hdul = None\n", - " handmod_hdul = None\n", - "\n", - " # Find the spectral extension (EXTRACT1D).\n", - " for extnum in range(len(original_hdul)):\n", - " hdu = original_hdul[extnum]\n", - " if hdu.name != \"EXTRACT1D\":\n", - " continue\n", - " slit_name = hdu.header[\"SLTNAME\"]\n", - "\n", - " if hdu.header[\"EXTVER\"] == ext_num:\n", - "\n", - " # Plot the original and cleaned spectra together.\n", - " fig, ax = plt.subplots(1, 2, figsize=(10, 5))\n", - " plt.suptitle(\n", - " f\"Compare 1D Spectra for {os.path.basename(spec_list[0])};\"\n", - " f\"EXP_TYPE/Slit = {slit_name}\",\n", - " fontsize=15,\n", - " )\n", - "\n", - " if \"nrs1\" in spec_list[0] and xlim_low is not None:\n", - " xlim_low = xlim_low\n", - " else:\n", - " xlim_low = 0\n", - "\n", - " ax[0].step(\n", - " hdu.data[\"WAVELENGTH\"][xlim_low:],\n", - " hdu.data[\"FLUX\"][xlim_low:],\n", - " linewidth=1,\n", - " label=\"Original\",\n", - " )\n", - " ax[0].step(\n", - " cleaned_hdul[extnum].data[\"WAVELENGTH\"][xlim_low:],\n", - " cleaned_hdul[extnum].data[\"FLUX\"][xlim_low:],\n", - " linewidth=1,\n", - " label=\"Cleaned\",\n", - " )\n", - " ax[0].set_xlabel(f\"Wavelength ({hdu.header['TUNIT1']})\", fontsize=12)\n", - " ax[0].set_ylabel(f\"Flux ({hdu.header['TUNIT2']})\", fontsize=12)\n", - " ax[0].set_title(\"1D Extracted Spectra\", fontsize=12)\n", - "\n", - " # Plot the difference between the spectra as a ratio.\n", - " diff = (\n", - " cleaned_hdul[extnum].data[\"FLUX\"][xlim_low:]\n", - " / hdu.data[\"FLUX\"][xlim_low:]\n", - " )\n", - " ax[1].step(\n", - " hdu.data[\"WAVELENGTH\"][xlim_low:],\n", - " diff,\n", - " linewidth=1,\n", - " label=\"Cleaned/Original\",\n", - " )\n", - " ax[1].set_xlabel(f\"Wavelength ({hdu.header['TUNIT1']})\", fontsize=12)\n", - " ax[1].set_ylabel(\"Cleaned/Original\", fontsize=12)\n", - " ax[1].set_title(\"Difference After Cleaning\", fontsize=12)\n", - " # print(\"Average Difference = {}\".format(np.nanmean(diff)))\n", - "\n", - " # If available, also plot the alternate spectra.\n", - " if alternate_hdul is not None:\n", - " ax[0].step(\n", - " alternate_hdul[extnum].data[\"WAVELENGTH\"][xlim_low:],\n", - " alternate_hdul[extnum].data[\"FLUX\"][xlim_low:],\n", - " linewidth=1,\n", - " label=\"Cleaned (alternate mask)\",\n", - " )\n", - " diff2 = (\n", - " alternate_hdul[extnum].data[\"FLUX\"][xlim_low:]\n", - " / hdu.data[\"FLUX\"][xlim_low:]\n", - " )\n", - " ax[1].step(\n", - " hdu.data[\"WAVELENGTH\"][xlim_low:],\n", - " diff2,\n", - " linewidth=1,\n", - " label=\"Cleaned (alternate mask)/Original\",\n", - " )\n", - " if handmod_hdul is not None:\n", - "\n", - " ax[0].step(\n", - " handmod_hdul[extnum].data[\"WAVELENGTH\"][xlim_low:],\n", - " handmod_hdul[extnum].data[\"FLUX\"][xlim_low:],\n", - " linewidth=1,\n", - " label=\"Cleaned (hand-modified mask)\",\n", - " color=\"Purple\"\n", - " )\n", - " diff3 = (\n", - " handmod_hdul[extnum].data[\"FLUX\"][xlim_low:]\n", - " / hdu.data[\"FLUX\"][xlim_low:]\n", - " )\n", - " ax[1].step(\n", - " hdu.data[\"WAVELENGTH\"][xlim_low:],\n", - " diff3,\n", - " linewidth=1,\n", - " label=\"Cleaned (hand-modified mask)/Original\",\n", - " )\n", - "\n", - " # Set the y-range of the plot if needed.\n", - " if flux_range is not None:\n", - " for key, y_range in flux_range.items():\n", - " if key.lower() in spec_list[0].lower() and y_range is not None:\n", - " if len(y_range) == 2:\n", - " ax[0].set_ylim(y_range)\n", - " ax[1].set_ylim(\n", - " [\n", - " np.nanpercentile(diff, scale_percent),\n", - " np.nanpercentile(diff, 100 - scale_percent),\n", - " ]\n", - " )\n", - " # ax[1].set_ylim(0.5, 1.5)\n", - "\n", - " else:\n", - " all_flux = [\n", - " hdu.data[\"FLUX\"],\n", - " cleaned_hdul[extnum].data[\"FLUX\"],\n", - " ]\n", - " y_range_ax0 = [\n", - " np.nanpercentile(all_flux[0], scale_percent),\n", - " np.nanpercentile(all_flux[0], 100 - scale_percent),\n", - " ]\n", - " ax[0].set_ylim(y_range_ax0)\n", - " # ax[1].set_ylim(0.5, 1.5)\n", - "\n", - " ax[1].set_ylim(\n", - " [\n", - " np.nanpercentile(diff, scale_percent),\n", - " np.nanpercentile(diff, 100 - scale_percent),\n", - " ]\n", - " )\n", - "\n", - " else:\n", - " all_flux = [hdu.data[\"FLUX\"], cleaned_hdul[extnum].data[\"FLUX\"]]\n", - " y_range_ax0 = [\n", - " np.nanpercentile(all_flux[0], scale_percent),\n", - " np.nanpercentile(all_flux[0], 100 - scale_percent),\n", - " ]\n", - " ax[0].set_ylim(y_range_ax0)\n", - " # ax[1].set_ylim(0.5, 1.5)\n", - "\n", - " ax[1].set_ylim(\n", - " [\n", - " np.nanpercentile(diff, scale_percent),\n", - " np.nanpercentile(diff, 100 - scale_percent),\n", - " ]\n", - " )\n", - "\n", - " # Set the x-range of the plot if needed.\n", - " for key in wavelength_range:\n", - " if key in spec_list[0] and wavelength_range[key] is not None:\n", - " ax[0].set_xlim(wavelength_range[key])\n", - " ax[1].set_xlim(wavelength_range[key])\n", - "\n", - " fig.tight_layout()\n", - "\n", - " ax[0].legend()\n", - " ax[1].legend()\n", - "\n", - " original_hdul.close()\n", - " cleaned_hdul.close()\n", - " if alternate_hdul is not None:\n", - " alternate_hdul.close()\n", - " handmod_hdul.close()" - ] - }, { "cell_type": "markdown", "id": "b557c042-77fb-4c76-99bc-6921a24a9a0c", "metadata": {}, "source": [ - "## 4. Download the Data \n", + "## 3. Download the Data \n", "
\n", " \n", "The input data for this notebook features a MOS observation from the CEERS program with grating/filter G395H/F290LP. The dataset is part of the [JWST Early Release Science program ERS-1345](https://www.stsci.edu/jwst/science-execution/approved-programs/dd-ers/program-1345), specifically observation 61. It consists of 3 integrations (3 dither points; 3-SHUTTER-SLITLET) with 9 groups each. This notebook focuses on the third dithered exposure (00003) as an example. However, it's important to note that before proceeding to the `Spec3Pipeline`, all exposures must first be processed through the `Spec2Pipeline`." @@ -727,7 +169,7 @@ "id": "a7af0dc0-4c21-41df-902d-8c4f179a7719", "metadata": {}, "source": [ - "## 5. Running `Spec2Pipeline` without NSClean (Original Data) \n", + "## 4. Running `Spec2Pipeline` without NSClean (Original Data) \n", "
\n", "\n", "The cell below executes the `Spec2Pipeline`, explicitly skipping the NSClean step during processing. The level 2 products generated will serve as a reference point to illustrate how the countrate images and final extracted spectra appear without the removal of 1/f noise." @@ -786,7 +228,7 @@ "id": "5a42d053-0288-4ed6-a962-a265b27ec7d6", "metadata": {}, "source": [ - "## 6. Clean up 1/f Noise with NSClean (Default Pipeline Mask) \n", + "## 5. Clean up 1/f Noise with NSClean (Default Pipeline Mask) \n", "
\n", "\n", "If a user-supplied mask file is not provided to the [NSClean step](https://jwst-pipeline.readthedocs.io/en/latest/jwst/nsclean/index.html) in the `Spec2Pipeline`, the pipeline will generate a mask based on default parameters. This mask will identify any pixel that is unilluminated. That is, the mask must contain True and False values, where True indicates that the pixel is dark, and False indicates that the pixel is illuminated (not dark).\n", @@ -871,7 +313,7 @@ "id": "3d9867fb-d83a-4cad-b95d-078286dce27a", "metadata": {}, "source": [ - "### 6.1 Verify the Mask (Default Pipeline Mask) \n", + "### 5.1 Verify the Mask (Default Pipeline Mask) \n", "
\n", "\n", "Check the mask against the rate data to make sure it keeps only dark areas of the detector.\n", @@ -904,7 +346,7 @@ "id": "58d0ce26-dbfb-4391-af9b-9eb206a76daa", "metadata": {}, "source": [ - "### 6.2 Comparing Original vs. Cleaned Data (Default Pipeline Mask) \n", + "### 5.2 Comparing Original vs. Cleaned Data (Default Pipeline Mask) \n", "
\n", "\n", "We can now compare the cleaned data (with the default pipeline mask) to the original rate file and verify that the 1/f noise has been reduced.\n", @@ -1050,7 +492,7 @@ "id": "dae7529b-587e-490d-8c1b-69568e472501", "metadata": {}, "source": [ - "## 7. Clean up 1/f Noise with NSClean (Alternate Mask) \n", + "## 6. Clean up 1/f Noise with NSClean (Alternate Mask) \n", "
\n", "\n", "For this data set, note that some regions of the detector are heavily masked, due to overlapping slit regions. For example, see the cleaned rate data for NRS2 above, around x=1000, y=600. \n", @@ -1131,7 +573,7 @@ "id": "dadf08d6-e425-4a79-be60-e6d7c94fb780", "metadata": {}, "source": [ - "### 7.1 Verify the Mask (Alternate Mask) \n", + "### 6.1 Verify the Mask (Alternate Mask) \n", "
\n", "\n", "Check the mask against the rate data to make sure it keeps only dark areas of the detector.\n" @@ -1162,7 +604,7 @@ "id": "35e3f088-f1d9-484c-8d8d-3ec056c52152", "metadata": {}, "source": [ - "### 7.2 Comparing Original vs. Cleaned Data (Alternate Mask) \n", + "### 6.2 Comparing Original vs. Cleaned Data (Alternate Mask) \n", "
" ] }, @@ -1235,7 +677,7 @@ "id": "f727eada", "metadata": {}, "source": [ - "## 8. Clean up 1/f Noise with NSClean (Hand-Modified Mask) \n", + "## 7. Clean up 1/f Noise with NSClean (Hand-Modified Mask) \n", "
\n", "In certain scenarios, manual generation of a mask may be required. Here, we present **one** approach to manually modify the mask (editing an overly masked region in NRS2 around x=1000, y=600 to include more background; excluding some large snowballs in NRS1), starting with the default mask output from the pipeline. It is worth noting that the mask modified using this method may not necessarily outperform the two previous options.\n" ] @@ -1307,7 +749,7 @@ "id": "91f5cf29", "metadata": {}, "source": [ - "### 8.1 Verify the Mask (Hand-Modified Mask) \n", + "### 7.1 Verify the Mask (Hand-Modified Mask) \n", "
\n", "\n", "Check the mask against the rate data to make sure it keeps only dark areas of the detector.\n" @@ -1400,7 +842,7 @@ "id": "6ba19f8a", "metadata": {}, "source": [ - "### 8.2 Comparing Original vs. Cleaned Data (Hand-Modified Mask) \n", + "### 7.2 Comparing Original vs. Cleaned Data (Hand-Modified Mask) \n", "
" ] }, @@ -1469,7 +911,7 @@ "id": "797f0c0d-e694-4103-8c4d-7ad3a3d84eaf", "metadata": {}, "source": [ - "## 9. Conclusion \n", + "## 8. Conclusion \n", "
\n", "\n", "The final plots below show the countrate images and the resulting 1D extracted spectra side-by-side to compare the different cleaning methods: the original (no NSClean applied), the cleaned countrate image (with the default pipeline mask), the cleaned countrate image (with an alternate pipeline mask), and finally, the cleaned countrate image (with the hand-modified mask).\n", diff --git a/notebooks/NIRSpec_NSClean/utils.py b/notebooks/NIRSpec_NSClean/utils.py new file mode 100644 index 000000000..e2d0d8f26 --- /dev/null +++ b/notebooks/NIRSpec_NSClean/utils.py @@ -0,0 +1,525 @@ +import os +import numpy as np +import requests + +# ------ Plotting/Stats Imports ------ +from matplotlib import pyplot as plt +from astropy.stats import sigma_clip +from astropy.io import fits + + +def get_jwst_file(name, mast_api_token=None, save_directory=".", redownload=False): + """ + Retrieve a JWST data file from MAST archive and save it to a specified directory. + + Parameters: + ---------- + name: str + File name. + mast_api_token: str + MAST authorization token. Get your MAST Token Here: https://auth.mast.stsci.edu/token. + save_directory: str + Save directory path. + redownload: bool + Redownload the data even if it exsits already? + """ + mast_url = "https://mast.stsci.edu/api/v0.1/Download/file" + params = dict(uri=f"mast:JWST/product/{name}") + + if mast_api_token: + headers = dict(Authorization=f"token {mast_api_token}") + else: + headers = {} + + file_path = os.path.join(save_directory, name) + + # Check if the file already exists in the save directory. + if os.path.exists(file_path) and not redownload: + print(f"The file {name} already exists in the directory. Skipping download.") + return + + r = requests.get(mast_url, params=params, headers=headers, stream=True) + r.raise_for_status() + + with open(file_path, "wb") as fobj: + for chunk in r.iter_content(chunk_size=1024000): + fobj.write(chunk) + + +# Define a helper function for plotting the dark areas (the mask). +def plot_dark_data( + rate_file, + mask_file, + cmap="viridis", + scale=0.2, + slice_index=1, + axis=0, + aspect="auto", + layout="rows", +): + """ + Plot the dark areas on the detector, masking out the illuminated areas. + This function can take 2D (_rate.fits) and 3D (_rateints.fits) files as input. + + Parameters: + ---------- + rate_file: str + File path to the FITS file containing countrate data. + mask_file: str + File path to the FITS file containing mask data. + slice_index: int + Index of the slice to be plotted, 2D data default is 1. + axis: int (0-2) + Axis along which the slice will be taken + (0 for x-axis, 1 for y-axis, 2 for z-axis). + aspect: int or 'auto' + Plot's aspect ratio. + cmap: str + Color map. + layout: str + Layout subplots in rows or columns? + scale: float + Scaling factor applied to determine the intensity range for visualization. + """ + + # Make a matplotlib figure. + if layout == "rows": + fig, ax = plt.subplots(3, 1, figsize=(10, 5)) + else: + fig, ax = plt.subplots(1, 3, figsize=(20, 5)) + plt.suptitle(f"Dark Areas for {os.path.basename(rate_file)}", fontsize=15) + + # Open the mask and science data. + with fits.open(mask_file) as hdulist: + mask = hdulist[1].data + + with fits.open(rate_file) as hdulist: + sci = hdulist["SCI"].data + + # Determine if countrate data is 2D or 3D. + if len(sci.shape) == 3: + dim = 3 + else: + dim = 2 + + # Get data limits from the dark data. + masked_sci = sci.copy() + masked_sci[mask == 0] = 0 + vmin = np.nanpercentile(masked_sci, scale) + vmax = np.nanpercentile(masked_sci, 100 - scale) + + # Plot the science image with limits from the dark data. + sci[np.isnan(sci)] = 0 + + # If the countrate data is 3D, determine the integration (slice) to plot. + if dim == 3: + sci = np.take(sci, slice_index, axis=axis) + mask = np.take(mask, slice_index, axis=axis) + masked_sci = np.take(masked_sci, slice_index, axis=axis) + + ax[0].set_title("Original Rate Data: Integration [{},:,:]".format(slice_index)) + ax[1].set_title( + "Illuminated Pixel Mask: Integration [{},:,:]".format(slice_index) + ) + ax[2].set_title( + "Dark Data (Illuminated Pixels Masked): Integration [{},:,:]".format( + slice_index + ) + ) + + else: + sci = sci + ax[0].set_title("Original Rate Data") + ax[1].set_title("Illuminated Pixel Mask") + ax[2].set_title("Dark Data (Illuminated Pixels Masked)") + + # Plot the science image with limits from the dark data. + ax[0].set_ylabel("Pixel Row", fontsize=12) + ax[0].set_xlabel("Pixel Column", fontsize=12) + im1 = ax[0].imshow( + sci, cmap=cmap, origin="lower", aspect=aspect, vmin=vmin, vmax=vmax + ) + cbar1 = fig.colorbar(im1, ax=ax[0], pad=0.05, shrink=0.7, label="DN/s") + + # Plot the mask: values are 1 or 0. + ax[1].set_ylabel("Pixel Row", fontsize=12) + ax[1].set_xlabel("Pixel Column", fontsize=12) + im2 = ax[1].imshow(mask, cmap=cmap, aspect=aspect, origin="lower", vmin=0, vmax=1) + cbar2 = fig.colorbar(im2, ax=ax[1], pad=0.05, shrink=0.7, ticks=[0, 1]) + cbar2.set_ticklabels(["Illuminated Pixel Masked", "Un-Illuminated Pixel"]) + + # Plot the dark data with the same limits as the science data. + ax[2].set_ylabel("Pixel Row", fontsize=12) + ax[2].set_xlabel("Pixel Column", fontsize=12) + im3 = ax[2].imshow( + masked_sci, cmap=cmap, aspect=aspect, origin="lower", vmin=vmin, vmax=vmax + ) + cbar3 = fig.colorbar(im3, ax=ax[2], pad=0.05, shrink=0.7, label="DN/s") + + if layout == "rows": # Adjust the multiplier as needed. + cbar1.ax.figure.set_size_inches(cbar1.ax.figure.get_size_inches() * 1.2) + cbar2.ax.figure.set_size_inches(cbar2.ax.figure.get_size_inches() * 1.2) + cbar3.ax.figure.set_size_inches(cbar3.ax.figure.get_size_inches() * 1.2) + + fig.tight_layout() # Adjusted tight_layout for better suptitle spacing. + + +# Define a helper function for plotting the cleaned data. +def plot_cleaned_data( + rate_file, + cleaned_file, + slice_index=1, + aspect="auto", + cmap="viridis", + scale=0.2, + layout="rows", +): + """ + Plot the 2D cleaned rate data (or 2D slice from a 3D cleaned data cube) + and compare against the original data (not applying NSClean). + + Parameters: + ---------- + rate_file: str + File path to the FITS file containing original countrate data. + cleaned_file: str + File path to the FITS file containing cleaned countrate data. + slice_index: int + Index of the slice to be plotted, 2D data default is 1. + aspect: int or 'auto' + Plot's aspect ratio. + cmap: str + Color map. + layout: str + Layout subplots in rows or columns? + scale: float + Scaling factor applied to determine the intensity range for visualization. + """ + + # Make a matplotlib figure. + if layout == "rows": + fig, ax = plt.subplots(4, 1, figsize=(10, 12)) + else: + fig, ax = plt.subplots(2, 2, figsize=(12, 10)) + ax = ax.flatten() + plt.suptitle( + f"Cleaned Data (1/f noise removed) for {os.path.basename(rate_file)}", + fontsize=15, + y=1, + ) + + # Open the original and cleaned data. + with fits.open(rate_file) as hdulist: + original = hdulist["SCI"].data + + with fits.open(cleaned_file) as hdulist: + cleaned = hdulist["SCI"].data + # Determine if rate data is 2D or 3D. + if len(cleaned.shape) == 3: + dim = 3 + else: + dim = 2 + + # Define image limits from the original data. + vmin = np.nanpercentile(original, scale) + vmax = np.nanpercentile(original, 100 - scale) + + original[np.isnan(original)] = 0 + cleaned[np.isnan(cleaned)] = 0 + + # If the rate data is 3D, determine the integration (slice) to plot. + if dim == 3: + original = original[slice_indx, :, :] + cleaned = cleaned[slice_indx, :, :] + ax[0].set_title( + "Original Rate Data: Integration [{},:,:]".format(slice_indx), fontsize=15 + ) + ax[1].set_title( + "Cleaned Rate Data: Integration [{},:,:]".format(slice_indx), fontsize=15 + ) + + else: + ax[0].set_title("Original Rate Data", fontsize=12) + ax[1].set_title("Cleaned Rate Data", fontsize=12) + + # Plot the original rate data. + ax[0].set_xlabel("Pixel Column", fontsize=10) + ax[0].set_ylabel("Pixel Row", fontsize=10) + fig.colorbar( + ax[0].imshow( + original, cmap=cmap, origin="lower", aspect=aspect, vmin=vmin, vmax=vmax + ), + ax=ax[0], + pad=0.05, + shrink=0.7, + label="DN/s", + ) + + # Plot the cleaned data with the same image limits. + ax[1].set_xlabel("Pixel Column", fontsize=10) + ax[1].set_ylabel("Pixel Row", fontsize=10) + fig.colorbar( + ax[1].imshow( + cleaned, cmap=cmap, origin="lower", aspect=aspect, vmin=vmin, vmax=vmax + ), + ax=ax[1], + pad=0.05, + shrink=0.7, + label="DN/s", + ) + + # Plot the relative difference between the original and cleaned data. + diff = (original - cleaned) / original + diff[~np.isfinite(diff)] = 0 + vmin_diff = np.nanpercentile(diff, scale) + vmax_diff = np.nanpercentile(diff, 100 - scale) + ax[2].set_title("Relative Difference (Original - Cleaned Rate Data)", fontsize=12) + ax[2].set_xlabel("Pixel Column", fontsize=10) + ax[2].set_ylabel("Pixel Row", fontsize=10) + fig.colorbar( + ax[2].imshow( + diff, + cmap=cmap, + origin="lower", + aspect=aspect, + vmin=vmin_diff, + vmax=vmax_diff, + ), + ax=ax[2], + pad=0.05, + shrink=0.7, + label="DN/s", + ) + + # Plot the relative difference again, + # this time hiding the outliers so that low-level + # background changes can be seen. + hide_outliers = np.ma.filled(sigma_clip(diff, masked=True), fill_value=0) + vmin_outliers = np.nanpercentile(hide_outliers, scale) + vmax_outliers = np.nanpercentile(hide_outliers, 100 - scale) + ax[3].set_title( + "Relative Difference (Original - Cleaned Rate Data) \n with Outliers Hidden", + fontsize=12, + ) + ax[3].set_xlabel("Pixel Column", fontsize=10) + ax[3].set_ylabel("Pixel Row", fontsize=10) + fig.colorbar( + ax[3].imshow( + hide_outliers, + cmap=cmap, + origin="lower", + aspect=aspect, + vmin=vmin_outliers, + vmax=vmax_outliers, + ), + ax=ax[3], + pad=0.05, + shrink=0.7, + label="DN/s", + ) + + fig.tight_layout() + + +# Define a helper function for plotting the 1D spectra. +def plot_spectra( + spec_list, + ext_num=1, + scale_percent=2.0, + wavelength_range=None, + flux_range=None, + xlim_low=None, +): + """ + Plot the cleaned extracted 1D spectrum and compare against + the original (not applying NSClean) extracted 1D spectrum. + + Parameters: + ---------- + spec_list: list + List of paths to the FITS files containing original/cleaned 1D extracted data. + scale_percent: float + Scaling factor applied to determine the intensity range for visualization. + ext_num: int + Index/extension of the slice to be plotted. + The EXTVER header value. The default is 1 for 2D data. + wavelength_range: dict + Wavelength range (x-axis) {'nrs1': [3.6, 3.65], 'nrs2': [1.65, 1.75]}. + flux_range: dict + Flux range (y-axis) {'nrs1': [1, 2], 'nrs2': [1, 2]}. + xlim_low: int + Define a lower wavelength end for the 1D spectrum. Helpful for BOTS data. + """ + + if wavelength_range is None: + wavelength_range = {} + + # Open the FITS files. + original_hdul = fits.open(spec_list[0]) + cleaned_hdul = fits.open(spec_list[1]) + if len(spec_list) == 3: + alternate_hdul = fits.open(spec_list[2]) + elif len(spec_list) == 4: + alternate_hdul = fits.open(spec_list[2]) + handmod_hdul = fits.open(spec_list[3]) + else: + alternate_hdul = None + handmod_hdul = None + + # Find the spectral extension (EXTRACT1D). + for extnum in range(len(original_hdul)): + hdu = original_hdul[extnum] + if hdu.name != "EXTRACT1D": + continue + slit_name = hdu.header["SLTNAME"] + + if hdu.header["EXTVER"] == ext_num: + + # Plot the original and cleaned spectra together. + fig, ax = plt.subplots(1, 2, figsize=(10, 5)) + plt.suptitle( + f"Compare 1D Spectra for {os.path.basename(spec_list[0])};" + f"EXP_TYPE/Slit = {slit_name}", + fontsize=15, + ) + + if "nrs1" in spec_list[0] and xlim_low is not None: + xlim_low = xlim_low + else: + xlim_low = 0 + + ax[0].step( + hdu.data["WAVELENGTH"][xlim_low:], + hdu.data["FLUX"][xlim_low:], + linewidth=1, + label="Original", + ) + ax[0].step( + cleaned_hdul[extnum].data["WAVELENGTH"][xlim_low:], + cleaned_hdul[extnum].data["FLUX"][xlim_low:], + linewidth=1, + label="Cleaned", + ) + ax[0].set_xlabel(f"Wavelength ({hdu.header['TUNIT1']})", fontsize=12) + ax[0].set_ylabel(f"Flux ({hdu.header['TUNIT2']})", fontsize=12) + ax[0].set_title("1D Extracted Spectra", fontsize=12) + + # Plot the difference between the spectra as a ratio. + diff = ( + cleaned_hdul[extnum].data["FLUX"][xlim_low:] + / hdu.data["FLUX"][xlim_low:] + ) + ax[1].step( + hdu.data["WAVELENGTH"][xlim_low:], + diff, + linewidth=1, + label="Cleaned/Original", + ) + ax[1].set_xlabel(f"Wavelength ({hdu.header['TUNIT1']})", fontsize=12) + ax[1].set_ylabel("Cleaned/Original", fontsize=12) + ax[1].set_title("Difference After Cleaning", fontsize=12) + # print("Average Difference = {}".format(np.nanmean(diff))) + + # If available, also plot the alternate spectra. + if alternate_hdul is not None: + ax[0].step( + alternate_hdul[extnum].data["WAVELENGTH"][xlim_low:], + alternate_hdul[extnum].data["FLUX"][xlim_low:], + linewidth=1, + label="Cleaned (alternate mask)", + ) + diff2 = ( + alternate_hdul[extnum].data["FLUX"][xlim_low:] + / hdu.data["FLUX"][xlim_low:] + ) + ax[1].step( + hdu.data["WAVELENGTH"][xlim_low:], + diff2, + linewidth=1, + label="Cleaned (alternate mask)/Original", + ) + if handmod_hdul is not None: + + ax[0].step( + handmod_hdul[extnum].data["WAVELENGTH"][xlim_low:], + handmod_hdul[extnum].data["FLUX"][xlim_low:], + linewidth=1, + label="Cleaned (hand-modified mask)", + color="Purple" + ) + diff3 = ( + handmod_hdul[extnum].data["FLUX"][xlim_low:] + / hdu.data["FLUX"][xlim_low:] + ) + ax[1].step( + hdu.data["WAVELENGTH"][xlim_low:], + diff3, + linewidth=1, + label="Cleaned (hand-modified mask)/Original", + ) + + # Set the y-range of the plot if needed. + if flux_range is not None: + for key, y_range in flux_range.items(): + if key.lower() in spec_list[0].lower() and y_range is not None: + if len(y_range) == 2: + ax[0].set_ylim(y_range) + ax[1].set_ylim( + [ + np.nanpercentile(diff, scale_percent), + np.nanpercentile(diff, 100 - scale_percent), + ] + ) + # ax[1].set_ylim(0.5, 1.5) + + else: + all_flux = [ + hdu.data["FLUX"], + cleaned_hdul[extnum].data["FLUX"], + ] + y_range_ax0 = [ + np.nanpercentile(all_flux[0], scale_percent), + np.nanpercentile(all_flux[0], 100 - scale_percent), + ] + ax[0].set_ylim(y_range_ax0) + # ax[1].set_ylim(0.5, 1.5) + + ax[1].set_ylim( + [ + np.nanpercentile(diff, scale_percent), + np.nanpercentile(diff, 100 - scale_percent), + ] + ) + + else: + all_flux = [hdu.data["FLUX"], cleaned_hdul[extnum].data["FLUX"]] + y_range_ax0 = [ + np.nanpercentile(all_flux[0], scale_percent), + np.nanpercentile(all_flux[0], 100 - scale_percent), + ] + ax[0].set_ylim(y_range_ax0) + # ax[1].set_ylim(0.5, 1.5) + + ax[1].set_ylim( + [ + np.nanpercentile(diff, scale_percent), + np.nanpercentile(diff, 100 - scale_percent), + ] + ) + + # Set the x-range of the plot if needed. + for key in wavelength_range: + if key in spec_list[0] and wavelength_range[key] is not None: + ax[0].set_xlim(wavelength_range[key]) + ax[1].set_xlim(wavelength_range[key]) + + fig.tight_layout() + + ax[0].legend() + ax[1].legend() + + original_hdul.close() + cleaned_hdul.close() + if alternate_hdul is not None: + alternate_hdul.close() + handmod_hdul.close()