From 021477a04a6cc6ce7f749246e817f069d1f90b70 Mon Sep 17 00:00:00 2001 From: haticekaratay Date: Fri, 18 Aug 2023 14:06:37 -0400 Subject: [PATCH 1/2] Add MIRI lrs notebook --- _toc.yml | 2 +- index.rst | 8 ++++---- 2 files changed, 5 insertions(+), 5 deletions(-) diff --git a/_toc.yml b/_toc.yml index 3c84ce2ae..cd0c5e33d 100644 --- a/_toc.yml +++ b/_toc.yml @@ -69,7 +69,7 @@ parts: title: MRS Mstar - Data Analysis - file: notebooks/MIRI_IFU_YSOs_in_the_LMC/isha_nayak_ysos_in_the_lmc.ipynb title: IFU of YSOs in LMC - - file: notebooks/MIRI_LRS_spectral_extraction/miri_lrs_spectral_extraction.ipynb + - file: notebooks/MIRI_LRS_spectral_extraction/miri_lrs_advanced_extraction_part1.ipynb title: LRS Optimal Spectral Extraction - caption: NIRCam chapters: diff --git a/index.rst b/index.rst index 2a4e5f2eb..cf9514208 100644 --- a/index.rst +++ b/index.rst @@ -87,10 +87,10 @@ location of the notebooks for you to :ref:`download `. - **Tools:** specutils, photutils, astropy. - **Cross-instrument:** * - `LRS Optimal Extraction `_ - - - **Use case:**Optimal spectral extraction. - - **Data:** MIRISim simulated LRS spectrum. - - **Tools:** jwst pipeline, gwcs. - - **Cross-instrument:** NIRSpec, NIRISS, MIRI + - - **Use case:**Spectral extraction of slit spectra with the JWST calibration pipeline. + - **Data:** Publicly available science data + - **Tools:** jwst, matplotlib, astropy. + - **Cross-instrument:** NIRSpec, MIRI * - NIRCam - * - `Point Source Aperture Photometry `_ From 35e88852030b0f52c881238c21a523c70f79cc1e Mon Sep 17 00:00:00 2001 From: haticekaratay Date: Fri, 18 Aug 2023 14:26:23 -0400 Subject: [PATCH 2/2] Delete the old notebook --- .../miri_lrs_spectral_extraction.ipynb | 864 ------------------ 1 file changed, 864 deletions(-) delete mode 100644 notebooks/MIRI_LRS_spectral_extraction/miri_lrs_spectral_extraction.ipynb diff --git a/notebooks/MIRI_LRS_spectral_extraction/miri_lrs_spectral_extraction.ipynb b/notebooks/MIRI_LRS_spectral_extraction/miri_lrs_spectral_extraction.ipynb deleted file mode 100644 index 5a705310e..000000000 --- a/notebooks/MIRI_LRS_spectral_extraction/miri_lrs_spectral_extraction.ipynb +++ /dev/null @@ -1,864 +0,0 @@ -{ - "cells": [ - { - "cell_type": "markdown", - "metadata": { - "slideshow": { - "slide_type": "slide" - } - }, - "source": [ - "# LRS Optimal Spectral Extraction" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "**Use case:** Extract spectra with different locations, extraction apertures, and techniques.
\n", - "**Data:** Simulated MIRI LRS spectrum.
\n", - "**Tools:** jwst, gwcs, matplotlib, astropy.
\n", - "**Cross-intrument:** NIRSpec, MIRI.
\n", - "**Documentation:** This notebook is part of a STScI's larger [post-pipeline Data Analysis Tools Ecosystem](https://jwst-docs.stsci.edu/jwst-post-pipeline-data-analysis).
\n", - "\n", - "## Introduction\n", - "\n", - "This notebook extracts a 1D spectra from a 2D MIRI LRS spectral observation (single image). The goal is to provide the ability to extract spectra with different locations, extraction apertures, and techniques than are done in the JWST pipeline.\n", - "\n", - "The simpliest spectral extraction is \"boxcar\" where all the pixels within some fixed width centered on the source position are summed at each wavelength. Background subtraction can be done using regions offset from the source center.\n", - "\n", - "For spectra taken with a diffraction limited telescope like JWST, a modification boxcar extraction is to vary the extraction width linearly with wavelength. Such a scaled boxcar extraction keeps the fraction of the source flux within the extraction region approximately constant with wavelength.\n", - "\n", - "For point sources, a PSF-weighted spectral extraction can be done. Using the PSF to weight the extraction uses the actual PSF as a function of wavelength to optimize the extraction to the pixels with the greatest signal. PSF-weighted extractions show the largest differences with boxcar extractions at lower S/N values." - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "**Note:** Corrections for the finite aperture used in all the extractions have not been applied. Thus, the physical flux densities of all the extracted spectra are lower than the actual values." - ] - }, - { - "cell_type": "markdown", - "metadata": { - "slideshow": { - "slide_type": "slide" - } - }, - "source": [ - "## Imports\n", - "\n", - "- *matplotlib.pyplot* for plotting data\n", - "- *numpy* to handle array functions\n", - "- *astropy.io fits* for accessing FITS files\n", - "- *astropy.visualization* for scaling image for display\n", - "- *astropy.table Table* for reading the pipeline 1d extractions\n", - "- *jwst datamodels* for reading/access the jwst data" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": { - "slideshow": { - "slide_type": "fragment" - } - }, - "outputs": [], - "source": [ - "import matplotlib.pyplot as plt\n", - "import matplotlib as mpl\n", - "%matplotlib inline\n", - "\n", - "import numpy as np\n", - "\n", - "from astropy.io import fits\n", - "from astropy.table import Table\n", - "from astropy.visualization import simple_norm\n", - "\n", - "from jwst import datamodels" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# useful function that work for boxcar, boxcar scaled with wavelength,\n", - "# and psf-weighted extractions\n", - "import numpy as np\n", - "from gwcs.wcstools import grid_from_bounding_box\n", - "\n", - "\n", - "def get_boxcar_weights(center, hwidth, npix):\n", - " \"\"\"\n", - " Compute the weights given an aperture center, half widths, and number of pixels\n", - " \"\"\"\n", - " weights = np.zeros((npix))\n", - " # pixels with full weight\n", - " fullpixels = [max(0, int(center - hwidth + 1)), min(int(center + hwidth), npix)]\n", - " weights[fullpixels[0] : fullpixels[1]] = 1.0\n", - "\n", - " # pixels at the edges of the boxcar with partial weight\n", - " if fullpixels[0] > 0:\n", - " weights[fullpixels[0] - 1] = hwidth - (center - fullpixels[0])\n", - " if fullpixels[1] < npix:\n", - " weights[fullpixels[1]] = hwidth - (fullpixels[1] - center)\n", - "\n", - " return weights\n", - "\n", - "\n", - "def ap_weight_images(\n", - " center, width, bkg_offset, bkg_width, image_size, waves, wavescale=None\n", - "):\n", - " \"\"\"\n", - " Create a weight image that defines the desired extraction aperture\n", - " and the weight image for the requested background regions\n", - "\n", - " Parameters\n", - " ----------\n", - " center : float\n", - " center of aperture in pixels\n", - " width : float\n", - " width of apeture in pixels\n", - " bkg_offset : float\n", - " offset from the extaction edge for the background\n", - " never scaled for wavelength\n", - " bkg_width : float\n", - " width of background region\n", - " never scaled with wavelength\n", - " image_size : tuple with 2 elements\n", - " size of image\n", - " waves : array\n", - " wavelegth values\n", - " wavescale : float\n", - " scale the width with wavelength (default=None)\n", - " wavescale gives the reference wavelenth for the width value\n", - "\n", - " Returns\n", - " -------\n", - " wimage, bkg_wimage : (2D image, 2D image)\n", - " wimage is the weight image defining the aperature\n", - " bkg_image is the weight image defining the background regions\n", - " \"\"\"\n", - " wimage = np.zeros(image_size)\n", - " bkg_wimage = np.zeros(image_size)\n", - " hwidth = 0.5 * width\n", - " # loop in dispersion direction and compute weights\n", - " for i in range(image_size[1]):\n", - " if wavescale is not None:\n", - " hwidth = 0.5 * width * (waves[i] / wavescale)\n", - "\n", - " wimage[:, i] = get_boxcar_weights(center, hwidth, image_size[0])\n", - "\n", - " # bkg regions\n", - " if (bkg_width is not None) & (bkg_offset is not None):\n", - " bkg_wimage[:, i] = get_boxcar_weights(\n", - " center - hwidth - bkg_offset, bkg_width, image_size[0]\n", - " )\n", - " bkg_wimage[:, i] += get_boxcar_weights(\n", - " center + hwidth + bkg_offset, bkg_width, image_size[0]\n", - " )\n", - " else:\n", - " bkg_wimage = None\n", - "\n", - " return (wimage, bkg_wimage)\n", - "\n", - "\n", - "def extract_1dspec(jdatamodel, center, width, bkg_offset, bkg_width, wavescale=None):\n", - " \"\"\"\n", - " Extract the 1D spectrum using the boxcar method.\n", - " Does a background subtraction as part of the extraction.\n", - "\n", - " Parameters\n", - " ----------\n", - " jdatamodel : jwst.DataModel\n", - " jwst datamodel with the 2d spectral image\n", - " center : float\n", - " center of aperture in pixels\n", - " width : float\n", - " width of apeture in pixels\n", - " bkg_offset : float\n", - " offset from the extaction edge for the background\n", - " never scaled for wavelength\n", - " bkg_width : float\n", - " width of background region\n", - " never scaled with wavelength\n", - " wavescale : float\n", - " scale the width with wavelength (default=None)\n", - " wavescale gives the reference wavelenth for the width value\n", - "\n", - " Returns\n", - " -------\n", - " waves, ext1d : (ndarray, ndarray)\n", - " 2D `float` array with wavelengths\n", - " 1D `float` array with extracted 1d spectrum in Jy\n", - " \"\"\"\n", - " # should be determined from the gWCS in cal.fits\n", - " image = np.transpose(jdatamodel.data)\n", - " grid = grid_from_bounding_box(jdatamodel.meta.wcs.bounding_box)\n", - " ra, dec, lam = jdatamodel.meta.wcs(*grid)\n", - " lam_image = np.transpose(lam)\n", - "\n", - " # compute a \"rough\" wavelength scale to allow for aperture to scale with wavelength\n", - " rough_waves = np.average(lam_image, axis=0)\n", - "\n", - " # images to use for extraction\n", - " wimage, bkg_wimage = ap_weight_images(\n", - " center,\n", - " width,\n", - " bkg_width,\n", - " bkg_offset,\n", - " image.shape,\n", - " rough_waves,\n", - " wavescale=wavescale,\n", - " )\n", - "\n", - " # extract the spectrum using the weight image\n", - " if bkg_wimage is not None:\n", - " ext1d_boxcar_bkg = np.average(image, weights=bkg_wimage, axis=0)\n", - " data_bkgsub = image - np.tile(ext1d_boxcar_bkg, (image.shape[0], 1))\n", - " else:\n", - " data_bkgsub = image\n", - "\n", - " ext1d = np.sum(data_bkgsub * wimage, axis=0)\n", - " # convert from MJy/sr to Jy\n", - " ext1d *= 1e6 * jdatamodel.meta.photometry.pixelarea_steradians\n", - "\n", - " # compute the average wavelength for each column using the weight image\n", - " # this should correspond directly with the extracted spectrum\n", - " # wavelengths account for any tiled spectra this way\n", - " waves = np.average(lam_image, weights=wimage, axis=0)\n", - "\n", - " return (waves, ext1d, data_bkgsub)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### Devloper notes\n", - "\n", - "The difference between the pipeline (x1d) and the extractions done in this notebook are quite large. Help in understanding the origin of these differences is needed.\n", - "\n", - "Not clear how to use the JWST pipeline `extract_1d` (quite complex) code.\n", - "Help to determine how to use the JWST pipeline code instead of the custom code for boxcar is needed. \n", - "\n", - "Applying aperture corrections for the finite extraction widths is needed. Help in how to get the needed informatinom for different (user set) extraction widths is needed. " - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## Download Files" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "from astropy.utils.data import download_file\n", - "\n", - "calfilename = \"det_image_seq5_MIRIMAGE_P750Lexp1_cal.fits\"\n", - "s2dfilename = \"det_image_seq5_MIRIMAGE_P750Lexp1_s2d.fits\"\n", - "x1dfilename = \"det_image_seq5_MIRIMAGE_P750Lexp1_x1d.fits\"\n", - "spatialprofilefilename = \"det_image_seq1_MIRIMAGE_P750Lexp1_s2d.fits\"\n", - "mainurl = \"https://data.science.stsci.edu/redirect/JWST/jwst-data_analysis_tools/MIRI_LRS_notebook/\"\n", - "\n", - "calfile_dld = download_file(mainurl + calfilename)\n", - "s2dfile_dld = download_file(mainurl + s2dfilename)\n", - "x1dfile_dld = download_file(mainurl + x1dfilename)\n", - "spatialprofilefile_dld = download_file(mainurl + spatialprofilefilename)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# rename files so that they have the right extensions\n", - "# required for the jwst datamodels to work\n", - "import os\n", - "calfile = calfile_dld + '_cal.fits'\n", - "os.rename(calfile_dld, calfile)\n", - "s2dfile = s2dfile_dld + '_s2d.fits'\n", - "os.rename(s2dfile_dld, s2dfile)\n", - "x1dfile = x1dfile_dld + '_x1d.fits'\n", - "os.rename(x1dfile_dld, x1dfile)\n", - "spatialprofilefile = spatialprofilefile_dld + '_s2d.fits'\n", - "os.rename(spatialprofilefile_dld, spatialprofilefile)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## File information\n", - "\n", - "The data used is a simulation of a LRS slit observation for a blackbody with a similar flux density to the star BD+60d1753, a flux calibration star. This simulation was created with MIRISim.\n", - "The simulated exposure was reduced using the JWST pipeline (v0.16.1) through the Detector1 and Spec2 stages.\n", - "\n", - "The cal file is one of the Spec2 products and is the calibration full frame image. It contains:\n", - "\n", - "1. (Primary): This HDU contains meta-data related to the observation and data reduction.\n", - "2. (SCI): The calibrated image. Units are MJy/sr.\n", - "3. (ERR): Uncertainty image. Units are MJy/sr.\n", - "4. (DQ): Data quality image.\n", - "5. (VAR_POISSON): Unc. component 1: Poisson uncertainty image. Units are (MJy/sr)^2.\n", - "6. (VAR_RNOISE): Unc. component 2: Read Noise uncertainty image. Units are (MJy/sr)^2.\n", - "7. (VAR_FLAT): Unc. component 3: Flat Field uncertainty image. Units are (MJy/sr)^2.\n", - "8. (ASDF_METADATA): Metadata.\n", - "\n", - "The s2d file is one of the Spec2 products and containes the calibrated rectified cutout of the LRS Slit region. It has:\n", - "\n", - "1. (Primary): This HDU contains meta-data related to the observation and data reduction.\n", - "2. (WGT): Weight.\n", - "3. (CON): ??\n", - "4. (ASDF_METADATA): Metadata." - ] - }, - { - "cell_type": "markdown", - "metadata": { - "slideshow": { - "slide_type": "slide" - } - }, - "source": [ - "## Loading data" - ] - }, - { - "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", - "cal = datamodels.open(calfile)\n", - "s2d = datamodels.open(s2dfile)" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Basic information about the image." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "print(\"cal image\")\n", - "print(cal.data.shape)\n", - "print(np.mean(cal.data))\n", - "print(np.amin(cal.data), np.amax(cal.data))\n", - "print(\"s2d image\")\n", - "print(s2d.data.shape)\n", - "print(np.mean(s2d.data))\n", - "print(np.amin(s2d.data), np.amax(s2d.data))" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Display the full 2D image" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "norm_data = simple_norm(cal.data, 'sqrt')\n", - "plt.figure(figsize=(6, 6))\n", - "plt.imshow(cal.data, norm=norm_data, origin=\"lower\")\n", - "plt.title(\"The full image from the MIRI IMAGER detector\")" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Display the LRS Slit region only (use s2d)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# transpose to make it display better\n", - "data_lrs_reg = np.transpose(s2d.data)\n", - "norm_data = simple_norm(data_lrs_reg, \"sqrt\")\n", - "plt.figure(figsize=(10, 3))\n", - "plt.imshow(data_lrs_reg, norm=norm_data, origin=\"lower\")\n", - "plt.title(\"The LRS region\")" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "JWST pipeline 1D extraction" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# for reference read in the JWST pipeline extracted spectrum\n", - "jpipe_x1d = Table.read(x1dfile, hdu=1)\n", - "print(jpipe_x1d.columns)\n", - "# plot\n", - "fig, ax = plt.subplots(figsize=(6, 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": [ - "## Boxcar Extraction\n", - "\n", - "Extract a 1D spectrum using a simple boxcar. Basically collapse the spectrum in the cross-dispersion direction over a specified number of pixels.\n", - "\n", - "Limitation: currently it is assumed there are no bad pixels." - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### Fixed width boxcar" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Define extraction parameters" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "ext_center = 30\n", - "ext_width = 8\n", - "bkg_offset = 4\n", - "bkg_width = 2" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Plot cross-disperion cut showing the extraction parameters" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "fig, ax = plt.subplots(figsize=(6, 6))\n", - "y = np.arange(data_lrs_reg.shape[0])\n", - "ax.plot(y, data_lrs_reg[:,300], 'k-')\n", - "mm = np.array([ext_center, ext_center])\n", - "mm_y = ax.get_ylim()\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", - "ax.set_title(\"Cross-dispersion Cut at Pixel=300\")" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Do the extraction" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# visualize the weight images used in the fixed boxcar extraction \n", - "wimage_fixedboxcar, wimage_fb_bkg = ap_weight_images(ext_center, ext_width, bkg_offset, \n", - " bkg_width, data_lrs_reg.shape, None)\n", - "\n", - "norm_data = simple_norm(wimage_fixedboxcar)\n", - "plt.figure(figsize=(10, 3))\n", - "plt.imshow(wimage_fixedboxcar, norm=norm_data, origin=\"lower\")\n", - "plt.title(\"Fixed boxcar weight image\")\n", - "\n", - "norm_data = simple_norm(wimage_fb_bkg)\n", - "plt.figure(figsize=(10, 3))\n", - "plt.imshow(wimage_fb_bkg, norm=norm_data, origin=\"lower\")\n", - "plt.title(\"Fixed boxcar backgound weight image\")" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# extract the spectrum using the weight images\n", - "\n", - "# without background subtraction\n", - "waves_boxcar, ext1d_boxcar, tmpval = extract_1dspec(s2d, ext_center, ext_width, None, None)\n", - "\n", - "# with background subtraction\n", - "waves_boxcar_bkgsub, ext1d_boxcar_bkgsub, tmpval = extract_1dspec(s2d, ext_center, ext_width, \n", - " bkg_offset, bkg_width)\n", - "\n", - "# plot\n", - "fig, ax = plt.subplots(figsize=(6, 6))\n", - "gpts = ext1d_boxcar_bkgsub > 0.\n", - "ax.plot(waves_boxcar[gpts], ext1d_boxcar[gpts], 'k-', label=\"boxcar\")\n", - "ax.plot(waves_boxcar_bkgsub[gpts], ext1d_boxcar_bkgsub[gpts], 'k:', label=\"boxcar (bkgsub)\")\n", - "ax.plot(jpipe_x1d['WAVELENGTH'], jpipe_x1d['FLUX'], 'k-', label=\"jpipe_x1d\")\n", - "ax.set_title(\"Fixed boxcar 1D extracted spectrum\")\n", - "ax.set_xlabel(r\"wavelength [$\\mu$m]\")\n", - "ax.set_ylabel(\"Flux Density [Jy]\")\n", - "ax.set_yscale(\"log\")\n", - "ax.legend()" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### Wavelength scaled width boxcar\n", - "\n", - "The LRS spatial profile changes as a function of wavelength as JWST is diffraction limited at these wavelengths. Nominally this means that the FWHM is changing linearly with wavelength. Scaling the width of the extraction aperture with wavelength accounts for the changing diffraction limit with wavelength to first order." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# visualize the weight images used in the scaled boxcar extraction \n", - "wimage_scaledboxcar, wimage_sb_bkg = ap_weight_images(ext_center, ext_width, bkg_offset, bkg_width, \n", - " data_lrs_reg.shape, waves_boxcar, wavescale=10.0)\n", - "\n", - "norm_data = simple_norm(wimage_scaledboxcar)\n", - "plt.figure(figsize=(10, 3))\n", - "plt.imshow(wimage_scaledboxcar, norm=norm_data, origin=\"lower\")\n", - "plt.title(\"Scaled boxcar weight image\")\n", - "\n", - "norm_data = simple_norm(wimage_sb_bkg)\n", - "plt.figure(figsize=(10, 3))\n", - "plt.imshow(wimage_sb_bkg, norm=norm_data, origin=\"lower\")\n", - "plt.title(\"Scaled boxcar backgound weight image\")" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# extract the spectrum using the weight image\n", - "\n", - "# with background subtraction\n", - "waves_sboxcar_bkgsub, ext1d_sboxcar_bkgsub, sboxcar_bkgsub_image = extract_1dspec(s2d, ext_center, \n", - " ext_width, bkg_offset, \n", - " bkg_width, wavescale=10)\n", - "\n", - "# plot\n", - "fig, ax = plt.subplots(figsize=(6, 6))\n", - "gpts = ext1d_boxcar_bkgsub > 0.\n", - "ax.plot(waves_boxcar_bkgsub[gpts], ext1d_boxcar_bkgsub[gpts], 'k:', label=\"fixed boxcar (bkgsub)\")\n", - "gpts = ext1d_sboxcar_bkgsub > 0.\n", - "ax.plot(waves_sboxcar_bkgsub[gpts], ext1d_sboxcar_bkgsub[gpts], 'k-', label=\"scaled boxcar (bkgsub)\")\n", - "ax.set_title(\"Scaled boxcar 1D extracted spectrum\")\n", - "ax.set_xlabel(\"wavelength [$\\mu$m]\")\n", - "ax.set_ylabel(\"Flux Density [Jy]\")\n", - "ax.set_yscale(\"log\")\n", - "ax.set_ylim(1e-3, 1e-1)\n", - "ax.legend()" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Note that the impact of the scaled boxcar is largest at shorter wavelengths. This is the result of using the same aperature at 10 microns for both the boxcar and scaled boxcar." - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "\n", - "## PSF based Extraction\n", - "\n", - "While to first order the PSF FHWM changes linearly with wavelength, this is an approximation. It is better to use the measured spatial profile as a function of wavelength to extract the spectrum. This tracks the actual variation with wavelength and optimizes the extraction to the higher S/N measurements. In general, PSF based extractions show the most improvements over boxcar extractions at lower the S/N.\n", - "\n", - "There are two PSF based extraction methods.\n", - "\n", - "1. PSF weighted: the spatial profile at each wavelength is used to weight the extraction.\n", - "2. PSF fitting: the spatial profile is fit at each wavelength with the scale parameter versus wavelength giving the spectrum.\n", - "\n", - "Only the PSF weighted technique is currently part of this notebook.\n", - "\n", - "Note 1: calibration reference file for the specific LRS slit position should be used.\n", - "\n", - "Note 2: Small shifts in the centering of the source in the slit should be investigated to see if they impact the PSF based extractions.\n", - "\n", - "Limitation: currently it is assumed there are no bad pixels." - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### PSF weighted extaction" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "#### Generate the PSF profile as a function of wavelength\n", - "For MIRI LRS slit observations, observations are made at two nod position in the slit after target acquisition. This means that the location of the sources in the slit is very well known. Hence, spatial profile (PSF) as a function of wavelength for the two nod positions is straightforward to measure using observations of a bright source.\n", - "\n", - "The next few steps generate the needed information for the nod position for which we are extracting spectra based on a simulation of a bright source at the same nod position." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# lrs spatial profile (PSF) as a function of wavelength\n", - "# currently, this is just a \"high\" S/N observation of a flat spectrum source at the same slit position\n", - "psf = datamodels.open(spatialprofilefile)\n", - "# transpose to make it display better\n", - "lrspsf = np.transpose(psf.data)\n", - "norm_data = simple_norm(lrspsf, \"sqrt\")\n", - "plt.figure(figsize=(10, 3))\n", - "plt.imshow(lrspsf, norm=norm_data, origin=\"lower\")\n", - "plt.title(\"The LRS Spatial Profile (PSF) Observation\")" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# Mock a LRS spectral profile reference file\n", - "# Sum along the spatial direction and normalize to 1\n", - "# assume there is no background (none was included in the MIRISim for the flat spectrum source observation)\n", - "# ignore regions far from the source using a scaled boxcar weight image\n", - "# the aperture (psf_width) used in the scaled boxcar weight image could be varied\n", - "psf_width = 12.0\n", - "(wimage_scaledboxcar, tmpvar) = ap_weight_images(ext_center, psf_width, bkg_offset, bkg_width, data_lrs_reg.shape, waves_boxcar, wavescale=10.0)\n", - "\n", - "psf_weightimage = lrspsf*wimage_scaledboxcar\n", - "\n", - "# generate a 2D image of the column sums for division\n", - "max_psf = np.max(psf_weightimage, axis=0)\n", - "div_image = np.tile(max_psf, (psf_weightimage.shape[0], 1))\n", - "div_image[div_image == 0.0] = 1.0 # avoid divide by zero issues\n", - "\n", - "# normalize \n", - "psf_weightimage /= div_image\n", - "\n", - "# display\n", - "norm_data = simple_norm(psf_weightimage, \"sqrt\")\n", - "plt.figure(figsize=(10, 3))\n", - "plt.imshow(psf_weightimage, norm=norm_data, origin=\"lower\")\n", - "plt.title(\"The LRS Spatial Profile Reference Image (Normalized)\")" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "fig, ax = plt.subplots(figsize=(6, 6))\n", - "y = np.arange(psf_weightimage.shape[0])\n", - "ax.plot(y, psf_weightimage[:,150], label=\"pixel=150\")\n", - "ax.plot(y, psf_weightimage[:,225], label=\"pixel=225\")\n", - "ax.plot(y, psf_weightimage[:,300], label=\"pixel=300\")\n", - "ax.plot(y, psf_weightimage[:,370], label=\"pixel=370\")\n", - "ax.set_title(\"Cross-dispersion Cuts\")\n", - "ax.set_xlim(ext_center-psf_width, ext_center+psf_width)\n", - "ax.legend()" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Note that the spatial profile becomes narrower as the pixel values increases as this corresponds to the wavelength decreasing." - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "#### Extract spectrum using wavelength dependent PSF profiles" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# use the normalized PSF weight image to extract the specrum\n", - "# use the background subtracted image from the scaled boxcar extraction\n", - "ext1d_psfweight = np.sum(sboxcar_bkgsub_image * psf_weightimage, axis=0)\n", - "ext1d_psfweight *= 1e6 * s2d.meta.photometry.pixelarea_steradians\n", - "\n", - "# plot\n", - "fig, ax = plt.subplots(figsize=(6, 6))\n", - "gpts = ext1d_psfweight > 0.\n", - "ax.plot(waves_boxcar_bkgsub[gpts], ext1d_psfweight[gpts], 'k-', label=\"psf weighted (bkgsub)\")\n", - "gpts = ext1d_sboxcar_bkgsub > 0.\n", - "ax.plot(waves_sboxcar_bkgsub[gpts], ext1d_sboxcar_bkgsub[gpts], 'k:', label=\"scaled boxcar (bkgsub)\")\n", - "ax.set_title(\"PSF weigthed extracted spectrum\")\n", - "ax.set_xlabel(\"wavelength [$\\mu$m]\")\n", - "ax.set_ylabel(\"Flux Density [Jy]\")\n", - "ax.set_yscale(\"log\")\n", - "ax.set_ylim(1e-3, 1e-1)\n", - "ax.legend()" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "Note that the psf weighted extraction has visabily higher S/N, especially at the longer wavelengths where the S/N is lowest overall." - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "#### Plotting in Rayleigh-Jeans units\n", - "\n", - "For sources that have stellar continuum, it can be useful to plot MIR spectra in Rayleigh-Jeans units. This just means removing the spectral shape expected for a blackbody with a peak at much shorter wavelengths than the MIR. This is easily done by multiplying the spectrum by lambda^4 or nu^2.\n", - "\n", - "An example of this is given below." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# Rayleigh-Jeans plot\n", - "fig, ax = plt.subplots(figsize=(6, 6))\n", - "gpts = ext1d_psfweight > 0.\n", - "ax.plot(waves_boxcar_bkgsub[gpts], (waves_boxcar_bkgsub[gpts]**4)*ext1d_psfweight[gpts], 'k-', label=\"psf weighted (bkgsub)\")\n", - "gpts = ext1d_sboxcar_bkgsub > 0.\n", - "ax.plot(waves_sboxcar_bkgsub[gpts], (waves_sboxcar_bkgsub[gpts]**4)*ext1d_sboxcar_bkgsub[gpts], 'k:', label=\"scaled boxcar (bkgsub)\")\n", - "gpts = ext1d_boxcar_bkgsub > 0.\n", - "ax.plot(waves_boxcar_bkgsub[gpts], (waves_boxcar_bkgsub[gpts]**4)*ext1d_boxcar_bkgsub[gpts], 'k--', label=\"fixed boxcar (bkgsub)\")\n", - "ax.set_title(\"Rayleigh-Jeans plot for all extractions\")\n", - "ax.set_xlabel(\"wavelength [$\\mu$m]\")\n", - "ax.set_ylabel(\"Rayleigh-Jeans Flux Density [$\\mu$m$^4$ Jy]\")\n", - "ax.set_yscale(\"log\")\n", - "ax.set_ylim(10, 100)\n", - "ax.legend()" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## Additional Resources\n", - "\n", - "- [MIRI LRS](https://jwst-docs.stsci.edu/mid-infrared-instrument/miri-observing-modes/miri-low-resolution-spectroscopy)\n", - "- [MIRISim](http://www.stsci.edu/jwst/science-planning/proposal-planning-toolbox/mirisim)\n", - "- [JWST pipeline](https://jwst-docs.stsci.edu/jwst-data-reduction-pipeline)\n", - "- PSF weighted extraction [Horne 1986, PASP, 98, 609](https://ui.adsabs.harvard.edu/abs/1986PASP...98..609H/abstract)." - ] - }, - { - "cell_type": "markdown", - "metadata": { - "slideshow": { - "slide_type": "slide" - } - }, - "source": [ - "## About this notebook\n", - "\n", - "**Author:** Karl Gordon, JWST\n", - "**Updated On:** 2020-07-07" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "***" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "[Top of Page](#top)\n", - "\"Space " - ] - } - ], - "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.10" - } - }, - "nbformat": 4, - "nbformat_minor": 2 -}