diff --git a/notebooks/ifu_optimal/ifu_optimal.ipynb b/notebooks/ifu_optimal/ifu_optimal.ipynb index 8b60f9c26..e27edc23b 100644 --- a/notebooks/ifu_optimal/ifu_optimal.ipynb +++ b/notebooks/ifu_optimal/ifu_optimal.ipynb @@ -4,42 +4,73 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "# IFU Optimal Spectral Extraction" + "# NIRSpec IFU Optimal Point Source Extraction" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "**Use case:** optimal spectral extraction; method by [Horne (1986)](https://ui.adsabs.harvard.edu/abs/1986PASP...98..609H/abstract).
\n", - "**Data:** JWST simulated NIRSpec IFU data; point sources.
\n", - "**Tools:** jwst, webbpsf, matplotlib, scipy, custom functions.
\n", - "**Cross-intrument:** any spectrograph.
\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).
" + "This notebook illustrates various extraction methods for a point source in JWST NIRSpec IFU data, utilizing the [Q3D](https://q3d.github.io/) (PID 1335) observation of quasar SDSS J165202.64+172852.3. The extraction techniques include subset extraction with Cubeviz, simple sum over spaxels, cylindrical aperture, conical aperture photometry, and optimal point source extraction using a WebbPSF model PSF. " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "## Imports " + "**Use case:** optimal spectral extraction; method by [Horne (1986)](https://ui.adsabs.harvard.edu/abs/1986PASP...98..609H/abstract).\n", + "\n", + "**Data:** JWST NIRSpec IFU data; point sources.\n", + "\n", + "**Tools:** jwst, webbpsf, matplotlib, scipy, custom functions.\n", + "\n", + "**Cross-intrument:** any spectrograph. \n", + "\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)." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Table of Contents\n", + "1. [Imports](#imports)\n", + "2. [Read in NIRSpec IFU Cube](#read)\n", + "3. [Visualize Science Data with Cubeviz](#visualize)\n", + "4. [Export Source and Good Data Regions from Cubeviz](#source)\n", + "5. [Extract Subset Spectrum and Background from Cubeviz Spectrum Viewer](#viewer)\n", + "6. [Extract Spectrum by Sum over Spaxels](#sum)\n", + "7. [Extract Spectrum in Constant Radius Circular Aperture (Cylinder)](#cylinder)\n", + "8. [Extract Spectrum in Linearly Expanding Circular Aperture (Cone)](#cone)\n", + "9. [Plot and Compare Non-optimal Spectral Extractions](#plot)\n", + "10. [WebbPSF Model for Optimal Extraction](#webbpsf)\n", + "11. [Align Model PSF Cube with Science Data](#align)\n", + "12. [Optimal Extraction Using WebbPSF Model](#optimal)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## 1. Imports " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "* _time_ for adding a delay in Cubeviz file upload\n", "* _numpy_ for array math\n", - "* _scipy_ for gaussian smoothing\n", + "* _scipy_ for ndimage shift\n", "* _specutils_ for Spectrum1D data model\n", - "* _jdaviz_ : Cubeviz data visualization tool\n", + "* _jdaviz_ : for data visualization\n", "* _photutils_ to define circular apertures\n", "* _astropy.io_ for reading and writing FITS cubes and images\n", "* _astropy.wcs, units, coordinates_ for defining and reading WCS\n", "* _astropy.stats_ for sigma_clipping\n", "* _astropy.utils_ for downloading files from URLs\n", - "* _matplotlib_ for plotting spectra and images" + "* _matplotlib_ for plotting spectra and images\n", + "* _os_ for file management\n", + "* _astroquery.mast_ to download the data" ] }, { @@ -48,16 +79,18 @@ "metadata": {}, "outputs": [], "source": [ - "import time\n", "import numpy as np\n", "import scipy\n", "from specutils import Spectrum1D\n", - "from jdaviz import Cubeviz\n", - "from photutils import CircularAperture, aperture_photometry \n", + "import jdaviz\n", + "from jdaviz import Cubeviz, Specviz\n", + "print(\"jdaviz Version={}\".format(jdaviz.__version__))\n", + "from photutils.aperture import CircularAperture, aperture_photometry \n", "from astropy.io import fits\n", "from astropy import wcs\n", - "from astropy.stats import sigma_clip\n", - "from astropy.utils.data import download_file" + "from astropy import units as u\n", + "import os\n", + "from astroquery.mast import Observations" ] }, { @@ -75,21 +108,32 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "## Introduction\n", + "## 2. Read in NIRSpec IFU Cube \n", "\n", - " This notebook illustrates various extraction methods for a point source in JWST NIRSpec IFU data. First we\n", - "demonstrate a number of regular extraction techniques, including subset extraction with Cubeviz, simple sum over spaxels, cylindrical aperture, and conical aperture photometry. Then we compare optimal extraction using a WebbPSF model PSF to optimal extraction using a reference star PSF. \n" + "The NIRSpec IFU observation of quasar SDSS J1652+1728 (redshift z=1.9) was taken using the G235H grating with F170LP filter, covering 1.66-3.17 microns at a spectral resolution of R~2700. The IFU spaxels are 0.1\" on a side. \n", + "The level-3 pipeline_processed datacube (s3d.fits, which combines all dithered exposures) is retrieved from MAST \n", + "in the next notebook cell below." ] }, { - "cell_type": "markdown", + "cell_type": "code", + "execution_count": null, "metadata": {}, + "outputs": [], "source": [ - "## Read in Simulated NIRSpec IFU Cube\n", + "# Download the data file\n", + "uri = \"mast:jwst/product/jw01335-o008_t007_nirspec_g235h-f170lp_s3d.fits\"\n", + "result = Observations.download_file(\n", + " uri, base_url=\"https://mast.stsci.edu/api/v0.1/Download/file\"\n", + ")\n", + "if result[0] == \"ERROR\":\n", + " raise RuntimeError(\"Error retrieving file: \" + result[1])\n", "\n", - "A faint (quasar) point source was simulated using the NIRSpec Instrument Performance Simulator (IPS), then run through the JWST Spec2 pipeline. We will use this for our science dataset.\n", + "# Construct the local filepath\n", + "filename = os.path.join(os.path.abspath(\".\"), uri.rsplit(\"/\", 1)[-1])\n", "\n", - "We read in the data both with fits.open and Spectrum1D.read, since the cube handling (slicing) we need to do is not implemented in specutils yet." + "# Optionally Replace MAST data with custom reprocessed data in the current directory\n", + "# filename=\"./jw01335-o008_t007_nirspec_g235h-f170lp_s3d.fits\"" ] }, { @@ -98,70 +142,30 @@ "metadata": {}, "outputs": [], "source": [ - "# NIRSpec IFU science data cube\n", - "BoxPath = \"https://data.science.stsci.edu/redirect/JWST/jwst-data_analysis_tools/IFU_optimal_extraction/\"\n", - "filename = BoxPath + \"NRS00001-faintQSO-F100LP-G140H-01_1_491_SE_2020-08-25T12h15m00_s3d.fits\"\n", - "\n", "# Open and inspect the file and WCS\n", - "# Load with astropy.fits.open\n", "with fits.open(filename, memmap=False) as hdulist:\n", - " sci = hdulist['SCI'].data\n", - " err = hdulist['ERR'].data\n", + " sci = hdulist[\"SCI\"].data\n", + " err = hdulist[\"ERR\"].data\n", " w = wcs.WCS(hdulist[1].header)\n", " hdr = hdulist[1].header\n", " hdulist.info()\n", " print(w)\n", - " \n", - "# Load with Spectrum1D \n", - "spec1d = Spectrum1D.read(filename)\n", - "\n", - "# Wavelengths\n", - "wavelength = np.array(spec1d.spectral_axis.value)\n", - "print(wavelength)\n", "\n", - "# Sum over spaxels\n", - "fnu_sum = np.sum(spec1d.flux, axis=(0, 1))\n", + "# Window the wavelength range to focus on Hbeta-[OIII]\n", + "spec1d = Spectrum1D.read(filename)\n", + "slice_range = range(500, 1100, 1)\n", + "wavelength = np.array(spec1d.spectral_axis.value)[slice_range[0]: slice_range[-1] + 1]\n", "\n", "# List of cube slices for aperture photometry\n", - "data = []\n", - "var = []\n", - "spec1d_len = len(spec1d.spectral_axis.value)\n", - "for idx in range(spec1d_len): \n", - " data.append(sci[idx, :, :])\n", - " var.append(err[idx, :, :]) # variance = err, not variance = err**2. Squaring the err gives noisy results. \n", + "sci_data = []\n", + "sci_var = []\n", + "for idx in slice_range:\n", + " sci_data.append(sci[idx, :, :])\n", + " sci_var.append(err[idx, :, :])\n", "\n", - "# Window data and variance (and replace NaNs)\n", - "# The existing JWST pipeline window is overgenerous (39x33 instead of the nominal 30x30 pixels)\n", - "data_win = np.nan_to_num(np.array(data)[:, 5:-4, 3:])\n", - "data_var = np.nan_to_num(np.array(var)[:, 5:-4, 3:]) " - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "*Developer Note:* Can we fix or suppress this AsdfWarning?" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "*Developer Note:* Is there a way to read in with only Spectrum1D to perform all of our cube operations, not using fits.open()?" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## Visualize Science Data with Cubeviz" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "*Developer Note:* Cubeviz incompatible with jupyter_client 6.1.6. Use jupyter_client 5.3.5 instead." + "data = np.nan_to_num(np.array(sci_data))\n", + "var = np.array(sci_var)\n", + "print(\"\\nTrimmed data shape:\", data.shape)" ] }, { @@ -171,88 +175,48 @@ "outputs": [], "source": [ "cubeviz = Cubeviz()\n", - "cubeviz.app" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "### UI Instructions:\n", - "* Load science datacube into Cubeviz using the next code cell below\n", - "* Go to the Hammer-and-Screwdriver icon: Gear icon: Layer in the leftmost image viewer \n", - "* In that tab, change the Linear stretch to 90 percentile to see the faint QSO target at (x,y) ~ (17, 21)\n", - "* Scrubbing through the cube also helps to locate the source\n", - "* Select a circular subset region centered on the source. \n", - "* Note that the region is pixelated and doesn't include fractional pixels\n", - "* Change the collapse method to \"Sum\" in spectrum viewer: Gear icon : Viewer \n", - "* --This \"Sum\" method yields our subset extraction\n", - "* Change the vertical zoom to see the spectral features in the Subset spectrum" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "*Developer Notes*: \n", + "cubeviz.load_data(filename)\n", + "cubeviz.show()\n", "\n", - "(1) Image viewer contrast settings change when you click on the side bar to expand/contract jupyter scroll window\n", + "# Set spectrum display limits\n", + "cubeviz.specviz.x_limits(1.65 * u.um, 2.4 * u.um)\n", + "cubeviz.specviz.y_limits(0.0, 5.0e3)\n", "\n", - "(2) Spectrum viewer: Viewer cube collapse method should default to Sum (not Maximum)\n", - "\n", - "(3) Spectrum viewer y scale returns to autoscale when the region is moved, and y-zoom has to be adjusted again\n", - "\n", - "(4) Region selection appears away from cursor after opening hammer-and-screwdriver to change cube viewer contrast" + "# Select slice to visualize\n", + "cubeviz.select_slice(714)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "## Load Cube into Cubeviz" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# Data from local directory\n", - "# cubeviz.app.load_data(filename)\n", - "\n", - "# Data from url:\n", - "url = filename\n", - "df = download_file(url)\n", - "time.sleep(2) # Sleep to avoid glue-jupyter timing issue\n", - "cubeviz.app.load_data(df)" + "## 3. Visualize Science Data with Cubeviz " ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "*Developer Note:* Spectral cube does not yet recognize JWST NIRSpec IFU datacubes, giving the above warning\n", - "for each FITS extension." + "\"Cubeviz," ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "## Export Region from Cubeviz\n", - "Export the region defined by the user in Cubeviz as an astropy CirclePixel Region, which has units of pixels." + "### UI Instructions:\n", + "* Scrub through the cube to a line-free region using the spectrum-viewer slice tool.\n", + "* In the flux-viewer, select one circular subset region centered on the quasar, and a square region to delimit the good area for spectral and background extraction\n", + "* Note that the regions are pixelated and don't include fractional pixels\n", + "* The default collapse method is \"Sum\" in the spectrum viewer (see Plot Options:Line). \"Median\" may also be useful for visualization but will not give an accurate measurement of the total flux.)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "*Developer Note:* cubeviz.app.get_subsets_from_viewer method doesn't work if there are more than 2 datasets selected in the spectrum viewer:\n", - "\n", - "#region1 = cubeviz.app.get_subsets_from_viewer('spectrum-viewer')\n", - "\n", - "#print(region1['Subset1'])\n" + "## 4. Export Source and Good Data Regions from Cubeviz \n", + "Export the region defined by the user in Cubeviz as astropy PixelRegions" ] }, { @@ -262,20 +226,63 @@ "outputs": [], "source": [ "cubeviz_data = cubeviz.app.data_collection[0]\n", + "\n", + "print(\"\\nSource Region\")\n", "try:\n", - " region1 = cubeviz_data.get_selection_definition(format='astropy-regions')\n", + " region1 = cubeviz_data.get_selection_definition(\n", + " \"Subset 1\", format=\"astropy-regions\"\n", + " )\n", " print(region1)\n", + " center_xy = [region1.center.x, region1.center.y]\n", + " r_pix = region1.radius\n", " region1_exists = True\n", "except Exception:\n", - " print(\"There are no regions selected in the cube viewer.\")\n", - " region1_exists = False" + " print(\"There is no Subset 1 selected in the cube viewer.\")\n", + " center_xy = [17.1, 20.0]\n", + " r_pix = 5.92\n", + " print(\"Using default pixel center and radius:\")\n", + " print(\"Center pixel:\", center_xy)\n", + " print(\"Radius (pixels):\", r_pix)\n", + " region1_exists = False\n", + " \n", + "print(\"\\nGood Data Region\")\n", + "try:\n", + " region2 = cubeviz_data.get_selection_definition(\n", + " \"Subset 2\", format=\"astropy-regions\"\n", + " )\n", + " print(region2)\n", + " region2_exists = True\n", + " data_xrange = [\n", + " round(region2.center.x - region2.width / 2),\n", + " round(region2.center.x + region2.width / 2),\n", + " ]\n", + " data_yrange = [\n", + " round(region2.center.y - region2.height / 2),\n", + " round(region2.center.y + region2.height / 2),\n", + " ]\n", + " print(\"Good data (xmin,xmax), (ymin,ymax):\", data_xrange, data_yrange)\n", + " good_data = np.nan_to_num(\n", + " data[:, data_yrange[0]: data_yrange[1], data_xrange[0]: data_xrange[1]]\n", + " )\n", + " good_var = var[:, data_yrange[0]: data_yrange[1], data_xrange[0]: data_xrange[1]]\n", + "\n", + "except Exception:\n", + " print(\"There is no Subset 2 selected in the cube viewer.\")\n", + " region1_exists = False\n", + " data_xrange = [7, 36]\n", + " data_yrange = [6, 33]\n", + " good_data = np.nan_to_num(\n", + " data[:, data_yrange[0]: data_yrange[1], data_xrange[0]: data_xrange[1]]\n", + " )\n", + " good_var = var[:, data_yrange[0]: data_yrange[1], data_xrange[0]: data_xrange[1]]\n", + " print(\"Using default good data (xmin,xmax), (ymin,ymax):\", data_xrange, data_yrange)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "## Extract Subset Spectrum in Cubeviz Spectrum Viewer\n", + "## 5. Extract Subset Spectrum and Background from Cubeviz Spectrum Viewer \n", "Retrieve the collapsed spectrum (Subset1) of the user-defined region from the Spectrum Viewer as a Spectrum1D object." ] }, @@ -285,25 +292,38 @@ "metadata": {}, "outputs": [], "source": [ + "subsets = cubeviz.specviz.get_spectra()\n", + "print(subsets.keys())\n", + "\n", + "print(\"\\nSource\")\n", "try:\n", - " spectrum_subset1 = cubeviz.app.get_data_from_viewer('spectrum-viewer')['Subset 1']\n", + " spectrum_subset1 = subsets[[i for i in subsets.keys() if \"Subset 1\" in i][0]]\n", " print(spectrum_subset1)\n", "except Exception:\n", - " print(\"There are no subsets selected in the spectrum viewer.\")" + " print(\"There is no Subset 1 selected in the spectrum viewer.\")\n", + "\n", + "print(\"\\nBackground\")\n", + "try:\n", + " spectrum_subset2 = spectrum_subset2 = subsets[\n", + " [i for i in subsets.keys() if \"Subset 2\" in i][0]\n", + " ]\n", + " print(spectrum_subset2)\n", + "except Exception:\n", + " print(\"There is no Subset 2 selected in the spectrum viewer.\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "*Developer Note:* Can we suppress or fix this glue/core warning?" + "*Developer Note:* The units of the Cubeviz \"Sum\" collapse method (and all other spectra below) need to be multiplied by the pixel area in sr to yield flux units (MJy) instead of surface brightness units (MJy/sr)." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "## Extract Spectrum by Sum Over Spaxels\n", + "## 6. Extract Spectrum by Sum Over Spaxels \n", "\n", "Perform a simple numpy sum over all spaxels in the cube as a rudimentary extraction method. Also sum over wavelength to collapse the cube." ] @@ -315,17 +335,26 @@ "outputs": [], "source": [ "# Sum over wavelength\n", - "cube_sum = np.sum(data_win, axis=0)\n", + "# Clip data for display purposes\n", + "clip_level = 4e4\n", + "data_clipped = np.clip(good_data, 0, clip_level)\n", + "cube_sum = np.sum(data_clipped, axis=0)\n", + "\n", + "# Extraction via sum over spaxels\n", + "fnu_sum = np.sum(good_data, axis=(1, 2))\n", + "fnu_sum_clipped = np.clip(fnu_sum, 0, clip_level)\n", + "flux_spaxsum = np.array(fnu_sum) * u.MJy / u.sr\n", + "spec1d_spaxsum = Spectrum1D(spectral_axis=wavelength * u.um, flux=flux_spaxsum)\n", "\n", "# Plots\n", - "f, (ax1, ax2) = plt.subplots(1, 2, figsize=(10, 5)) \n", - "ax1.plot(wavelength, fnu_sum) \n", - "ax1.set_xlim(0.95, 1.5)\n", - "ax1.set_title(\"Spaxel sums\")\n", - "ax1.set_xlabel(\"Wavelength (um)\") \n", - "ax1.set_ylabel(\"SCA491 Flux Density\")\n", + "f, (ax1, ax2) = plt.subplots(1, 2, figsize=(10, 4))\n", "\n", - "ax2.imshow(cube_sum, norm=LogNorm())\n", + "ax1.plot(wavelength, fnu_sum)\n", + "ax1.set_title(\"Spaxel sums\")\n", + "ax1.set_xlabel(\"Wavelength (um)\")\n", + "ax1.set_ylabel(\"Flux (MJy/sr)\")\n", + "ax1.set_ylim(0, 5e3)\n", + "ax2.imshow(cube_sum, norm=LogNorm(vmin=100, vmax=clip_level), origin=\"lower\")\n", "ax2.set_title(\"Slice sums\")\n", "\n", "plt.show()" @@ -335,7 +364,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "## Extract Spectrum in Constant Radius Circular Aperture (Cylinder)\n", + "## 7. Extract Spectrum in Constant Radius Circular Aperture (Cylinder) \n", "This method is appropriate for an extended source." ] }, @@ -345,24 +374,17 @@ "metadata": {}, "outputs": [], "source": [ - "# IFU pixel scale\n", - "pixelscale = 0.1 # arcsec/pixel\n", - "\n", "# CircularAperture uses xy pixels\n", - "center_xy = [17.1, 20.]\n", - "r_pix = 5.92\n", - "if region1_exists:\n", - " center_xy = [region1.center.x, region1.center.y] \n", - " r_pix = region1.radius\n", - "\n", "aperture = CircularAperture(center_xy, r=r_pix)\n", "print(aperture)\n", "\n", "cylinder_sum = []\n", "for slice2d in data:\n", - " phot_table = aperture_photometry(slice2d, aperture, wcs=w.celestial, method='exact')\n", " phot_table = aperture_photometry(slice2d, aperture)\n", - " cylinder_sum.append(phot_table['aperture_sum'][0])" + " cylinder_sum.append(phot_table[\"aperture_sum\"][0])\n", + "\n", + "flux_cylinder = np.array(cylinder_sum) * u.MJy / u.sr\n", + "spec1d_cylinder = Spectrum1D(spectral_axis=wavelength * u.um, flux=flux_cylinder)" ] }, { @@ -376,7 +398,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "## Extract Spectrum in Linearly Expanding Circular Aperture (Cone)\n", + "## 8. Extract Spectrum in Linearly Expanding Circular Aperture (Cone) \n", "This method is appropriate for a point source PSF with width proportional to wavelength" ] }, @@ -388,23 +410,29 @@ "source": [ "# Reference wavelength for expanding aperture\n", "lambda0 = wavelength[0]\n", - "print('Reference wavelength:', lambda0)\n", + "print(\"Reference wavelength:\", lambda0)\n", "\n", "cone_sum = []\n", "idx = -1\n", "for (slice2d, wave) in zip(data, wavelength):\n", " idx = idx + 1\n", " r_cone = r_pix * wave / lambda0\n", + "\n", " aperture_cone = CircularAperture(center_xy, r=r_cone)\n", - " phot_table = aperture_photometry(slice2d, aperture_cone, wcs=w.celestial, method='exact')\n", - " cone_sum.append(phot_table['aperture_sum'][0])" + " phot_table = aperture_photometry(\n", + " slice2d, aperture_cone, wcs=w.celestial, method=\"exact\"\n", + " )\n", + " cone_sum.append(phot_table[\"aperture_sum\"][0])\n", + "\n", + "flux_cone = np.array(cone_sum) * u.MJy / u.sr\n", + "spec1d_cone = Spectrum1D(spectral_axis=wavelength * u.um, flux=flux_cone)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "## Plot and Compare Non-optimal Spectral Extractions\n", + "## 9. Plot and Compare Non-optimal Spectral Extractions \n", "Compare spectra extracted in cylinder, cone, Cubeviz subset." ] }, @@ -414,21 +442,27 @@ "metadata": {}, "outputs": [], "source": [ - "f, (ax1) = plt.subplots(1, 1, figsize=(15, 5)) \n", + "f, (ax1) = plt.subplots(1, 1, figsize=(15, 5))\n", "\n", - "ax1.set_title(\"Non-optimal spectral extractions\")\n", - "ax1.set_xlabel(\"Observed Wavelength (microns)\") \n", - "ax1.set_ylabel(\"Flux Density\")\n", - "ax1.set_xlim(0.95, 1.5)\n", - "ax1.set_ylim(0, 0.6)\n", - "ax1.plot(wavelength, np.array(cylinder_sum), label=\"Cylinder\", c='b')\n", - "ax1.plot(wavelength, np.array(cone_sum), label=\"Cone\", c='darkorange', alpha=0.5)\n", + "# ax1.plot(wavelength, flux_spaxsum.value, label=\"All spaxels\", c='k')\n", + "ax1.plot(wavelength, flux_cylinder.value, label=\"Cylinder\", c=\"b\")\n", + "ax1.plot(wavelength, flux_cone.value, label=\"Cone\", c=\"darkorange\", alpha=0.5)\n", "try:\n", - " ax1.plot(wavelength, spectrum_subset1.flux.value, c='r', label=\"Subset1\", alpha=0.4)\n", + " ax1.plot(\n", + " wavelength,\n", + " spectrum_subset1.flux.value[slice_range[0]: slice_range[-1] + 1],\n", + " c=\"r\",\n", + " label=\"Subset1\",\n", + " alpha=0.4,\n", + " )\n", "except Exception:\n", " print(\"There is no Cubeviz Subset1 spectrum to plot.\")\n", - "ax1.legend()\n", "\n", + "ax1.set_title(\"Non-optimal spectral extractions\")\n", + "ax1.set_xlabel(\"Observed Wavelength (microns)\")\n", + "ax1.set_ylabel(\"Flux Density\")\n", + "ax1.set_ylim(0, 5.0e3)\n", + "ax1.legend()\n", "plt.show()" ] }, @@ -436,8 +470,8 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "Comparison of the (non-optimal) cylindrical, conical, and Cubeviz subset spectral extractions. \n", - "The conical extraction captures slightly more flux but is noisier than the other spectra at long wavelengths.\n", + "The non-optimal cylindrical, conical, and CubeViz subset spectral extractions are quite similar. \n", + "The conical extraction captures more flux at long wavelengths.\n", "Red-shifted Broad H-beta and narrow [O III] lines are visible in the quasar spectra. " ] }, @@ -445,7 +479,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "## WebbPSF Model PSF for Optimal Extraction\n", + "## 10. WebbPSF Model PSF for Optimal Extraction \n", "Generate PSF model cube using WebbPSF for NIRSpec IFU, or read in precomputed PSF model cube." ] }, @@ -453,7 +487,9 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "Caution! The WebbPSF model takes about 10 hr to run. Uncomment the following cell to do so. Otherwise, read in the precomputed WebbPSF model, which covers the full F100LP/G140H wavelength range (blue and red). For other filter/grating combinations, uncomment and run the cell below using the wavelengths from the science data set." + "WebbPSF installation instructions can be found in [ReadTheDocs](https://webbpsf.readthedocs.io/en/latest/).\n", + "\n", + "Caution! The WebbPSF model takes about 10 hr to run. Uncomment the following cell to do so. Otherwise, read in the precomputed WebbPSF model, which covers the full wavelength range of the present NIRSpec G235H dataset. For other filter/grating combinations, uncomment and run the cell below using the wavelengths from the science data set." ] }, { @@ -462,18 +498,17 @@ "metadata": {}, "outputs": [], "source": [ - "'''\n", - "#WebbPSF imports\n", - "%pylab inline\n", - "import webbpsf\n", + "# #WebbPSF imports\n", + "# %pylab inline\n", + "# import webbpsf\n", + "#\n", + "# #WebbPSF commands used to create PSF model cube\n", + "# ns = webbpsf.NIRSpec()\n", + "# ns.image_mask = \"IFU\" # Sets to 3x3 arcsec square mask\n", "\n", - "#WebbPSF commands used to create PSF model cube\n", - "ns.image_mask = \"IFU\" # Sets to 3x3 arcsec square mask\n", - "ns = webbpsf.NIRSpec()\n", - "wavelengths = wavelength*1.0E-6\n", - "psfcube = ns.calc_datacube(wavelengths, fov_pixels=30, oversample=4, add_distortion=True)\n", - "psfcube.writeto(\"Webbpsf_ifucube.fits\")\n", - "'''" + "# wavelengths = wavelength*1.0E-6\n", + "# psfcube = ns.calc_datacube(wavelengths, fov_pixels=30, oversample=4, add_distortion=True)\n", + "# psfcube.writeto(\"Webbpsf_ifucube.fits\")" ] }, { @@ -483,20 +518,26 @@ "outputs": [], "source": [ "BoxPath = \"https://data.science.stsci.edu/redirect/JWST/jwst-data_analysis_tools/IFU_optimal_extraction/\"\n", - "psf_filename = BoxPath+\"Webbpsf_ifucube.fits\"\n", + "psf_filename = BoxPath + \"Webbpsf_ifucube.fits\"\n", "\n", - "# Load with astropy.fits.open\n", + "# Open WebbPSF data cube\n", "with fits.open(psf_filename, memmap=False) as hdulist:\n", - " psf_model = hdulist['DET_SAMP'].data\n", - " psf_hdr = hdulist['DET_SAMP'].header\n", - " hdulist.info() \n", - "print(psf_model.shape)\n", + " psf_model = hdulist[\"DET_SAMP\"].data\n", + " psf_hdr = hdulist[\"DET_SAMP\"].header\n", + " hdulist.info()\n", + "\n", + "# Pad PSF model cube with zeros to match the present dataset\n", + "# (Different padding may be needed for your particular dataset)\n", + "print(sci.shape, psf_model.shape)\n", + "psf_model_padded = np.pad(psf_model, ((0, 0), (4, 5), (6, 7)), \"constant\")\n", "\n", "# Sum over wavelength\n", - "psf_model_sum = np.sum(psf_model, axis=0)\n", + "psf_model_sum = np.sum(psf_model_padded[slice_range[0]: slice_range[-1] + 1], axis=0)\n", "\n", "# Sum over spaxels\n", - "psf_model_fnusum = np.sum(psf_model, axis=(1, 2))" + "psf_model_fnusum = np.sum(\n", + " psf_model_padded[slice_range[0]: slice_range[-1] + 1], axis=(1, 2)\n", + ")" ] }, { @@ -511,15 +552,12 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "## Align Model PSF Cube with Science Data\n", - "Flip, smooth, and shift the model PSF cube to align with the simulated data. Trim the simulated data. " - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "*Developer Note:* Automate this by finding the (x,y) offset between the Model and simulated PSF peaks. Currently the shift is determined empirically by eye." + "## 11. Align Model PSF Cube with Science Data \n", + "Flip, smooth, and shift the model PSF cube to align with the simulated data. Trim the simulated data. \n", + "\n", + "Important Note 1: this PSF will likely be rotated with respect to your dataset, depending on telescope roll angle. You can either rotate it to match your data or reprocess your data using the ifualign keyword to align the WCS with the instrumental coordinate frame.\n", + "\n", + "Important Note 2: this PSF will likely be shifted with respect to your dataset. It would be beneficial to automatically find the (x,y) offset between the data and simulated PSF peaks. Currently the shift is determined empirically by eye." ] }, { @@ -528,56 +566,53 @@ "metadata": {}, "outputs": [], "source": [ - "# Flip model PSF left-right. For some unknown reason, WebbPSF is flipped with respect to the IPS simulation.\n", - "psf_model_fliplr = psf_model[:, ::-1, :]\n", - "\n", - "# Smooth model\n", - "# EMSM smoothing for G140H grating\n", - "scalerad = 0.046 # sigma (arcsec)\n", - "pixelscale = 0.1\n", - "scalerad_pix = scalerad / pixelscale\n", - "psf_model_smoothed = scipy.ndimage.filters.gaussian_filter(psf_model_fliplr, \n", - " (0.0, scalerad_pix, scalerad_pix), \n", - " order=0, mode='reflect', cval=0.0, \n", - " truncate=10.0)\n", + "# Flip model PSF left-right to match data.\n", + "psf_model_fliplr = psf_model_padded[:, ::-1, :]\n", "\n", "# Empirically (chi-by-eye) determined shift\n", - "shiftx = 1.75 \n", - "shifty = 0.\n", + "shiftx = 1.5 # 2\n", + "shifty = 1.5 # 1.5\n", "\n", "# Shift model PSF using linear interpolation\n", - "psf_model_aligned = scipy.ndimage.shift(psf_model_smoothed, (0.0, shiftx, shifty), order=1, \n", - " mode='constant', cval=0.0, prefilter=True)\n", + "psf_model_aligned = scipy.ndimage.shift(\n", + " psf_model_fliplr,\n", + " (0.0, shiftx, shifty),\n", + " order=1,\n", + " mode=\"constant\",\n", + " cval=0.0,\n", + " prefilter=True,\n", + ")\n", + "\n", + "good_psf_model = psf_model_aligned[\n", + " slice_range[0]: slice_range[-1] + 1,\n", + " data_yrange[0]: data_yrange[1],\n", + " data_xrange[0]: data_xrange[1],\n", + "]\n", "\n", "# Sum over wavelength\n", - "psf_model_sum = np.sum(psf_model_aligned, axis=0)\n", + "psf_model_sum = np.sum(good_psf_model, axis=0)\n", "\n", "# Scale factor for PSF subtraction\n", - "psf_sum_min = np.amin(psf_model_sum)\n", "psf_sum_max = np.amax(psf_model_sum)\n", "scalefactor = np.amax(cube_sum) / psf_sum_max\n", "\n", "# Plots\n", - "f, ([ax1, ax2, ax3], [ax4, ax5, ax6]) = plt.subplots(2, 3, figsize=(10, 10)) \n", + "f, ([ax1, ax2], [ax3, ax4]) = plt.subplots(2, 2, figsize=(10, 10))\n", "\n", "ax1.set_title(\"PSF slice sum\")\n", - "ax1.imshow(psf_model_sum, norm=LogNorm())\n", + "ax1.imshow(psf_model_sum, norm=LogNorm(), origin=\"lower\")\n", "\n", "ax2.set_title(\"Science Data slice sum\")\n", - "ax2.imshow(cube_sum, norm=LogNorm()) \n", + "ax2.imshow(cube_sum, norm=LogNorm(), origin=\"lower\")\n", "\n", "ax3.set_title(\"Data / PSF Ratio\")\n", - "ax3.imshow(cube_sum / psf_model_sum, norm=LogNorm())\n", + "ax3.imshow(cube_sum / psf_model_sum, norm=LogNorm(vmin=1, vmax=1e6), origin=\"lower\")\n", "\n", - "ax4.set_title(\"PSF Model integrated flux\")\n", - "ax4.plot(psf_model_fnusum)\n", - "\n", - "ax5.set_title(\"Data - PSF\")\n", - "ax5.imshow(cube_sum - scalefactor * psf_model_sum)\n", - "\n", - "im6 = ax6.imshow(np.log10(np.absolute(cube_sum - scalefactor * psf_model_sum)))\n", - "plt.colorbar(im6)\n", - "ax6.set_title(\"log abs(Data - PSF)\")\n", + "im4 = ax4.imshow(\n", + " np.log10(np.absolute(cube_sum - 0.75 * scalefactor * psf_model_sum)), origin=\"lower\"\n", + ")\n", + "plt.colorbar(im4)\n", + "ax4.set_title(\"log abs(Data - PSF)\")\n", "\n", "plt.show()" ] @@ -586,20 +621,19 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "_Figure top row_: Comparison of smoothed, aligned WebbPSF PSF (left) to IPS simulation (center). \n", + "_Figure top row_: Comparison of shifted WebbPSF PSF (left) to science data (right). The NIRSpec IFU PSF is signficantly different from the WebbPSF simulation for the telescope. This is partly due to the IFU optics and perhaps partly due to the cube building algorithm. There is also some real extended emission from the QSO host and surrounding galaies.\n", "\n", - "_Figure bottom row_: Integrated WebbPSF model flux (left) decreases with wavelength as PSF expands outside of the FOV. \n", - "Differences (center, right) between the model PSF and IPS-simulated PSF will translate to inaccuracy in the optimally-extracted spectrum." + "_Figure bottom row_: Differences between the model PSF and the observed PSF will result in a sub-optimal extraction, not quite attaining the maximum SNR, but still better than a sum over all spaxels." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "## Optimal Extraction using WebbPSF Model\n", - "Optimal extraction (Horne 1986, PASP, 98, 609) weights the flux contributions to a spectrum by their signal-to-noise ratio (SNR). Dividing the simulated data by the model PSF gives an estimate of the total flux density spectrum in each spaxel. A weighted average of these estimates over all spaxels yields the optimally extracted spectrum over the cube. In the faint source limit, where the noise is background-dominated, optimal extraction inside a 3-sigma radius can increase the effective exposure time by a factor of 1.69 (Horne et al. 1986). In the bright source limit, where the noise is dominated by the Poisson statistics of the source, optimal extraction is formally identical to a straight sum over spaxels for a perfect PSF model. \n", + "## 12. Optimal Extraction using WebbPSF Model \n", + "Optimal extraction ([Horne 1986, PASP, 98, 609](https://ui.adsabs.harvard.edu/abs/1986PASP...98..609H/abstract)) weights the flux contributions to a spectrum by their signal-to-noise ratio (SNR). Dividing the simulated data by the model PSF gives an estimate of the total flux density spectrum in each spaxel. A weighted average of these estimates over all spaxels yields the optimally extracted spectrum over the cube. In the faint source limit, where the noise is background-dominated, optimal extraction inside a 3-sigma radius can increase the effective exposure time by a factor of 1.69 (Horne et al. 1986). In the bright source limit, where the noise is dominated by the Poisson statistics of the source, optimal extraction is formally identical to a straight sum over spaxels for a perfect PSF model. \n", "\n", - "We use the WebbPSF PSF model for this first attempt at optimal extraction." + "We use the precomputed WebbPSF PSF model for optimal extraction here." ] }, { @@ -609,183 +643,75 @@ "outputs": [], "source": [ "# Window PSF model (and replace NaNs)\n", - "profile = np.nan_to_num(psf_model_aligned[0:2059, :, :]) \n", + "good_profile = np.nan_to_num(good_psf_model)\n", + "var_clean = np.nan_to_num(good_var, nan=1e12, posinf=1e12, neginf=1e12)\n", + "zerovar = np.where(var_clean == 0)\n", + "var_clean[zerovar] = 1e12\n", + "var_clean_sum = np.sum(var_clean, axis=(0))\n", + "snr_clean = np.nan_to_num(good_data / var_clean)\n", "\n", "# Divide data by PSF model\n", - "data_norm = np.nan_to_num(data_win / profile)\n", - "data_norm_sum = np.sum(data_norm, axis=0) \n", - "\n", - "# Mask out bad data using 3-sigma clipping in each slice\n", - "data_norm_clipped = sigma_clip(data_norm, sigma=3.0, maxiters=5, axis=(1, 2))\n", - "data_norm_clipped_sum = np.sum(data_norm_clipped, axis=0) \n", - "badvoxel = np.where(data_norm_clipped == 0)\n", - "data_clean = 1.0 * data_win\n", + "data_norm = np.nan_to_num(good_data / good_profile, posinf=0, neginf=0)\n", + "data_norm_sum = np.sum(data_norm, axis=0)\n", + "\n", + "# Mask out bad data\n", + "# data_norm_clipped = sigma_clip(data_norm, sigma=3.0, maxiters=5, axis=(1, 2))\n", + "data_norm_clipped = data_norm\n", + "data_norm_clipped_sum = np.sum(data_norm_clipped, axis=0)\n", + "snr_thresh = 1.0\n", + "badvoxel = np.where((data_norm_clipped == 0) | (snr_clean < snr_thresh))\n", + "data_clean = 1.0 * good_data\n", "data_clean[badvoxel] = 0.0\n", + "data_clean_sum = np.sum(data_clean, axis=0)\n", "\n", "# Optimal extraction, using model profile weight and variance cube from the simulated data\n", - "optimal_weight = profile ** 2 / data_var\n", + "optimal_weight = np.nan_to_num(\n", + " good_profile**2 / var_clean, posinf=0, neginf=0\n", + ") # Replace nans and infs with 0\n", + "optimal_weight_sum = np.sum(optimal_weight, axis=(0))\n", "optimal_weight_norm = np.sum(optimal_weight, axis=(1, 2))\n", - "spectrum_optimal = np.sum(profile * data_clean / data_var, axis=(1, 2)) / optimal_weight_norm\n", - "\n", - "opt_scalefactor = np.median(np.nan_to_num(cone_sum / spectrum_optimal)) # = 1.33, not ~1.0 because PSF model isn't perfect\n", + "spectrum_optimal = (\n", + " np.sum(good_profile * data_clean / var_clean, axis=(1, 2)) / optimal_weight_norm\n", + ")\n", "\n", "# Plots\n", - "f, (ax1) = plt.subplots(1, 1, figsize=(12, 6)) \n", + "f, (ax1) = plt.subplots(1, 1, figsize=(12, 6))\n", "ax1.set_title(\"Optimal Extraction Comparison\")\n", - "ax1.set_xlabel(\"Observed Wavelength (microns)\") \n", + "ax1.set_xlabel(\"Observed Wavelength (microns)\")\n", "ax1.set_ylabel(\"Flux Density\")\n", - "ax1.set_ylim(0, 0.5)\n", + "ax1.set_ylim(0, 5000)\n", "ax1.plot(wavelength, cone_sum, label=\"Conical Extraction\", alpha=0.5)\n", "ax1.plot(wavelength, spectrum_optimal, label=\"Optimal\")\n", "ax1.legend()\n", - "\n", - "plt.show()" - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "The optimally extracted spectrum is less noisy than the aperture extraction, but the flux density is low by a factor of ~1.33 because the PSF model doesn't match the science data perfectly." - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## Optimal Extraction with (Simulated) Reference Star PSF\n", - "A real (or simulated in this case) IFU observation of a star may be used for the PSF model rather than WebbPSF. We employ a NIRSpec IPS simulated PSF, which matches our data better than the WebbPSF model. We don't have to shift or smooth the PSF model because it was simulated at the same dither/detector position as the data. When using a real observation of a star for the PSF model, make sure it was observed at the same dither positions. It is also beneficial to reduce and extract both simulated datasets in the 'ifualign' detector coordinate system, so that we don't have to rotate the PSF star to match the science data." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "BoxPath = \"https://data.science.stsci.edu/redirect/JWST/jwst-data_analysis_tools/IFU_optimal_extraction/\"\n", - "filename_star = BoxPath + \"NRS00001-brightQSO-F100LP-G140H-01_1_491_SE_2020-08-26T12h15m00_s3d.fits\"\n", - "\n", - "# Open and inspect the file and WCS\n", - "with fits.open(filename_star, memmap=False) as hdulist:\n", - " sci_star = hdulist['SCI'].data\n", - " err_star = hdulist['ERR'].data\n", - " w_star = wcs.WCS(hdulist[1].header)\n", - " hdr_star = hdulist[1].header\n", - " hdulist.info()\n", - " print(w_star)\n", - " \n", - "# Load with Spectrum1D \n", - "spec1d_star = Spectrum1D.read(filename_star)\n", - "\n", - "# Wavelengths\n", - "wavelength_star = np.array(spec1d_star.spectral_axis.value)\n", - "\n", - "# Window reference star to match science data (and replace NaNs)\n", - "ref_star = np.nan_to_num(sci_star[:, 5:-4, 3:])\n", - "\n", - "# Sum over spaxels\n", - "ref_star_fnusum = np.sum(ref_star, axis=(1, 2))\n", - "\n", - "# Normalize PSF star profile to unity. (The flux will still be slightly off. Please see Developer's Note below.)\n", - "ref_star_norm = []\n", - "for idx, norm in zip(range(len(wavelength_star)), ref_star_fnusum):\n", - " ref_star_norm.append(ref_star[idx] / norm)\n", - "profile_star = np.array(ref_star_norm)\n", - " \n", - "# Sum over spaxels \n", - "profile_star_fnusum = np.sum(profile_star, axis=(1, 2))\n", - "\n", - "# Sum over wavelength\n", - "profile_star_sum = np.sum(profile_star, axis=0)\n", - "\n", - "# Scale factor for PSF subtraction\n", - "profile_star_sum_max = np.amax(profile_star_sum)\n", - "star_scalefactor = np.amax(cube_sum) / profile_star_sum_max\n", - "\n", - "# Make slight adjustment to scale factor\n", - "star_scalefactor = 0.175\n", - "\n", - "# Plots\n", - "f, ([ax1, ax2, ax3], [ax4, ax5, ax6]) = plt.subplots(2, 3, figsize=(10, 10)) \n", - "\n", - "ax1.imshow(profile_star_sum, norm=LogNorm())\n", - "ax1.set_title(\"PSF Star Slice sum\")\n", - "\n", - "ax2.imshow(cube_sum, norm=LogNorm()) \n", - "ax2.set_title(\"Science Data Slice sum\")\n", - "\n", - "ax3.imshow(cube_sum / profile_star_sum, norm=LogNorm())\n", - "ax3.set_title(\"Data/Star_PSF Ratio\")\n", - "\n", - "star_model_ratio = profile_star_sum / psf_model_sum\n", - "ax4.imshow(star_model_ratio, norm=LogNorm())\n", - "ax4.set_title(\"Star PSF/WebbPSF\")\n", - "\n", - "ax5.imshow(cube_sum - star_scalefactor * profile_star_sum)\n", - "ax5.set_title(\"Data - Star PSF\")\n", - "\n", - "ax6.imshow(np.log10(np.absolute(cube_sum - star_scalefactor * profile_star_sum)))\n", - "ax6.set_title(\"log abs(Data - Star PSF)\")\n", - "\n", "plt.show()" ] }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "_Figure top row_: Comparison of PSF star and science data. Bottom left: ratio of PSF star to WebbPSF model shows \n", - "significant differences that can affect the quality of the optimal extraction. Bottom right:\n", - "Difference of PSF star from science data shows they are well matched, with a scale factor of 0.175. " - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "*Developer Note:* It would be good to renormalize the PSF profile to account for the fraction of flux lost outside of the detector. Otherwise the extracted flux will be low by a factor of roughly 0.972 to 0.980." - ] - }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ - "# Mask out bad data using 3-sigma clipping in each slice\n", - "data_norm = np.nan_to_num(data_win / profile_star)\n", - "data_norm_clipped = sigma_clip(data_norm, sigma=3.0, maxiters=5, axis=(1, 2))\n", - "data_norm_clipped_sum = np.sum(data_norm_clipped, axis=0) \n", - "badvoxel = np.where(data_norm_clipped == 0)[0]\n", - "data_clean = 1.0 * data_win\n", - "data_clean[badvoxel] = 0.0\n", - "\n", - "# Optimal extraction, using model profile weight and variance cube from the simulated data\n", - "optimal_weight = profile_star**2 / data_var\n", - "optimal_weight_norm = np.sum(optimal_weight, axis=(1, 2))\n", - "spectrum_optimal_star = np.sum(profile_star * data_clean / data_var, axis=(1, 2)) / optimal_weight_norm\n", + "specviz = Specviz()\n", "\n", - "# Plots\n", - "f, (ax1) = plt.subplots(1, 1, figsize=(12, 6)) \n", - "ax1.set_title(\"Optimal Extraction Comparison\")\n", - "ax1.set_xlabel(\"Observed Wavelength (microns)\") \n", - "ax1.set_ylabel(\"Flux Density\")\n", - "ax1.set_ylim(0, 0.5)\n", + "flux_opt = spectrum_optimal * u.MJy / u.sr\n", + "spec1d_opt = Spectrum1D(spectral_axis=wavelength * u.um, flux=flux_opt)\n", "\n", - "ax1.plot(wavelength, cone_sum, label=\"Conical Extraction\", alpha=0.5)\n", - "ax1.plot(wavelength, spectrum_optimal * opt_scalefactor, label=\"1.3 * Optimal with WebbPSF model\", alpha=0.5)\n", - "ax1.plot(wavelength, spectrum_optimal_star, label=\"Optimal with ref. star\")\n", - "ax1.legend()\n", + "# specviz.load_data(spec1d_spaxsum, data_label=\"collapse spec\")\n", + "specviz.load_data(spec1d_opt, data_label=\"optimal spec\")\n", + "specviz.load_data(spec1d_cone, data_label=\"cone spec\")\n", + "specviz.show()\n", "\n", - "plt.show()" + "# set spectrum display limits\n", + "# specviz.x_limits()\n", + "specviz.y_limits(0.0, clip_level / 7)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "The optimal extraction with the perfectly matched PSF star is less noisy than that achieved with WebbPSF, and unlike the latter, doesn't need to be rescaled. The scaling can be off if the PSF of the reference star is not a good match to the PSF of the science data." + "The optimally extracted spectrum is less noisy than the aperture extraction and incorporates fewer bad pixels and cosmic ray events. The OIII line profile is different because extended emission is downweighted with respect to the unresolved quasar nucleus." ] }, { @@ -799,15 +725,8 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "Notebook created by Patrick Ogle and James Davies." + "Notebook created by Patrick Ogle and James Davies. Last update: December 2023." ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [] } ], "metadata": { diff --git a/notebooks/ifu_optimal/pre-requirements.txt b/notebooks/ifu_optimal/pre-requirements.txt deleted file mode 100644 index dc9fd9ab0..000000000 --- a/notebooks/ifu_optimal/pre-requirements.txt +++ /dev/null @@ -1,2 +0,0 @@ -astropy-helpers>=2.0.11 - diff --git a/notebooks/ifu_optimal/requirements.txt b/notebooks/ifu_optimal/requirements.txt index d338557ed..fc4aff4ff 100644 --- a/notebooks/ifu_optimal/requirements.txt +++ b/notebooks/ifu_optimal/requirements.txt @@ -1,17 +1,3 @@ -notebook>=6.1.5 -jupyter_client>=5.3.5 -numpy>=1.19.2 -scipy>=1.10.0 -specutils>=1.1.1 -glue-astronomy>=0.1 -glue-core>=1.0.0 -glue-jupyter>=0.2.1 -glue-vispy-viewers>=1.0.1 -jdaviz>=1.0.3 -regions>=0.7 -photutils>=1.0.1 -astropy>=4.2 -matplotlib>=3.3.2 -ipyvuetify>=1.6.2 -ipyvue>=1.5.0 -radio-beam>=0.3.6 +numpy>=1.26.2 +scipy>=1.11.4 +jdaviz>=3.8.0 diff --git a/notebooks/ifu_optimal/sdssj1652_nirspec_ifu_cubeviz.png b/notebooks/ifu_optimal/sdssj1652_nirspec_ifu_cubeviz.png new file mode 100644 index 000000000..e75819210 Binary files /dev/null and b/notebooks/ifu_optimal/sdssj1652_nirspec_ifu_cubeviz.png differ