Boxcar extraction using Trace class

notebook_sandbox/jwst_boxcar/boxcar_extraction.ipynb

@@ -0,0 +1,281 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "# Boxcar extraction with specreduce"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "This notebook demonstrates a simplified boxcar extraction scenario, using a rectified spectrum in a 2D spectral image where the dispersion axis is the second (X) axis. \n",
+    "\n",
+    "The extraction algorithm in `specreduce` is a plain adaptation of the algorithm presented in the MIRI LRS extraction notebook at https://github.com/spacetelescope/jdat_notebooks/tree/main/notebooks/MIRI_LRS_spectral_extraction This algorithm is demonstrated separately from the `specreduce` package, in the acompanying notebook http://localhost:8888/notebooks/notebook_sandbox/jwst_boxcar/jwst_boxcar_algorithm.ipynb\n",
+    "\n",
+    "Note that the extraction algorithm in the MIRI LRS extraction notebook performs simultaneous background subtraction when extracting from the source. Although it can only subtract the average of the two background extraction boxes. The original extraction code in `specreduce`'s `BoxcarExtract` class was capable of modeling the background with a polynomial along the cross-dispersion direction. This feature was lost in the conversion. \n",
+    "\n",
+    "Note also that we cannot yet perform offline, after-the-fact background subtractions from spectra extracted with `BoxcarExtract`. The reason is that `BoxcarExtract` delivers its `Spectrum1D` product with `spectral_axis` in units of `pixel`. Subsequent operations with such spectra, such as subtraction, are severely limited due to the `pixel` units."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "slideshow": {
+     "slide_type": "fragment"
+    }
+   },
+   "outputs": [],
+   "source": [
+    "from matplotlib.colors import LogNorm\n",
+    "%matplotlib inline\n",
+    "\n",
+    "from pathlib import Path\n",
+    "from zipfile import ZipFile\n",
+    "\n",
+    "from astropy.io import fits\n",
+    "from astropy.table import Table\n",
+    "from astropy.visualization import simple_norm\n",
+    "from astropy.utils.data import download_file\n",
+    "\n",
+    "from jwst import datamodels\n",
+    "\n",
+    "from specreduce.extract import BoxcarExtract\n",
+    "from specreduce.tracing import FlatTrace\n",
+    "\n",
+    "import os\n",
+    "import tempfile\n",
+    "\n",
+    "import numpy as np\n",
+    "import ccdproc\n",
+    "\n",
+    "import matplotlib.pyplot as plt\n",
+    "import matplotlib as mpl"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Ingest s2d data"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# data is taken from s2d file. x1d is used for comparison with pipeline extraction.\n",
+    "zipped_datapath = Path(download_file('https://stsci.box.com/shared/static/qdj79y7rv99wuui2hot0l3kg5ohq0ah9.zip', cache=True))\n",
+    "\n",
+    "data_dir = Path(tempfile.gettempdir())\n",
+    "\n",
+    "with ZipFile(zipped_datapath, 'r') as sample_data_zip:\n",
+    "    sample_data_zip.extractall(data_dir)\n",
+    "\n",
+    "s2dfile = str(data_dir / \"nirspec_fssim_d1_s2d.fits\")\n",
+    "x1dfile = str(data_dir / \"nirspec_fssim_d1_x1d.fits\")"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# use a jwst datamodel to provide a good interface to the data and wcs info\n",
+    "s2d = datamodels.open(s2dfile)\n",
+    "image = s2d.slits[0].data"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# display s2d image\n",
+    "norm_data = simple_norm(image, \"sqrt\")\n",
+    "plt.figure(figsize=(15, 15))\n",
+    "plt.imshow(image, norm=norm_data, origin=\"lower\")\n",
+    "plt.title(\"slit[0]\")"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "scrolled": false
+   },
+   "outputs": [],
+   "source": [
+    "# pipeline 1d extraction (for comparison)\n",
+    "jpipe_x1d = Table.read(x1dfile, hdu=1)\n",
+    "print(jpipe_x1d.columns)\n",
+    "# plot\n",
+    "fig, ax = plt.subplots(figsize=(10, 6))\n",
+    "ax.plot(jpipe_x1d['WAVELENGTH'], jpipe_x1d['FLUX'], 'k-', label=\"jpipe_x1d\")\n",
+    "ax.set_title(\"JWST Pipeline x1d extracted spectrum\")\n",
+    "ax.set_xlabel(\"wavelength\")\n",
+    "ax.set_ylabel(\"Flux Density [Jy]\")\n",
+    "ax.set_yscale(\"log\")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Define region to be extracted"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "At this point, `specreduce` doesn't provide a tool for finding the trace. Since this is rectified data, we can use the `FlatTrace` class and find the spectrum position and width by eye."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# blow up of the region to be extracted\n",
+    "plt.figure(figsize=(15, 15))\n",
+    "plt.imshow(image[::,0:100], norm=norm_data, origin=\"lower\")\n",
+    "plt.title(\"slit[0] slice\")"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# extraction parameters based on image above\n",
+    "ext_center = 27\n",
+    "ext_width = 4"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Plot along cross-disperion cut showing the extraction parameters\n",
+    "fig, ax = plt.subplots(figsize=(10, 6))\n",
+    "y = np.arange(image.shape[0])\n",
+    "ax.plot(y, image[:,70], 'k-')\n",
+    "mm = np.array([ext_center, ext_center])\n",
+    "mm_y = ax.get_ylim()\n",
+    "\n",
+    "ax.plot(mm, mm_y, 'b--')\n",
+    "ax.plot(mm - ext_width/2., mm_y, 'g:')\n",
+    "ax.plot(mm + ext_width/2., mm_y, 'g:')\n",
+    "\n",
+    "ax.set_title(\"Cross-dispersion Cut at Pixel=70\")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Extract"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "scrolled": false
+   },
+   "outputs": [],
+   "source": [
+    "# define the Trace\n",
+    "trace = FlatTrace(image, ext_center)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# extract\n",
+    "boxcar = BoxcarExtract()\n",
+    "spectrum = boxcar(image, trace, width=ext_width)"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# plot \n",
+    "f, ax = plt.subplots(figsize=(10, 6))\n",
+    "ax.plot(spectrum.flux.value, 'k-')\n",
+    "ax.set_title(\"Boxcar 1D extracted spectrum\")\n",
+    "ax.set_xlabel(r\"pixel\")\n",
+    "ax.set_ylabel(\"Flux Density [Jy]\")\n",
+    "ax.set_yscale(\"log\")"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {
+    "slideshow": {
+     "slide_type": "slide"
+    }
+   },
+   "source": [
+    "## About this notebook\n",
+    "\n",
+    "**Author:** Ivo Busko, JWST\n",
+    "**Updated On:** 2022-02-18"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "***"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "[Top of Page](#top)\n",
+    "<img style=\"float: right;\" src=\"https://raw.githubusercontent.com/spacetelescope/notebooks/master/assets/stsci_pri_combo_mark_horizonal_white_bkgd.png\" alt=\"Space Telescope Logo\" width=\"200px\"/> "
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3 (ipykernel)",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.8.12"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 2
+}