diff --git a/.circleci/config.yml b/.circleci/config.yml index 1ecfe313..343d8de4 100644 --- a/.circleci/config.yml +++ b/.circleci/config.yml @@ -271,24 +271,6 @@ jobs: no_output_timeout: 120m - store_artifacts: path: /tmp/nbcollection-ci-artifacts - Publish Website: - executor: nbcollection-builder - steps: - - checkout - - run: - command: bash ./.circleci/setup-env.sh - name: Setup Environment - no_output_timeout: 120m - - run: - command: nbcollection-ci merge-artifacts -c jdat_notebooks -o $CIRCLE_PROJECT_USERNAME - -r $CIRCLE_PROJECT_REPONAME - name: Publish Website - no_output_timeout: 120m - - add_ssh_keys - - run: - command: nbcollection-ci site-deployment -r origin -b gh-pages - name: Deploy Website - no_output_timeout: 120m Pull Request: executor: nbcollection-builder steps: @@ -326,25 +308,4 @@ workflows: - Jdat Notebooks-Nircam Photometry - Jdat Notebooks-Optimal Extraction Dynamic - Pull Request - - Publish Website: - requires: - - Jdat Notebooks-Niriss Wfss Postpipeline - - Jdat Notebooks-Mos-Spectroscopy - - Jdat Notebooks-Mrs Mstar Analysis - - Jdat Notebooks-Cube Fitting - - Jdat Notebooks-Ifu Cube Continuum Fit - - Jdat Notebooks-Nircam Psf-Matched Photometry - - Jdat Notebooks-Miri Ifu Ysos In The Lmc - - Jdat Notebooks-Optimal Extraction - - Jdat Notebooks-Ifu Optimal - - Jdat Notebooks-Background Estimation Imaging - - Jdat Notebooks-Psf Photometry - - Jdat Notebooks-Specviz Notebookgui Interaction - - Jdat Notebooks-Asdf Example - - Jdat Notebooks-Aperture Photometry - - Jdat Notebooks-Nirspec Mast Query - - Jdat Notebooks-Miri Lrs Spectral Extraction - - Jdat Notebooks-Transit Spectroscopy Notebook - - Jdat Notebooks-Nircam Photometry - - Jdat Notebooks-Optimal Extraction Dynamic version: '2.1' diff --git a/.circleci/setup-env.sh b/.circleci/setup-env.sh index 70a87ddd..37bddce8 100644 --- a/.circleci/setup-env.sh +++ b/.circleci/setup-env.sh @@ -4,7 +4,8 @@ set -e export LANG=C.UTF-8 export LC_ALL=C.UTF-8 -conda install python=3.7.11 + +conda install python=3.8.10 apt-get update apt-get install -y git make build-essential libssl-dev zlib1g-dev libbz2-dev libreadline-dev libsqlite3-dev wget curl \ llvm libncurses5-dev xz-utils tk-dev libxml2-dev libxmlsec1-dev libffi-dev liblzma-dev libfreetype6-dev @@ -14,4 +15,4 @@ cd nbcollection pip install -U pip setuptools pip install -r ci_requirements.txt python setup.py install -cd - \ No newline at end of file +cd - diff --git a/jdat_notebooks/IFU_cube_continuum_fit/NGC4151_FeII_ContinuumFit.ipynb b/jdat_notebooks/IFU_cube_continuum_fit/NGC4151_FeII_ContinuumFit.ipynb index c324c190..de90a782 100644 --- a/jdat_notebooks/IFU_cube_continuum_fit/NGC4151_FeII_ContinuumFit.ipynb +++ b/jdat_notebooks/IFU_cube_continuum_fit/NGC4151_FeII_ContinuumFit.ipynb @@ -4,7 +4,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "# IFU Data Modeling, NGC 4151 Notebook #1 - Isolating Line Emission\n" + "# IFU Data Modeling\n" ] }, { @@ -15,7 +15,7 @@ "**Data:** [NIFS](https://www.gemini.edu/instrumentation/current-instruments/nifs) on Gemini; NGC 4151.
\n", "**Tools:** specutils, jdaviz/cubeviz, astropy, matplotlib, bottleneck.
\n", "**Cross-intrument:** NIRSpec; potentially 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", + "**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) and can be [downloaded](https://github.com/spacetelescope/jdat_notebooks/tree/main/notebooks/IFU_cube_continuum_fit/) directly from the [JDAT Notebook Github directory](https://github.com/spacetelescope/jdat_notebooks).
\n", "\n", "## Introduction\n", "\n", @@ -25,7 +25,7 @@ "\n", "**Note:** This notebook is designed to analyze the 1.6440 [Fe II] emission but the wavelengths can be altered to fit and remove continuum around any emission line of interest.\n", "\n", - "## Imports\n", + "## Import Packages\n", " - time for timing\n", " - numpy for array processing and math\n", " - matplotlib.pyplot for plotting images and spectra\n", @@ -44,22 +44,24 @@ "metadata": {}, "outputs": [], "source": [ - "\n", "# load important packages\n", "import time\n", + "import os\n", "from copy import copy\n", - "\n", - "import numpy as np\n", + "from IPython.display import HTML, YouTubeVideo\n", "\n", "import astropy\n", + "import numpy as np\n", "from astropy.io import fits, ascii\n", "from astropy import units as u\n", "from astropy.modeling import models\n", "from astropy.utils.data import download_file\n", "from specutils.fitting import fit_lines\n", "from specutils import Spectrum1D\n", - "\n", - "from jdaviz.app import Application\n" + "from jdaviz import Cubeviz\n", + "from jdaviz.app import Application\n", + "from specutils.manipulation import extract_region\n", + "from specutils.spectra import SpectralRegion" ] }, { @@ -68,11 +70,30 @@ "metadata": {}, "outputs": [], "source": [ - "\n", "# load and configure matplotlib\n", "%matplotlib inline\n", "import matplotlib.pyplot as plt\n", - "plt.rcParams.update({'figure.max_open_warning': 0})\n" + "plt.rcParams.update({'figure.max_open_warning': 0})" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Save and Load Objects Using Pickle\n", + "import pickle\n", + "\n", + "\n", + "def save_obj(obj, name):\n", + " with open(name, 'wb') as f:\n", + " pickle.dump(obj, f, pickle.HIGHEST_PROTOCOL)\n", + "\n", + "\n", + "def load_obj(name):\n", + " with open(name, 'rb') as f:\n", + " return pickle.load(f)" ] }, { @@ -83,8 +104,7 @@ }, "outputs": [], "source": [ - "\n", - "#This cell access the datacube file, defines the wavelength grid from header information and then plots a simple\n", + "# This cell accesses the datacube file, defines the wavelength grid from header information and then plots a simple\n", "# 1-D collapsed spectrum of the IFU data.\n", "\n", "# Read in a 3-D IFU datacube of interest, and header.\n", @@ -93,25 +113,25 @@ "cube = fits.getdata(cube_file)\n", "header_cube = fits.getheader(cube_file)\n", "\n", - "#grab data information and wavelength definitions.\n", + "# grab data information and wavelength definitions.\n", "nz, ny, nx = cube.shape\n", "crdelt3 = header_cube['CDELT3']\n", "crval3 = header_cube['CRVAL3']\n", "\n", - "#define the wavelength grid (microns) from the header (Angstroms)\n", + "# define the wavelength grid (microns) from the header (Angstroms)\n", "# and the AGN redshift and the emission line of interest.\n", - "wave =((crdelt3 * (np.arange(0,nz,1))) + crval3)/10000.0\n", + "wave = ((crdelt3 * (np.arange(0, nz, 1))) + crval3)/10000.0\n", "redshift = 0.00332\n", "emission_line = 1.64400*(1 + redshift)\n", "emission_line_index = (np.abs(wave-emission_line)).argmin()\n", "\n", "# make a simple summed 1d spectrum of the full cube\n", - "flux1 = np.sum(cube, axis=(1,2))\n", + "flux1 = np.sum(cube, axis=(1, 2))\n", "\n", "# plot the full 1-D spectrum.\n", "plt.figure(0)\n", "plt.plot(wave, flux1)\n", - "plt.show()\n" + "plt.show()" ] }, { @@ -139,9 +159,8 @@ }, "outputs": [], "source": [ - "\n", - "#This cell defines the wavelength regions of interest: around the emission line, and the location\n", - "#where you want to fit and remove the continuum very accurately. Make a plot that shows the regions.\n", + "# This cell defines the wavelength regions of interest: around the emission line, and the location\n", + "# where you want to fit and remove the continuum very accurately. Make a plot that shows the regions.\n", "\n", "# Here we select a region that includes the emission line\n", "# wavelength plus a small range of continuum around it. \n", @@ -162,7 +181,7 @@ "continuum_limit1 = 1.656\n", "continuum_limit2 = 1.673\n", " \n", - "#Define the wavelength region around the emission - indices\n", + "# Define the wavelength region around the emission - indices\n", "wavemin = (np.abs(wave-wave_emission_limit1)).argmin()\n", "wavemax = (np.abs(wave-wave_emission_limit2)).argmin()\n", "\n", @@ -170,23 +189,178 @@ "continuummin = (np.abs(wave-continuum_limit1)).argmin()\n", "continuummax = (np.abs(wave-continuum_limit2)).argmin()\n", "\n", - "#show the region used for the emission line and continuum fit. Alter the wavelengths \n", + "# Show the region used for the emission line and continuum fit. Alter the wavelengths \n", "# above if this doesn't look good. \n", "plt.figure(1)\n", "plt.plot(wave, flux1)\n", "plt.plot(wave[wavemin:wavemax], flux1[wavemin:wavemax])\n", - "plt.plot(wave[continuummin:continuummax], flux1[continuummin:continuummax],color='r')\n", - "plt.show()\n" + "plt.plot(wave[continuummin:continuummax], flux1[continuummin:continuummax], color='r')\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Cubeviz Visualization\n", + "You can also visualize images inside a Jupyter notebook using [Cubeviz](https://jdaviz.readthedocs.io/en/latest/cubeviz/index.html)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Video1: \n", + "\n", + "This Cubeviz Instructional Demo is from STScI's official YouTube channel and provides an introduction to Cubeviz." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "vid = YouTubeVideo(\"HMSYwiH3Gl4\")\n", + "display(vid)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "cubeviz = Cubeviz()\n", + "cubeviz.app" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Here, we load the data into the Cubeviz app.\n", + "cubeviz.load_data(fn) " + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "A video is shown below illustrating the procedure. The following steps are applied:\n", + "\n", + "When you load your cube, you will see a collapsed spectrum of all the spaxels in the spectral viewer at the bottom.\n", + "\n", + "If you draw a region (circle or square) in the flux viewer, you will see a collapsed spectrum of that particular region in the spectral viewer, too. To do this, use the expandable menu at the left of the flux viewer window and select the 'define circular region of interest' icon. For this example, we want to first define a circular region at the central position of the bright AGN flux, which is at approximately the cube center position.\n", + "\n", + "Now, use the flux viewer and again use the 'define circular region of interest' icon to make spectra at two positions associated with the outflow emission in [Fe II]. The redshifted outflow is at approximate x position = 12, y position = 36. This will be 'Subset 2' and will show up in green in the display. The blueshifted outflow is at approximately x position = 50, y position = 28 in pixel index units. This will be 'Subset 3' and will show up in blue in the display.\n", + "\n", + "(If the notebook is being run non-interactively, automatically make two datasets that mimic the AGN outflow red/blueshifted spectra)." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Defining your Spectral Regions\n", + "\n", + "Next, you will want to define the wavelengths of interest in your spectral viewer for both your line and continuum analysis. To do this, you will similarly click the 'define region of interest' icon in your spectral viewer and drag a box over the wavelengths you desire. Again, a video is show below illustrating the process.\n", + "\n", + "There is no option to set the spectral regions to user input, so we recommend zooming in and drawing by eye. The line emission (;Subset 4' ) should span approximately 1.630 - 1.665 um, and the continuum emission ('Subset 5') should span approximately 1.656 - 1.673 um." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Some Notes:\n", + "\n", + "*If your cell window requires you to scroll to see the different displays in cubeviz, you can toggle the scroll window in the main menu of the notebook: Cell -> Current Outputs -> Toggle Scrolling\n", + "\n", + "*In the datacube viewing panel, you can select the 'layer' tab (horizontal line settings icon) in the upper panel within the viewer and change the display scaling. Decreasing the maximum display value by 10x brings out the low level extended emission in this dataset. In this cube, data from slice ~1060 to ~1090 shows the extended [Fe II] emission. " + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Video2: \n", + "\n", + "Here is a video illustrating how to load and manipulate data for this particular notebook. A description of how to manipulate the data and use the \"Collapse\" tool is given above." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "IFrame(\"https://data.science.stsci.edu/redirect/JWST/jwst-data_analysis_tools/IFU_cube_continuum_fit/cubeviz1.mov\", width=700, height=500)\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "HTML('')" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Video3:\n", + "\n", + "Here is a video illustrating how to define your regions of interest and extract them into the notebook for more analysis below." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# PART 1\n", + "HTML('')" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# PART 2\n", + "HTML('')" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Extract Subset Spectrum in Cubeviz Spectrum Viewer" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - " For this particular dataset, this continuum region looks very good.\n", - " if you have a more structured continuum you can define additional\n", - " regions and append them into a larger set of wave / flux arrays to \n", - " derive a more accurate fit in the below poly-fit analysis." + "Retrieve the spectrum (Subset1) of the user-defined region from the Spectrum Viewer as a Spectrum1D object." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Extract spectra corresponding to the colored regions in cubeviz\n", + "spectrum1 = cubeviz.app.get_data_from_viewer('spectrum-viewer', 'Subset 1') # AGN Center\n", + "spectrum2 = cubeviz.app.get_data_from_viewer('spectrum-viewer', 'Subset 2') # Red shifted component\n", + "spectrum3 = cubeviz.app.get_data_from_viewer('spectrum-viewer', 'Subset 3') # Blue shifted component\n", + "spectrum1" ] }, { @@ -197,10 +371,13 @@ }, "outputs": [], "source": [ - "\n", - "# Use cubeviz to look at the data in the cube. Here, we configure and bring up the app.\n", - "app = Application(configuration='cubeviz')\n", - "app\n" + "# Extract the line region defined in the spectral viewer\n", + "regions = cubeviz.specviz.get_spectral_regions()\n", + " \n", + "if regions and \"Subset 4\" in regions.keys():\n", + " line_region = regions[\"Subset 4\"]\n", + "else:\n", + " line_region = SpectralRegion(1.630*u.um, 1.665*u.um)" ] }, { @@ -211,25 +388,98 @@ }, "outputs": [], "source": [ - "\n", - "# Here, we load the data into the cubeviz app.\n", - "app.load_data(fn)\n", - " " + "# Extract the continuum region defined in the spectral viewer\n", + "if regions and \"Subset 5\" in regions.keys():\n", + " continuum_region = regions[\"Subset 5\"]\n", + "else:\n", + " continuum_region = SpectralRegion(1.656*u.um, 1.673*u.um)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Apply the spectral region\n", + "# (creates new collapsed spectra if user did not in jdaviz)\n", + "if not spectrum1:\n", + " flux_agn = np.sum(cube[:, (ny//2)-3:(ny//2)+3, (nx//2)-3:(nx//2)+3], axis=(1, 2))\n", + " tmpspec = Spectrum1D(flux=flux_agn*u.Unit('count'), spectral_axis=wave*u.micron) \n", + " spec_agn = extract_region(tmpspec, line_region)\n", + " spec_agn_continuum = extract_region(tmpspec, continuum_region) \n", + "else: \n", + " spec_agn = extract_region(spectrum1, line_region)\n", + " spec_agn_continuum = extract_region(spectrum1, continuum_region)\n", + "\n", + "if not spectrum2:\n", + " flux_feii_red = np.sum(cube[:, (36)-3:(36)+3, (12)-3:(12)+3], axis=(1, 2))\n", + " tmpspec = Spectrum1D(flux=flux_feii_red*u.Unit('count'), spectral_axis=wave*u.micron) \n", + " spec_feii_red = extract_region(tmpspec, line_region)\n", + " spec_feii_red_continuum = extract_region(tmpspec, continuum_region)\n", + "else: \n", + " spec_feii_red = extract_region(spectrum2, line_region)\n", + " spec_feii_red_continuum = extract_region(spectrum2, continuum_region)\n", + "\n", + "if not spectrum3:\n", + " flux_feii_blue = np.sum(cube[:, (28)-3:(28)+3, (50)-3:(50)+3], axis=(1, 2))\n", + " tmpspec = Spectrum1D(flux=flux_feii_blue*u.Unit('count'), spectral_axis=wave*u.micron) \n", + " spec_feii_blue = extract_region(tmpspec, line_region)\n", + " spec_feii_blue_continuum = extract_region(tmpspec, continuum_region)\n", + "else: \n", + " spec_feii_blue = extract_region(spectrum3, line_region)\n", + " spec_feii_blue_continuum = extract_region(spectrum3, continuum_region)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Visualize new subsets\n", + "plt.figure()\n", + "plt.plot(spec_agn.spectral_axis, spec_agn.flux, color='black')" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Visualize new subsets\n", + "plt.figure()\n", + "plt.plot(spec_feii_blue.spectral_axis, spec_feii_blue.flux, color='b')\n", + "plt.plot(spec_feii_red.spectral_axis, spec_feii_red.flux, color='r')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "\n", - "Set the left cubeviz panel to view our cube by linking the dataset to the visualization. This is done in the cube visualization panel that is to the upper left; it is an expandible menu that shows up as lines at the left display view. When you click the gear icon in this expanded menu, the dataset should show up as a weird name next to a blank checkbox. Check the box.\n", - "\n", - "In the same way, select the datacube in the spectral viewer at the bottom of the cubeviz pane.\n", - "Use the expandible menu gear icon and select the loaded datacube. This is the summed 1-D spectrum from the full spatial field of the data cube. \n", - "\n", - "In the datacube viewing panel, you can also select the 'layer' tab in the gear (data) icon and change \n", - "display scaling. Decreasing the maximum display value by 10x brings out the low level extended emission\n", - "in this dataset. In this cube, data from slice ~1060 to ~1090 shows the extended [Fe II] emission. " + "## Fit the Continuum at the Spectral Region Location\n", + "\n", + "A video is shown below illustrating the procedure for fitting the continuum for 'Subset 5' throughout the entire cube. The following steps are applied to this particular example:
\n", + "\n", + "Open up Model Fitting Plugin. There are a number of fields to fill in and drop down menus to select from. It is important to keep in mind that the Data menu will provide only spectra to model, while the Spectral Region menu will provide only spectral region subsets to choose. In other words, you can fit the spectra in specific spectral regions. If no spectral region is selected, the entire wavelength array will be fit by the mode.
\n", + "\n", + "Select Data Cube: contentsSCI
\n", + "Select Spectral Region Subset 5
\n", + "Model: Linear1D
\n", + "ModelID: L1
\n", + "Model Parameters: Leave Default
\n", + "Model Equation Editor: L1
\n", + "Model Label: LinFitCont
\n", + "
\n", + "Hit Fit, which fits the collapsed spectrum.
\n", + "View the fit in the spectral viewer and confirm you are happy with it.
\n", + "Then hit Apply to Cube.
\n", + "\n", + "This will create two models that can now be accessed within the Data Dropdown menus:
\n", + "A 1D linear fit of the continuum for the collapsed cube.
\n", + "A 3D linear fit of the continuum for each spaxel in the cube.
\n", + "\n" ] }, { @@ -238,28 +488,33 @@ "metadata": {}, "outputs": [], "source": [ - "# show the 1-D spectrum by grabbing the cubeviz spectral viewer app within the notebook.\n", - "# If the 1-D spectrum hasn't been selected in the cubeviz window, automatically collapse\n", - "# and plot the loaded data cube.\n", - "\n", - "spectrum_viewer = app.get_viewer('spectrum-viewer')\n", - "if not spectrum_viewer.data():\n", - " app.add_data_to_viewer('spectrum-viewer',app.data_collection[0].label)\n", - "# spectrum_viewer.add_data(app.data_collection[0].label)\n", - " \n", - "spectrum_viewer.show()" + "# VIDEO OF CONTINUUM FITTING\n", + "# PART 1\n", + "HTML('')" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# PART 2\n", + "HTML('')" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "Use the flux viewer in cubeviz to extract a spectrum located at the central AGN position.\n", - "It should show in the above spectral viewer, too.\n", + "## Pulling Data from Viewers\n", "\n", - "To do this, use the expandable menu at the left of the flux viewer window and select the 'define circular region of interest' icon. Make a circular region at the central position of the bright AGN flux, which is at approximately the cube center position.\n", + "Note, in cubeviz, you have 4 viewers from which you can pull data. Make sure your data are properly loaded into each viewer before executing the get_data_from_viewer command:
\n", "\n", - "(If the notebook is being run non-interactively, this cell will automatically make a dataset that mimics the AGN specrum)." + "Top Left: flux-viewer
\n", + "Center: uncert-viewer
\n", + "Top Right: mask-viewer
\n", + "Bottom: spectrum-viewer
" ] }, { @@ -268,33 +523,26 @@ "metadata": {}, "outputs": [], "source": [ - "# Grab the layer spectrum that corresponds to the central AGN position as an array in the notebook, then plot it.\n", - "# If there is no subset in cubeviz, make one to plot.\n", - "\n", - "spec_agn = app.get_data_from_viewer('spectrum-viewer', 'Subset 1')\n", - "spec_agn\n", - "\n", - "if not spec_agn:\n", - " flux_agn = np.sum(cube[:,(ny//2)-3:(ny//2)+3,(nx//2)-3:(nx//2)+3], axis=(1,2))\n", - " spec_agn = Spectrum1D(flux=flux_agn*u.Unit('count'), spectral_axis=wave*u.micron) \n", - " \n", - "# plot the 1-D spectrum of this smaller sub-region.\n", - "plt.figure(2)\n", - "plt.plot(wave[wavemin:wavemax], spec_agn.flux[wavemin:wavemax])\n", - "plt.show()" + "# List data from viewer\n", + "regions = cubeviz.app.get_data_from_viewer(\"uncert-viewer\")\n", + "regions" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Extract continuum model from Cubeviz above\n", + "cont_psf_cube = cubeviz.app.get_data_from_viewer(\"uncert-viewer\", \"LinFitCont [Cube] 1\")" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "Now, use the flux viewer and again use the 'define circular region of interest' icon to make spectra at two positions associated with the outflow emission in [Fe II].\n", - "\n", - "The redshifted outflow is at approximate x position = 12, y position = 36. This will be 'Subset 2' and will show up in green in the display.\n", - "\n", - "The blueshifted outflow is at approximately x position = 50, y position = 28 in pixel index units. This will be 'Subset 3' and will show up in blue in the display.\n", - "\n", - "(If the notebook is being run non-interactively, automatically make two datasets that mimic the AGN outflow red/blueshifted spectra)." + "## Important Note: Always save data and variables when possible." ] }, { @@ -303,34 +551,16 @@ "metadata": {}, "outputs": [], "source": [ - "\n", - "#grab the redshifted and blueshifted outflow spectra from the spectral viewer window.\n", - "\n", - "spec_feii_red = app.get_data_from_viewer('spectrum-viewer', 'Subset 2')\n", - "spec_feii_red\n", - "\n", - "if not spec_feii_red:\n", - " flux_feii_red = np.sum(cube[:,(36)-3:(36)+3,(12)-3:(12)+3], axis=(1,2))\n", - " spec_feii_red = Spectrum1D(flux=flux_feii_red*u.Unit('count'), spectral_axis=wave*u.micron) \n", - "\n", - "spec_feii_blue = app.get_data_from_viewer('spectrum-viewer', 'Subset 3')\n", - "spec_feii_blue\n", - "\n", - "if not spec_feii_blue:\n", - " flux_feii_blue = np.sum(cube[:,(28)-3:(28)+3,(50)-3:(50)+3], axis=(1,2))\n", - " spec_feii_blue = Spectrum1D(flux=flux_feii_blue*u.Unit('count'), spectral_axis=wave*u.micron) \n", - "\n", - "#plot a zoomed view of the spectra grabbed from cubeviz\n", - "\n", - "# plot the 1-D redshifted outflow spectrum.\n", - "plt.figure(3)\n", - "plt.plot(wave[wavemin:wavemax], spec_feii_blue.flux[wavemin:wavemax])\n", - "plt.show()\n", - "\n", - "# plot the 1-D spectrum blueshifted outflow spectrum.\n", - "plt.figure(3)\n", - "plt.plot(wave[wavemin:wavemax], spec_feii_red.flux[wavemin:wavemax], color='r')\n", - "plt.show()\n" + "# Delete any existing output in current directory\n", + "if os.path.exists(\"NGC4151_Hband_ContinuumSubtract.fits\"):\n", + " os.remove(\"NGC4151_Hband_ContinuumSubtract.fits\")\n", + "else:\n", + " print(\"The file does not exist\")\n", + "\n", + "if os.path.exists(\"NGC4151_Hband_ContinuumPSF.fits\"):\n", + " os.remove(\"NGC4151_Hband_ContinuumPSF.fits\")\n", + "else:\n", + " print(\"The file does not exist\")" ] }, { @@ -339,25 +569,46 @@ "metadata": {}, "outputs": [], "source": [ + "# Subtract Continuum\n", "\n", - "#demonstration of a linear fit to the continuum flux level\n", - "# this method uses simple functions in numpy to fit the continuum. \n", - "\n", - "cont_fit = np.polyfit(wave[continuummin:continuummax], flux1[continuummin:continuummax], 1)\n", - "fitval = np.poly1d(cont_fit)\n", - "continuum = fitval(wave)\n", - "\n", - "plt.figure(4)\n", - "plt.plot(wave, flux1)\n", - "plt.show()\n", - "\n", - "plt.figure(4)\n", - "plt.plot(wave[continuummin:continuummax], flux1[continuummin:continuummax])\n", - "plt.show()\n", - "\n", - "plt.figure(4)\n", - "plt.plot(wave, continuum)\n", - "plt.show()\n" + "# Re-read in original IFU cube for manipulation\n", + "cube_file = 'https://data.science.stsci.edu/redirect/JWST/jwst-data_analysis_tools/IFU_cube_continuum_fit/NGC4151_Hband.fits'\n", + "newfn = download_file('https://data.science.stsci.edu/redirect/JWST/jwst-data_analysis_tools/IFU_cube_continuum_fit/NGC4151_Hband.fits', cache=True)\n", + "newheader_cube = fits.getheader(cube_file)\n", + "\n", + "# Check to see if user made a continuum fit in Cubeviz, make continuum subtraction, and save output\n", + "if not cont_psf_cube:\n", + " start_time = time.time()\n", + "\n", + " cont_sub_cube = np.zeros([nz, ny, nx])\n", + " cont_psf_cube = np.zeros([nz, ny, nx])\n", + "\n", + " for i in range(1, nx-2):\n", + " for j in range(1, ny-2):\n", + " flux1 = cube[:, j, i] \n", + " cont_fit = np.polyfit(wave[continuummin:continuummax], flux1[continuummin:continuummax], 1)\n", + " fitval = np.poly1d(cont_fit)\n", + " continuum = fitval(wave) \n", + " cont_sub_cube[:, j, i] = flux1 - continuum\n", + " cont_psf_cube[:, j, i] = continuum \n", + "\n", + " del newheader_cube['MODE']\n", + " fits.writeto('NGC4151_Hband_ContinuumSubtract.fits', cont_sub_cube, newheader_cube, overwrite=True)\n", + " fits.writeto('NGC4151_Hband_ContinuumPSF.fits', cont_psf_cube, newheader_cube, overwrite=True)\n", + " print('Continuum subtracted cube saved. PSF continuum cube saved.')\n", + "else:\n", + " with fits.open(newfn, memmap=False) as cont_sub_cube:\n", + " sci = cont_sub_cube['SCI'].data\n", + "\n", + " # Get List of different viewers\n", + " continuumflux = cont_psf_cube[\"flux\"]\n", + "\n", + " sci_contsub = sci-continuumflux\n", + " cont_sub_cube['SCI'].data = sci_contsub \n", + " del cont_sub_cube['PRIMARY'].header['MODE']\n", + " cont_sub_cube.writeto('NGC4151_Hband_ContinuumSubtract.fits')\n", + " del newheader_cube['MODE']\n", + " fits.writeto('NGC4151_Hband_ContinuumPSF.fits', continuumflux, newheader_cube, overwrite=True)" ] }, { @@ -366,44 +617,23 @@ "metadata": {}, "outputs": [], "source": [ - "# That looks pretty good.\n", - "\n", - "# This cell continues a demonstration of a linear fit to the continuum flux level.\n", - "# This method uses the specutils and modeling packages available\n", - "# in astropy. It does the same thing as the prior cell, but in an astropy \n", - "# spectral format.\n", - "\n", - "# Here we use the polynomial routine to do a linear fit - this is so it's easy\n", - "# to update the fit order, if necessary.\n", - "\n", - "full_spectrum=Spectrum1D(\n", - " flux=flux1* u.Unit('count'), spectral_axis=wave* u.micron)\n", - "\n", - "# Make Spectrum1D specutils version of data for the continuum segment.\n", - "spectrum = Spectrum1D(\n", - " flux=flux1[continuummin:continuummax]*u.Unit('count'),\n", - " spectral_axis=wave[continuummin:continuummax]*u.micron)\n", - "\n", - "# Make an empty model with no initial guess of the coefficients\n", - "# If guessing make sure to only pass values (no units)\n", - "m = models.Polynomial1D(degree=1) # You can place a guess by setting c0 and c1\n", - "# Fit model and save fitted model containing the fitted params\n", - "fitted_model = fit_lines(spectrum, m)\n", - "# Just to showcase how to access the fitted params\n", - "cont_fit = [fitted_model.c1.value, fitted_model.c0.value]\n", - "\n", - "# notice I dont have to use cont_fit to get my continuum\n", - "# I just call the model with my spectral values\n", - "continuum = fitted_model(full_spectrum.spectral_axis)\n", - "\n", - "#plot the results - identical to above numpy code - looks good.\n", - "plt.figure(5)\n", - "plt.plot(wave, flux1)\n", - "plt.show()\n", - "\n", - "plt.figure(5)\n", - "plt.plot(wave, continuum)\n", - "plt.show()\n" + "# You can also read out your model fit parameters \n", + "params = cubeviz.get_model_parameters(model_label=\"LinFitCont\")\n", + "params" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "*(The following cell is just an example of how you can acces the model fit parameter values):*" + ] + }, + { + "cell_type": "raw", + "metadata": {}, + "source": [ + "params['LinFitCont_3d']['slope']" ] }, { @@ -412,37 +642,9 @@ "metadata": {}, "outputs": [], "source": [ - "\n", - "# This cell makes a data sub-cube around the emission feature that has the continuum flux \n", - "# level subtracted off. then make another datacube that is only the continuum flux, to serve\n", - "# as a PSF model to correct out the bright, central HI features.\n", - "\n", - "# Here I'm using the numpy functions for continuum fitting and subtraction.\n", - "\n", - "# This is done in nested for loops, looping over the spatial axes in the cube.\n", - "\n", - "start_time = time.time()\n", - "\n", - "cont_sub_cube=np.zeros([nz,ny,nx])\n", - "cont_psf_cube=np.zeros([nz,ny,nx])\n", - "\n", - "for i in range(1, nx-2):\n", - " for j in range(1, ny-2):\n", - " flux1 = cube[:,j,i] \n", - " cont_fit = np.polyfit(wave[continuummin:continuummax], flux1[continuummin:continuummax], 1)\n", - " fitval = np.poly1d(cont_fit)\n", - " continuum = fitval(wave) \n", - " cont_sub_cube[:,j,i]= flux1 - continuum\n", - " cont_psf_cube[:,j,i]= continuum \n", - "\n", - "del header_cube['MODE']\n", - "\n", - "fits.writeto('NGC4151_Hband_ContinuumSubtract.fits', cont_sub_cube, header_cube, overwrite=True)\n", - "fits.writeto('NGC4151_Hband_ContinuumPSF.fits', cont_psf_cube, header_cube, overwrite=True)\n", - "print('Continuum subtracted cube saved. PSF continuum cube saved.')\n", - "\n", - "print('Time count')\n", - "print(\"--- %s seconds ---\" % (time.time() - start_time))\n" + "# Open up a new instance of Cubeviz to visualize continuum subtracted data\n", + "cubeviz2 = Cubeviz()\n", + "cubeviz2.app" ] }, { @@ -451,34 +653,45 @@ "metadata": {}, "outputs": [], "source": [ + "cont_sub_cube = 'NGC4151_Hband_ContinuumSubtract.fits'\n", + "cubeviz2.app.load_data(cont_sub_cube, data_label='Continuum Subtracted')" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Fitting your multiple component Gaussian Model\n", "\n", - "# This cell does the same thing as the prior cell, but uses specutils models polyfit instead of numpy \n", - "# to do the continuum fitting. This uses the poly fit options instead of line fitting so that you can \n", - "# change the polynomial order for the continuum fitting, if necessary.\n", + "Now we want to investigate an initial fit to the Br 12 emission feature, which is a pesky contaminant nearby in wavelength to our target [Fe II] emission. The Br 12 is centrally compact and arises from only from the nucleus of the AGN, not from the outflow. Make a plot of the fit results.
\n", "\n", - "# This is done in nested for loops, looping over the RA, Dec axes in the cube.\n", + "A video is shown below illustrating the procedure. The following steps are applied:
\n", "\n", - "# This does the same thing but is slow compared to the prior cell...\n", + "First, select the wavelength region of interest following a similar procedure as performed at the top. There is no option to set the spectral regions to a user input, so we recommend zooming in and drawing by eye. The line emission ('Subset 1' ) should again span approximately 1.630 - 1.665 um.\n", "\n", - "start_time = time.time()\n", + "For this example, we recommend setting up a 3 component gaussian model with the following inputs
\n", + "Open up Model Fitting Plugin. There are a number of fields to fill in and drop down menus to select from. It is important to keep in mind that the Data menu will provide only spectra to model, while the Spectral Region menu will provide only spectral region subsets to choose. In other words, you can fit the spectra in specific spectral regions. If no spectral region is selected, the entire wavelength array will be fit by the mode.
\n", "\n", - "cont_sub_cube_specutils=np.zeros([nz,ny,nx])\n", - "cont_psf_cube_specutils=np.zeros([nz,ny,nx])\n", + "Data: Continuum Subtracted
\n", + "Spectral region or subset covering about 1.632 to 1.656 um
\n", + "Model: Three different Gaussians with ModelID's set to G1, G2, and G3
\n", + "Model Parameters:
\n", "\n", - "for i in range(1, nx-2):\n", - " for j in range(1, ny-2):\n", - " flux1 = Spectrum1D(flux = cube[:,j,i]*u.Unit('count'), spectral_axis=wave*u.micron)\n", - " m = models.Polynomial1D(degree=1)\n", - " fitted_model = fit_lines(flux1[continuummin:continuummax], m)\n", - " cont_fit = [fitted_model.c1.value, fitted_model.c0.value]\n", - " continuum = fitted_model(flux1.spectral_axis)\n", - " cont_sub_cube_specutils[:,j,i]= flux1.flux - continuum\n", - " cont_psf_cube_specutils[:,j,i]= continuum\n", + "G1: stdev=0.0008E-6, mean=1.641E-6
\n", + "G2: stdev=0.0007E-6, mean=1.648E-6
\n", + "G3: stdev=0.005E-6, mean=1.646E-6
\n", + "You can turn on the 'Fixed' option if you need to, but these numbers should provide a good starting guess for the fit.
\n", "\n", - "print('Done')\n", + "Model Equation Editor: G1+G2+g3
\n", + "Model Label: GaussAll
\n", + "
\n", + "Hit Fit, which fits the collapsed spectrum.
\n", + "View the fit in the spectral viewer and confirm you are happy with it. Modify if necessary.
\n", + "Then hit Apply to Cube.
\n", "\n", - "print('Time count')\n", - "print(\"--- %s seconds ---\" % (time.time() - start_time))\n" + "This will again create two models that can now be accessed within the Data Dropdown menus:
\n", + "A 1D linear fit of the lines in the collapsed cube.
\n", + "A 3D linear fit of the lines for each spaxel in the cube.
" ] }, { @@ -487,76 +700,27 @@ "metadata": {}, "outputs": [], "source": [ - "\n", - "# Now, let's take a look at the emission line in the continuum by plotting some of the results\n", - "# of the subtracted cube.\n", - "flux1=np.sum(cont_sub_cube, axis=(1,2))\n", - "\n", - "# plot the 1-D spectrum of the full continuum subtracted cube\n", - "plt.figure(6)\n", - "plt.plot(wave[wavemin:wavemax], flux1[wavemin:wavemax])\n", - "plt.show()\n", - "\n", - "# plot the 1-D spectrum of a single spaxel.\n", - "plt.figure(7)\n", - "plt.plot(wave[wavemin:wavemax], cont_sub_cube[wavemin:wavemax,30,30])\n", - "plt.show()\n" + "# VIDEO PART1\n", + "HTML('')" ] }, { "cell_type": "code", "execution_count": null, - "metadata": { - "scrolled": true - }, + "metadata": {}, "outputs": [], "source": [ + "# VIDEO PART2\n", + "HTML('')" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Wow, that multi-component fit looks great. Good deal.\n", "\n", - "# This cell builds the spectrum at the central location into the format\n", - "# needed for use with specutils. It then investigates an initial fit to the\n", - "# Br 12 emission feature, which is a pesky contaminant nearby in wavelength\n", - "# to our target [Fe II] emission. The Br 12 is centrally compact and arises from\n", - "# only from the nucleus of the AGN, not from the outflow. Make a plot of the fit\n", - "# results.\n", - "\n", - "#zoom in wavelength into the region of interest, create subcube and subwave arrays.\n", - "flux = (cont_sub_cube[wavemin:wavemax,30,30])\n", - "minwave = wave[wavemin:wavemax]\n", - "\n", - "# put the flux spectrum into the spec utils expected format.\n", - "spectrum = Spectrum1D(flux=(flux)*u.Unit('count'), spectral_axis=minwave*u.micron)\n", - "\n", - "#define the fit the line for the Brackett emission (position was found by hand @ pix 1023):\n", - "# the central emission is best fit by two gaussian components: one @ br12, one @ [Fe II].\n", - "# Here we fit a third component too: the [Fe II] outflow emission.\n", - "l1 = models.Gaussian1D(amplitude = (flux[1023-wavemin])*u.Unit('count'), mean = minwave[1023-wavemin]*u.micron, stddev = 0.0009*u.micron)\n", - "l2 = models.Gaussian1D(amplitude = (flux[emission_line_index-wavemin])*u.Unit('count'), mean = minwave[emission_line_index-wavemin]*u.micron, stddev = 0.005*u.micron)\n", - "#define and fit the line for the outflow [Fe II] emission:\n", - "l3 = models.Gaussian1D(amplitude = (flux[emission_line_index-wavemin])*u.Unit('count'), mean = minwave[emission_line_index-wavemin]*u.micron, stddev = 0.0008*u.micron)\n", - "\n", - "#run the lfit - this tweaks the above parameters to optimize the fits of the three components.\n", - "lfit = fit_lines(spectrum, l1 + l2 + l3)\n", - "#make the yfit\n", - "y_fit = lfit(minwave*u.micron)\n", - "\n", - "# Build the fits from the fit_lines function into specutils format for plotting.\n", - "lineflux = (lfit[0](minwave*u.micron))\n", - "linemodel = Spectrum1D(spectral_axis=minwave*u.micron, flux=lineflux*u.Unit('count'))\n", - "\n", - "component1 = lfit[0](minwave*u.micron)\n", - "component2 = lfit[1](minwave*u.micron)\n", - "component3 = lfit[2](minwave*u.micron)\n", - "\n", - "plt.figure(8)\n", - "plt.plot(minwave, flux)\n", - "plt.plot(minwave, component1)\n", - "plt.plot(minwave, component2)\n", - "plt.plot(minwave, component3)\n", - "plt.plot(minwave, component1 + component2 + component3)\n", - "plt.show()\n", - "\n", - "#we want to isolate just the [Fe II] outflow emission, so subtract off the central compact flux sources\n", - "central_flux_model_only = component1 + component2\n" + "Now we're going to use the continuum psf cube from a prior cell with the Brackett model created in the above cell to create a full 3-D model of the central emission that isn't caused by the outflow [Fe II]." ] }, { @@ -565,59 +729,259 @@ "metadata": {}, "outputs": [], "source": [ + "# Extract the spectral regions defined in the spectral viewer\n", + "regions = cubeviz2.specviz.get_spectral_regions()\n", "\n", - "# Wow, that multi-component fit looks great. Good deal.\n", - "\n", - "#now we're going to use the continuum psf cube from a prior cell \n", - "# with the Brackett model created in the above cell to create a full\n", - "# 3-D model of the central emission that isn't caused by the outflow [Fe II].\n", - "\n", - "continuum_subcube = cont_psf_cube[wavemin:wavemax,:,:]\n", - "nz, ny, nx = continuum_subcube.shape\n", - "\n", - "model_cube=np.zeros([nz,ny,nx])\n", - "\n", - "#construct the scaled Brackett flux model\n", - "model_cube[0,:,:] = continuum_subcube[0,:,:] * (central_flux_model_only[0]/continuum_subcube[0, 30, 30])\n", - "for i in range(1, nz-2):\n", - " model_cube[i,:,:] = continuum_subcube[i,:,:] * (central_flux_model_only[i] / continuum_subcube[i, 30, 30])\n", - "model_cube[nz-1,:,:] = continuum_subcube[nz-1,:,:] * (central_flux_model_only[nz-1] / continuum_subcube[nz-1,30,30])\n", - "\n", - "# the full model of the AGN central emission is the continuum plus Brackett line.\n", - "full_model = continuum_subcube + model_cube\n", - "\n", - "# subtract the model to create the final cube where the [Fe II] emission\n", - "# is isolated.\n", - "final_sub_cube = cube[wavemin:wavemax,:,:] - full_model\n", + "if regions and \"Subset 1\" in regions.keys():\n", + " line_region = regions[\"Subset 1\"]\n", + "else:\n", + " line_region = SpectralRegion(1.6322*u.um, 1.6563*u.um)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# List spectra available in spectrum-viewer\n", + "spec = cubeviz2.app.get_data_from_viewer('spectrum-viewer') \n", + "spec" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Get gauss model spectrum and model cube\n", + "all_spec = cubeviz2.app.get_data_from_viewer('spectrum-viewer', 'Continuum Subtracted[SCI]') # AGN Center Data Cube\n", + "gauss_spec = cubeviz2.app.get_data_from_viewer('spectrum-viewer', 'GaussAll') # AGN Center Model Spec\n", + "gauss_cube = cubeviz2.app.get_data_from_viewer(\"uncert-viewer\", \"GaussAll [Cube] 1\") # AGN Center Model Cube" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Check to see if user used Cubeviz (above), and, if not, read in premade data\n", + "if not gauss_cube:\n", + " fn = download_file('https://data.science.stsci.edu/redirect/JWST/jwst-data_analysis_tools/IFU_cube_continuum_fit/gauss_model_cube.fits', cache=False)\n", + " gauss_cube = fits.getdata(fn)\n", + "else:\n", + " gauss_cube = gauss_cube[\"flux\"]\n", + " \n", + "if not all_spec:\n", + " fn = download_file('https://data.science.stsci.edu/redirect/JWST/jwst-data_analysis_tools/IFU_cube_continuum_fit/all_spec.fits', cache=False)\n", + " all_spec = Spectrum1D.read(fn)\n", + " \n", + "if not gauss_spec:\n", + " fn = download_file('https://data.science.stsci.edu/redirect/JWST/jwst-data_analysis_tools/IFU_cube_continuum_fit/gauss_spec.fits', cache=False)\n", + " gauss_spec = Spectrum1D.read(fn)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Save your gauss model cube if necessary\n", + "if os.path.exists(\"gauss_model_cube.fits\"):\n", + " os.remove(\"gauss_model_cube.fits\")\n", + "else:\n", + " print(\"The file does not already exist\")\n", + "\n", + "fits.writeto('gauss_model_cube.fits', gauss_cube, overwrite=True)\n", + "print(type(gauss_cube))\n", + "print(type(all_spec.flux))" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "*(Developer Note: if saving with 'Count' doesn't work, try the following code.)*" + ] + }, + { + "cell_type": "raw", + "metadata": {}, + "source": [ + "print(all_spec.flux)\n", + "os.remove(\"all_spec.fits\")\n", + "specref=Spectrum1D(flux=all_spec.flux*u.Unit('Jy'), spectral_axis=all_spec.spectral_axis)\n", + "specref.write(\"all_spec.fits\")\n", + "\n", + "os.remove(\"gauss_spec.fits\")\n", + "specref=Spectrum1D(flux=gauss_spec.flux*u.Unit('Jy'), spectral_axis=gauss_spec.spectral_axis)\n", + "specref.write(\"gauss_spec.fits\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# You can also read out your model fit parameters \n", + "params = cubeviz2.get_model_parameters(model_label=\"GaussAll\")\n", + "params" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Save parameters as a pickle file if necessary\n", + "save_obj(params, \"gauss_params.pkl\")\n", + "params = load_obj(\"gauss_params.pkl\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Check to see if user used Cubeviz (above), and, if not, read in premade data\n", + "if not params:\n", + " fn = download_file('https://data.science.stsci.edu/redirect/JWST/jwst-data_analysis_tools/IFU_cube_continuum_fit/gauss_params.pkl', cache=True)\n", + " params = load_obj(fn)\n", + " print(\"Loaded\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Overwrite gauss model with only 2 of the components of interest\n", + "gauss_cube_2component = gauss_cube*0\n", + "\n", + "nz, ny, nx = gauss_cube_2component.shape\n", + "for i in range(0, nx-1):\n", + " for j in range(0, ny-1):\n", + " amp1 = params['GaussAll_3d']['amplitude_0'][i][j]\n", + " amp2 = params['GaussAll_3d']['amplitude_2'][i][j]\n", + " m1 = params['GaussAll_3d']['mean_0'][i][j]\n", + " m2 = params['GaussAll_3d']['mean_2'][i][j]\n", + " stdev1 = params['GaussAll_3d']['stddev_0'][i][j]\n", + " stdev2 = params['GaussAll_3d']['stddev_2'][i][j]\n", + " g1 = models.Gaussian1D(amplitude=amp1*u.Unit('count'), mean=m1*u.m, stddev=stdev1*u.m)\n", + " g2 = models.Gaussian1D(amplitude=amp2*u.Unit('count'), mean=m2*u.m, stddev=stdev2*u.m)\n", + " gauss_cube_2component[:, i, j] = g1(all_spec.spectral_axis)+g2(all_spec.spectral_axis)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Add the continuum cube to the new model cube\n", + "continuum_file = 'NGC4151_Hband_ContinuumPSF.fits'\n", + "newfull_header = fits.getheader(continuum_file)\n", "\n", - "# make an appropriate header for the output sub-cube\n", - "header_cube_small = copy(header_cube)\n", - "del header_cube_small['CRVAL3']\n", - "header_cube_small['CRVAL3'] = wave[wavemin] * 10000.0\n", - "del header_cube_small['CRPIX3']\n", - "header_cube_small['CRPIX3'] = 1\n", + "with fits.open(continuum_file, memmap=False) as continuum_cube: \n", + " continuum_data = continuum_cube[0].data\n", + " full_model = gauss_cube_2component+continuum_data" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Subtract the model to create the final cube where the [Fe II] emission is isolated.\n", + "# Re-read in original IFU cube for manipulation\n", + "cube_file = 'https://data.science.stsci.edu/redirect/JWST/jwst-data_analysis_tools/IFU_cube_continuum_fit/NGC4151_Hband.fits'\n", + "newfinalsub_header = fits.getheader(cube_file)\n", "\n", - "# Save the .fits data sub-cube that has the continuum and Br model subtracted off of the\n", - "# [Fe II] emission, and the datacube that is the continuum+Br model.\n", - "fits.writeto('NGC4151_Hband_FinalSubtract.fits', final_sub_cube, header_cube_small, overwrite=True)\n", - "fits.writeto('NGC4151_Hband_ContinuumandBrackettModel.fits', full_model, header_cube_small, overwrite=True)\n", - "print('Continuum and Brackett subtracted cube saved. Full model cube saved.')\n", + "with fits.open(cube_file, memmap=False) as original_cube:\n", + " original_data = original_cube['SCI'].data\n", + " final_sub_cube = original_data - full_model" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Delete any existing output in current directory\n", + "if os.path.exists(\"NGC4151_Hband_FinalSubtract.fits\"):\n", + " os.remove(\"NGC4151_Hband_FinalSubtract.fits\")\n", + "else:\n", + " print(\"The file does not exist\")\n", + "\n", + "if os.path.exists(\"NGC4151_Hband_ContinuumandBrackettModel.fits\"):\n", + " os.remove(\"NGC4151_Hband_ContinuumandBrackettModel.fits\")\n", + "else:\n", + " print(\"The file does not exist\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "del newfinalsub_header['MODE']\n", "\n", - "#make a plot of the central spectrum, the full model and the continuum.\n", - "plt.figure(9)\n", - "plt.plot(minwave, continuum_subcube[:,30,30])\n", - "plt.plot(minwave, cube[wavemin:wavemax,30,30])\n", - "plt.plot(minwave, full_model[:,30,30])\n", - "plt.show()\n" + "fits.writeto('NGC4151_Hband_ContinuumandBrackettModel.fits', full_model, newfull_header, overwrite=True)\n", + "fits.writeto('NGC4151_Hband_FinalSubtract.fits', final_sub_cube, newfinalsub_header, overwrite=True)\n", + "print('Continuum subtracted cube saved. PSF continuum cube saved.')" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Make the final plots to illustrate the original spectrum, model fits, and final continuum+gassian subtracted cube\n", + "plt.figure()\n", + "plt.xlim([1.630E-6, 1.665E-6])\n", + "plt.ylim([600, 900])\n", + "plt.plot(all_spec.spectral_axis, continuum_data[:, 30, 30], label='Continuum')\n", + "plt.plot(all_spec.spectral_axis, original_data[:, 30, 30], label='Original Data')\n", + "plt.plot(all_spec.spectral_axis, full_model[:, 30, 30], label='2 Component Model')\n", + "plt.plot(all_spec.spectral_axis, final_sub_cube[:, 30, 30]+700, label='Model Subtraction+Offset')\n", + "plt.legend()\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Open up a new instance of Cubeviz to visualize continuum subtracted data\n", + "cubeviz3 = Cubeviz()\n", + "cubeviz3.app" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "cont_sub_cube = 'NGC4151_Hband_FinalSubtract.fits'\n", + "cubeviz3.app.load_data(cont_sub_cube, data_label='Red/Blue Shift')" ] } ], "metadata": { - "kernelspec": { - "display_name": "Python 3", - "language": "python", - "name": "python3" - }, "language_info": { "codemirror_mode": { "name": "ipython", @@ -628,7 +992,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.7.1" + "version": "3.8.10" } }, "nbformat": 4, diff --git a/jdat_notebooks/IFU_cube_continuum_fit/requirements.txt b/jdat_notebooks/IFU_cube_continuum_fit/requirements.txt index cbfc4525..669df510 100644 --- a/jdat_notebooks/IFU_cube_continuum_fit/requirements.txt +++ b/jdat_notebooks/IFU_cube_continuum_fit/requirements.txt @@ -1,7 +1,7 @@ -astropy>=4.1 -specutils>=1.0 -matplotlib>=3.2 +astropy>=4.3.1 +specutils>=1.4.0 +matplotlib>=3.4.3 Bottleneck>=1.3.2 -ipyvue>=1.3.4 -ipyvuetify>=1.4.1 -git+https://github.com/spacetelescope/jdaviz.git@625acca8562249c2f4d9041c566c1820a4d480b2#egg=jdaviz +ipyvue>=1.5.0 +ipyvuetify>=1.8.1 +git+https://github.com/spacetelescope/jdaviz.git@main