diff --git a/notebooks/MIRI/MRS_Mstar_analysis/JWST_Mstar_dataAnalysis_analysis.ipynb b/notebooks/MIRI/MRS_Mstar_analysis/JWST_Mstar_dataAnalysis_analysis.ipynb index 09594f6ce..866a536ed 100644 --- a/notebooks/MIRI/MRS_Mstar_analysis/JWST_Mstar_dataAnalysis_analysis.ipynb +++ b/notebooks/MIRI/MRS_Mstar_analysis/JWST_Mstar_dataAnalysis_analysis.ipynb @@ -17,7 +17,7 @@ "source": [ "**Use case:** Extract spatial-spectral features from IFU cube and measure their attributes.
\n", "**Data:** Simulated [MIRI MRS](https://jwst-docs.stsci.edu/mid-infrared-instrument/miri-observing-modes/miri-medium-resolution-spectroscopy) spectrum of AGB star.
\n", - "**Tools:** specutils, jwst, photutils, astropy, scipy.
\n", + "**Tools:** specutils, astropy, scipy.
\n", "**Cross-intrument:** NIRSpec, MIRI.
\n", "**Documentation:** This notebook is part of a STScI's larger [post-pipeline Data Analysis Tools Ecosystem](https://jwst-docs.stsci.edu/jwst-post-pipeline-data-analysis) and can be [downloaded](https://github.com/spacetelescope/dat_pyinthesky/tree/main/jdat_notebooks/MRS_Mstar_analysis) directly from the [JDAT Notebook Github directory](https://github.com/spacetelescope/jdat_notebooks).
\n", "**Source of Simulations:** [MIRISim](https://www.stsci.edu/jwst/science-planning/proposal-planning-toolbox/mirisim)
\n", @@ -94,7 +94,7 @@ }, "outputs": [], "source": [ - "# Import astropy packages \n", + "# Import astropy packages\n", "from astropy import units as u\n", "from astropy.io import ascii\n", "from astropy.nddata import StdDevUncertainty\n", @@ -110,7 +110,7 @@ "from jdaviz import Cubeviz\n", "\n", "# Display the video\n", - "from IPython.display import HTML, YouTubeVideo" + "from IPython.display import YouTubeVideo" ] }, { @@ -127,7 +127,7 @@ " with open(name, 'wb') as f:\n", " pickle.dump(obj, f, pickle.HIGHEST_PROTOCOL)\n", "\n", - " \n", + "\n", "def load_obj(name):\n", " with open(name, 'rb') as f:\n", " return pickle.load(f)" @@ -140,7 +140,7 @@ "outputs": [], "source": [ "def checkKey(dict, key):\n", - " \n", + "\n", " if key in dict.keys():\n", " print(\"Present, \", end=\" \")\n", " print(\"value =\", dict[key])\n", @@ -169,25 +169,28 @@ "# Check if Pipeline 3 Reduced data exists and, if not, download it\n", "import os\n", "import urllib.request\n", + "import tarfile\n", "\n", "if os.path.exists(\"combine_dithers_all_exposures_ch1-long_s3d.fits\"):\n", " print(\"Pipeline 3 Data Exists\")\n", "else:\n", " url = 'https://data.science.stsci.edu/redirect/JWST/jwst-data_analysis_tools/MRS_Mstar_analysis/reduced.tar.gz'\n", " urllib.request.urlretrieve(url, './reduced.tar.gz')\n", - " # Unzip Tar Files\n", "\n", - " import tarfile\n", + " base_extract_to = os.path.abspath(\".\") # Current directory\n", "\n", - " # Unzip files if they haven't already been unzipped\n", - " if os.path.exists(\"reduced/\"):\n", - " print(\"Pipeline 3 Data Exists\")\n", - " else:\n", - " tar = tarfile.open('./reduced.tar.gz', \"r:gz\")\n", - " tar.extractall()\n", - " tar.close()\n", - " \n", - " # Move Files \n", + " # Open and securely extract files from the tar archive\n", + " with tarfile.open('./reduced.tar.gz', \"r:gz\") as tar:\n", + " for member in tar.getmembers():\n", + " # Calculate the absolute path of where the file will be extracted\n", + " member_path = os.path.abspath(os.path.join(base_extract_to, member.name))\n", + "\n", + " # Check if the file path is within the base extraction directory\n", + " if member_path.startswith(base_extract_to):\n", + " # Extract only safe files, directly to the base directory\n", + " tar.extract(member, path=base_extract_to)\n", + " else:\n", + " print(f\"Skipped {member.name} due to potential security risk\")\n", " os.system('mv reduced/*fits .')" ] }, @@ -297,7 +300,7 @@ "metadata": {}, "outputs": [], "source": [ - "# Load in the spectrum list from above. \n", + "# Load in the spectrum list from above.\n", "specviz.load_data(splist)" ] }, @@ -352,7 +355,7 @@ "metadata": {}, "outputs": [], "source": [ - "# Here, we load the data into the Cubeviz app for visual inspection. \n", + "# Here, we load the data into the Cubeviz app for visual inspection.\n", "# In this case, we're just looking at a single channel because, unlike Specviz, Cubeviz can only load a single cube at a time.\n", "\n", "ch1short_cubefile = 'combine_dithers_all_exposures_ch1-long_s3d.fits'\n", @@ -363,26 +366,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "Next, you want to define a pixel region subset that is specific to the AGB star. You can do this with the regions utility button shown in the video below and drawing a circular region around the AGB star at approximate pixels x=20, y=30." - ] - }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## Video 3:\n", - " \n", - "Here is a video that quickly shows how to select a spatial region in Cubeviz." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "metadata": {}, - "outputs": [], - "source": [ - "# Video showing the selection of the star with a circular region of interest\n", - "HTML('')" + "Next, you want to define a pixel region subset that is specific to the AGB star. You can do this with the regions utility button and drawing a circular region around the AGB star at approximate pixels x=20, y=30." ] }, { @@ -392,16 +376,16 @@ "outputs": [], "source": [ "# Now extract spectrum from your spectral viewer\n", + "# NEED TO show how to use Spectral Extraction plugin and calculate mean instead of sum spectra\n", "try:\n", - " spec_agb = cubeviz.get_data('combine_dithers_all_exposures_ch1-long_s3d[SCI]',\n", - " spatial_subset='Subset 1', function='mean') # AGB star only\n", + " spec_agb = cubeviz.get_data('Spectrum (Subset 1, sum)') # AGB star only\n", " print(spec_agb)\n", " spec_agb_exists = True\n", "except Exception:\n", " print(\"There are no subsets selected.\")\n", " spec_agb_exists = False\n", - " spec_agb = cubeviz.get_data('combine_dithers_all_exposures_ch1-long_s3d[SCI]',\n", - " function='mean') # Whole field of view" + " spec_agb = cubeviz.get_data('Spectrum (sum)') # Whole field of view\n", + " print(spec_agb)" ] }, { @@ -424,12 +408,13 @@ "metadata": {}, "outputs": [], "source": [ - "wav = wlall*u.micron # Wavelength: microns\n", - "fl = fnuall*u.Jy # Fnu: Jy\n", - "efl = dfnuall*u.Jy # Error flux: Jy\n", + "wav = wlall*u.micron # Wavelength: microns\n", + "fl = fnuall*u.Jy # Fnu: Jy\n", + "efl = dfnuall*u.Jy # Error flux: Jy\n", "\n", "# Make a 1D spectrum object\n", - "spec = Spectrum1D(spectral_axis=wav, flux=fl, uncertainty=StdDevUncertainty(efl))" + "spec = Spectrum1D(spectral_axis=wav, flux=fl,\n", + " uncertainty=StdDevUncertainty(efl))" ] }, { @@ -439,9 +424,9 @@ "outputs": [], "source": [ "# Apply a 5 pixel boxcar smoothing to the spectrum\n", - "spec_bsmooth = box_smooth(spec, width=5) \n", + "spec_bsmooth = box_smooth(spec, width=5)\n", "\n", - "# Plot the spectrum & smoothed spectrum to inspect features \n", + "# Plot the spectrum & smoothed spectrum to inspect features\n", "plt.figure(figsize=(8, 4))\n", "plt.plot(spec.spectral_axis, spec.flux, label='Source')\n", "plt.plot(spec.spectral_axis, spec_bsmooth.flux, label='Smoothed')\n", @@ -450,11 +435,14 @@ "plt.ylim(-0.05, 0.15)\n", "\n", "# Overplot the original input spectrum for comparison\n", - "origspecfile = fn = download_file('https://data.science.stsci.edu/redirect/JWST/jwst-data_analysis_tools/MRS_Mstar_analysis/63702662.txt', cache=True)\n", + "origspecfile = fn = download_file(\n", + " 'https://data.science.stsci.edu/redirect/JWST/jwst-data_analysis_tools/MRS_Mstar_analysis/63702662.txt', cache=True)\n", "origdata = ascii.read(origspecfile)\n", "wlorig = origdata['col1']\n", - "fnujyorig = origdata['col2']*0.001 # comes in as mJy, change to Jy to compare with pipeline output\n", - "plt.plot(wlorig, fnujyorig, '.', color='grey', markersize=1, label='Original Input')\n", + "# comes in as mJy, change to Jy to compare with pipeline output\n", + "fnujyorig = origdata['col2']*0.001\n", + "plt.plot(wlorig, fnujyorig, '.', color='grey',\n", + " markersize=1, label='Original Input')\n", "\n", "plt.legend(frameon=False, fontsize='medium')\n", "plt.tight_layout()\n", @@ -482,26 +470,13 @@ "specviz.show()" ] }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "#### Developer Note: Cannot currently open a spectrum1d output from cubeviz in specviz. https://jira.stsci.edu/browse/JDAT-1791\n", - "\n", - "#specviz.load_data(spec_agb)" - ] - }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ - "# For now, you must create your own spectrum1d object from your extracted cubeviz spectrum. \n", - "flux = spec_agb.flux\n", - "wavelength = spec_agb.spectral_axis\n", - "spec1d = Spectrum1D(spectral_axis=wavelength, flux=flux)\n", - "specviz.load_data(spec1d)" + "specviz.load_data(spec_agb)" ] }, { @@ -512,23 +487,13 @@ "specviz.load_data(spec)" ] }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## Video 4:\n", - " \n", - "Here is a video that quickly shows how to smooth your spectrum in Specviz." - ] - }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ - "# Video showing how to smooth a spectrum in Specviz\n", - "HTML('')" + "# Make new video to show how to smooth spectrum in Specviz" ] }, { @@ -554,23 +519,13 @@ "For now switching to a blackbody." ] }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## Video 5:\n", - " \n", - "Here is a video that shows how to fit a blackbody model to the spectrum" - ] - }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ - "# Video showing how to fit a blackbody \n", - "HTML('')" + "# Make new video to show how to fit a blackbody model to the spectrum" ] }, { @@ -580,14 +535,15 @@ "outputs": [], "source": [ "spectra = specviz.get_spectra()\n", - " \n", + "\n", "a = checkKey(spectra, \"BB1\")\n", "if a is True:\n", " # Extract Blackbody fit from Specviz\n", " blackbody = spectra[\"BB1\"]\n", "else:\n", " print(\"No Blackbody\")\n", - " fn = download_file('https://data.science.stsci.edu/redirect/JWST/jwst-data_analysis_tools/MRS_Mstar_analysis/blackbody.fits', cache=True)\n", + " fn = download_file(\n", + " 'https://data.science.stsci.edu/redirect/JWST/jwst-data_analysis_tools/MRS_Mstar_analysis/blackbody.fits', cache=True)\n", " blackbody = Spectrum1D.read(fn)" ] }, @@ -667,7 +623,8 @@ "\n", "# Now subtract the BB and plot the underlying dust continuum\n", "plt.figure(figsize=(8, 4))\n", - "plt.plot(spec.spectral_axis, spec.flux.value - ybest.value, color='purple', label='Dust spectra')\n", + "plt.plot(spec.spectral_axis, spec.flux.value -\n", + " ybest.value, color='purple', label='Dust spectra')\n", "plt.axhline(0, color='r', linestyle='dashdot', alpha=0.5)\n", "plt.xlabel('Wavelength (microns)')\n", "plt.ylabel(\"Flux ({:latex})\".format(spec.flux.unit))\n", @@ -745,23 +702,13 @@ "specviz.load_data(bbsub_spectra)" ] }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## Video 6:\n", - " \n", - "Here is a video that shows how to fit a polynomial to two separate spectral regions within a single subset to remove more underlying continuum" - ] - }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ - "# Video showing how to fit a polynomial to two separate spectral regions within a single subset\n", - "HTML('')" + "# Make new video to show how to fit a polynomial to two separate spectral regions within a single subset" ] }, { @@ -771,14 +718,15 @@ "outputs": [], "source": [ "spectra = specviz.get_spectra()\n", - " \n", + "\n", "a = checkKey(spectra, \"PolyFit\")\n", "if a is True:\n", " # Extract polynomial fit from Specviz\n", " poly = spectra[\"PolyFit\"]\n", "else:\n", " print(\"No Polyfit\")\n", - " fn = download_file('https://data.science.stsci.edu/redirect/JWST/jwst-data_analysis_tools/MRS_Mstar_analysis/poly.fits', cache=True)\n", + " fn = download_file(\n", + " 'https://data.science.stsci.edu/redirect/JWST/jwst-data_analysis_tools/MRS_Mstar_analysis/poly.fits', cache=True)\n", " poly = Spectrum1D.read(fn)" ] }, @@ -828,8 +776,10 @@ "# -----------------------------------------------------------------\n", "# Generate a continuum subtracted and continuum normalised spectra\n", "\n", - "line_spec_norm = Spectrum1D(spectral_axis=line_spec.spectral_axis, flux=line_spec.flux/line_y_continuum, uncertainty=StdDevUncertainty(np.zeros(len(line_spec.spectral_axis))))\n", - "line_spec_consub = Spectrum1D(spectral_axis=line_spec.spectral_axis, flux=line_spec.flux - line_y_continuum, uncertainty=StdDevUncertainty(np.zeros(len(line_spec.spectral_axis))))\n", + "line_spec_norm = Spectrum1D(spectral_axis=line_spec.spectral_axis, flux=line_spec.flux /\n", + " line_y_continuum, uncertainty=StdDevUncertainty(np.zeros(len(line_spec.spectral_axis))))\n", + "line_spec_consub = Spectrum1D(spectral_axis=line_spec.spectral_axis, flux=line_spec.flux -\n", + " line_y_continuum, uncertainty=StdDevUncertainty(np.zeros(len(line_spec.spectral_axis))))\n", "\n", "# -----------------------------------------------------------------\n", "# Plot the dust feature & continuum fit to the region\n", @@ -886,25 +836,13 @@ "specviz.load_data(line_spec_norm, data_label='Normalized')" ] }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "## Video7:\n", - " \n", - "Here is a video that shows how to make line analysis measurements within specviz (i.e., line flux, line centroid, equivalent width)
\n", - "Note: You want to calculate your equivalent width on the normalized spectrum
\n", - "Note: You can also hack to convert the line flux value into more conventional units
" - ] - }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ - "# Video showing how to measure lines within specviz\n", - "HTML('')" + "# Make new video to show how to measure lines within specviz" ] }, { @@ -915,8 +853,10 @@ "source": [ "# Alternative method to analyze the 10um line within the notebook. Calculate the Line flux; Line Centroid; Equivalent width\n", "\n", - "line_centroid = centroid(line_spec_consub, SpectralRegion(sw_line*u.um, lw_line*u.um))\n", - "line_flux_val = line_flux(line_spec_consub, SpectralRegion(sw_line*u.um, lw_line*u.um))\n", + "line_centroid = centroid(\n", + " line_spec_consub, SpectralRegion(sw_line*u.um, lw_line*u.um))\n", + "line_flux_val = line_flux(\n", + " line_spec_consub, SpectralRegion(sw_line*u.um, lw_line*u.um))\n", "\n", "equivalent_width_val = equivalent_width(line_spec_norm)\n", "\n", @@ -953,7 +893,7 @@ "plt.figure(figsize=(10, 6))\n", "plt.plot(optdepth_spec.spectral_axis, optdepth_spec.flux)\n", "plt.xlabel(\"Wavelength ({:latex})\".format(spec.spectral_axis.unit))\n", - "plt.ylabel('Tau') \n", + "plt.ylabel('Tau')\n", "plt.tight_layout()\n", "plt.show()\n", "plt.close()" @@ -996,7 +936,8 @@ "**Author:** Olivia Jones, Project Scientist, UK ATC.
\n", "**Updated On:** 2020-08-11
\n", "**Updated On:** 2021-09-06 by B. Sargent, STScI Scientist, Space Telescope Science Institute (added MRS Simulated Data)
\n", - "**Updated On:** 2021-12-12 by O. Fox, STScI Scientist (added blackbody and polynomial fitting within the notebook)
" + "**Updated On:** 2021-12-12 by O. Fox, STScI Scientist (added blackbody and polynomial fitting within the notebook)
\n", + "**Updated On:** 2024-10-29 by C. Pacifici, STScI Data Scientist, adapt to Jdaviz 4.0 (still need to update videos)
" ] }, { @@ -1030,7 +971,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.9.13" + "version": "3.12.7" } }, "nbformat": 4,