From 613c7da5085d9b0824853f7ceb959dae99bc4d13 Mon Sep 17 00:00:00 2001 From: Kyle Conroy Date: Fri, 17 Dec 2021 13:26:10 -0500 Subject: [PATCH 1/2] update notebook and requirements --- .../JWST_Mstar_dataAnalysis_analysis.ipynb | 568 ++++++++++++------ .../JWST_Mstar_dataAnalysis_runpipeline.ipynb | 26 - .../MRS_Mstar_analysis/requirements.txt | 1 - 3 files changed, 394 insertions(+), 201 deletions(-) diff --git a/03_JDatNotebooks/MRS_Mstar_analysis/JWST_Mstar_dataAnalysis_analysis.ipynb b/03_JDatNotebooks/MRS_Mstar_analysis/JWST_Mstar_dataAnalysis_analysis.ipynb index 1e0e055..ef12dbe 100644 --- a/03_JDatNotebooks/MRS_Mstar_analysis/JWST_Mstar_dataAnalysis_analysis.ipynb +++ b/03_JDatNotebooks/MRS_Mstar_analysis/JWST_Mstar_dataAnalysis_analysis.ipynb @@ -20,7 +20,7 @@ "**Data:** Simulated [MIRI MRS](https://jwst-docs.stsci.edu/mid-infrared-instrument/miri-observing-modes/miri-medium-resolution-spectroscopy) spectrum of AGB star.
\n", "**Source of Simulations:** [MIRISim](https://www.stsci.edu/jwst/science-planning/proposal-planning-toolbox/mirisim)
\n", "**Pipeline Version:** [JWST Pipeline](https://jwst-docs.stsci.edu/jwst-data-reduction-pipeline)
\n", - "**Tools:** specutils, jwst, photutils, astropy, aplpy, scipy.
\n", + "**Tools:** specutils, jwst, photutils, 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", "**Note**: This notebook includes MIRI simulated data cubes obtained using MIRISim (https://wiki.miricle.org//bin/view/Public/MIRISim_Public)\n", @@ -34,7 +34,7 @@ "\n", "This is the second notebook, which will perform data analysis using `specutils`. Specifically, it will fit a model photosphere/blackbody to the spectra. Then it will calculate the centroids, line integrated flux and equivalent width for each dust and molecular feature. \n", "\n", - "**This notebook assumes you have run the first notebook.**\n", + "**This notebook utilizes reduced data created in the first notebook (JWST_Mstar_dataAnalysis_runpipeline.ipynb), although you don't have to run that notebook to use this notebook. All data created in notebook number 1 are available for download within this noteobok.**\n", "\n", "## To Do:\n", "- Make function to extract spectra from datacube using a different means of extraction.\n", @@ -73,8 +73,8 @@ "\n", "# Set general plotting options\n", "params = {'legend.fontsize': '18', 'axes.labelsize': '18',\n", - " 'axes.titlesize': '18', 'xtick.labelsize': '18',\n", - " 'ytick.labelsize': '18', 'lines.linewidth': 2, 'axes.linewidth': 2, 'animation.html': 'html5'}\n", + " 'axes.titlesize': '18', 'xtick.labelsize': '18',\n", + " 'ytick.labelsize': '18', 'lines.linewidth': 2, 'axes.linewidth': 2, 'animation.html': 'html5'}\n", "plt.rcParams.update(params)\n", "plt.rcParams.update({'figure.max_open_warning': 0})" ] @@ -113,12 +113,14 @@ "from specutils.analysis import line_flux, centroid, equivalent_width\n", "from specutils.spectra import SpectralRegion\n", "from specutils import SpectrumList\n", - "\n", - "# To make nice plots with WCS axis\n", - "import aplpy\n", + "from jdaviz import Specviz\n", + "from jdaviz import Cubeviz\n", "\n", "# To fit a curve to the data\n", - "from scipy.optimize import curve_fit" + "from scipy.optimize import curve_fit\n", + "\n", + "# Display the video\n", + "from IPython.display import HTML, YouTubeVideo" ] }, { @@ -127,8 +129,18 @@ "metadata": {}, "outputs": [], "source": [ - "import warnings\n", - "warnings.simplefilter('ignore')" + "# 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)" ] }, { @@ -137,16 +149,15 @@ "metadata": {}, "outputs": [], "source": [ - "# Check if Pipeline 3 Reduced data exists and, if not, download it\n", - "import os\n", - "import urllib.request\n", - "\n", - "if os.path.exists(\"combine_dithers_all_exposures_ch1-long_s3d\"):\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", - " " + "def checkKey(dict, key):\n", + " \n", + " if key in dict.keys():\n", + " print(\"Present, \", end=\" \")\n", + " print(\"value =\", dict[key])\n", + " return(True)\n", + " else:\n", + " print(\"Not present\")\n", + " return(False)" ] }, { @@ -155,17 +166,8 @@ "metadata": {}, "outputs": [], "source": [ - "# Unzip Tar Files\n", - "\n", - "import tarfile\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()" + "import warnings\n", + "warnings.simplefilter('ignore')" ] }, { @@ -174,9 +176,29 @@ "metadata": {}, "outputs": [], "source": [ - "# Move Files \n", + "# Check if Pipeline 3 Reduced data exists and, if not, download it\n", + "import os\n", + "import urllib.request\n", "\n", - "os.system('mv reduced/*fits .')" + "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", + "\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", + " os.system('mv reduced/*fits .')" ] }, { @@ -230,7 +252,18 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "#### Developer Note: At this point, the 12 extracted 1D spectra need to get spliced together with a specialty function written for MRS. " + "## Specviz Visualization of the SpectrumList\n", + "\n", + "You can visualize the spectrum list inside a Jupyter notebook using Specviz" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Video 1:\n", + " \n", + "This Specviz Instructional Demo is from STScI's official YouTube channel and provides an introduction to Cubeviz." ] }, { @@ -239,35 +272,45 @@ "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", - "\n", - "# Make a 1D spectrum object\n", - "spec = Spectrum1D(spectral_axis=wav, flux=fl, uncertainty=StdDevUncertainty(efl))\n", - "\n", - "# Apply a 5 pixel boxcar smoothing to the spectrum\n", - "spec_bsmooth = box_smooth(spec, width=5) \n", - "\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", - "plt.xlabel('Wavelength (microns)')\n", - "plt.ylabel(\"Flux ({:latex})\".format(spec.flux.unit))\n", - "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=False)\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", - "\n", - "plt.legend(frameon=False, fontsize='medium')\n", - "plt.tight_layout()\n", - "plt.show()\n", - "plt.close()" + "vid = YouTubeVideo(\"zLyRnfG32Bo\")\n", + "display(vid)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Read in the SpectrumList (12 unique spectra)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "scrolled": false + }, + "outputs": [], + "source": [ + "# Open these spectra up in Specviz\n", + "specviz = Specviz()\n", + "specviz.app" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "#### Load in the spectrum list from above. Note, only the first spectrum in your list is displayed automatically. You will need to turn on the remaining spectra in the \"DATA\" drop-down, then hit the \"Home\" button in the toolbar, and scale our plot accordingly to see the other spectra." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Load in the spectrum list from above. \n", + "specviz.load_spectrum(splist)" ] }, { @@ -283,7 +326,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "## Video1:\n", + "## Video 2:\n", " \n", "This Cubeviz Instructional Demo is from STScI's official YouTube channel and provides an introduction to Cubeviz." ] @@ -294,8 +337,6 @@ "metadata": {}, "outputs": [], "source": [ - "from IPython.display import HTML, YouTubeVideo\n", - "\n", "vid = YouTubeVideo(\"HMSYwiH3Gl4\")\n", "display(vid)" ] @@ -308,8 +349,6 @@ }, "outputs": [], "source": [ - "from jdaviz import Cubeviz\n", - "\n", "cubeviz = Cubeviz()\n", "cubeviz.app" ] @@ -327,7 +366,8 @@ "metadata": {}, "outputs": [], "source": [ - "# Here, we load the data into the Cubeviz app for visual inspection. In this case, we're just looking at a single channel, because Cubeviz can only load a single cube at a time.\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", "cubeviz.load_data(ch1short_cubefile)" @@ -344,7 +384,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "## Video2:\n", + "## Video 3:\n", " \n", "Here is a video that quickly shows how to select a spatial region in Cubeviz." ] @@ -358,7 +398,6 @@ "outputs": [], "source": [ "# Video showing the selection of the star with a circular region of interest\n", - "from IPython.display import HTML\n", "HTML('')" ] }, @@ -372,36 +411,74 @@ "spec_agb = cubeviz.app.get_data_from_viewer('spectrum-viewer', 'Subset 1') # AGB star only" ] }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "#### Developer Note: Since Cubeviz only displays a single cube at a time, you can't extract a full spectrum at the current time. So, you should use the spectrum defined above ('spec')" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Make a 1D spectrum object" + ] + }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ - "# Check to see if user drew a subset in cubeviz. If not, define spec_agb as just the original spectrum.\n", - "if not spec_agb:\n", - " spec_agb = spec\n", - " print(\"defining spec_agb=spec\")\n", - "else: \n", - " print(\"spec_agb defined.\")\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))" ] }, { - "cell_type": "markdown", + "cell_type": "code", + "execution_count": null, "metadata": {}, + "outputs": [], "source": [ - "## Specviz Visualization\n", + "# Apply a 5 pixel boxcar smoothing to the spectrum\n", + "spec_bsmooth = box_smooth(spec, width=5) \n", "\n", - "You can visualize the extracted spectrum inside Specviz" + "# 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", + "plt.xlabel('Wavelength (microns)')\n", + "plt.ylabel(\"Flux ({:latex})\".format(spec.flux.unit))\n", + "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=False)\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", + "\n", + "plt.legend(frameon=False, fontsize='medium')\n", + "plt.tight_layout()\n", + "plt.show()\n", + "plt.close()" ] }, { - "cell_type": "markdown", + "cell_type": "raw", "metadata": {}, "source": [ - "## Video3:\n", - " \n", - "This Specviz Instructional Demo is from STScI's official YouTube channel and provides an introduction to Cubeviz." + "# Check to see if user drew a subset in cubeviz. If not, define spec_agb as just the original spectrum.\n", + "if not spec_agb:\n", + " spec_agb = spec\n", + " print(\"defining spec_agb=spec\")\n", + "else: \n", + " print(\"spec_agb defined.\")\n" ] }, { @@ -410,10 +487,16 @@ "metadata": {}, "outputs": [], "source": [ - "from IPython.display import HTML, YouTubeVideo\n", + "spec_agb = spec" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Visualize for Analysis the Single Spectrum1D Object Created Above from All 12 Individual Spectra\n", "\n", - "vid = YouTubeVideo(\"zLyRnfG32Bo\")\n", - "display(vid)" + "You can visualize the extracted spectrum inside Specviz" ] }, { @@ -425,8 +508,6 @@ "outputs": [], "source": [ "# Open these spectra up in Specviz\n", - "from jdaviz import Specviz\n", - "\n", "specviz = Specviz()\n", "specviz.app" ] @@ -465,7 +546,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "## Video4:\n", + "## Video 4:\n", " \n", "Here is a video that quickly shows how to smooth your spectrum in Specviz." ] @@ -486,41 +567,60 @@ "source": [ "## Data analysis \n", "\n", - "Analyze the spectrum1d object (spec) created from the spectrumlist list above.\n", - "\n", - "#### Developer Note: Add blackbody fitting to the modeling. And, for that matter, more of the astropy models: https://docs.astropy.org/en/stable/modeling/physical_models.html Can't do the next step in Specviz because don't have a blackbody fitter." + "Analyze the spectrum1d object (spec) created from the spectrumlist list above." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ - "### Fit a continuum - find the best-fitting template (stellar photosphere model or blackbody)\n", + "#### Fit a continuum - find the best-fitting template (stellar photosphere model or blackbody)\n", "\n", "For AGB stars with a photosphere component fit a stellar photosphere model or a blackbody to short wavelength end of the spectra.\n", "\n", "#### Developer Note: Would idealy like to fit the photosphere with a set of Phoenix Models.\n", "I think `template_comparison` may be a good function here to work with the Phoenix Models which have been setup to\n", - "interface with `pysynphot`.\n", - "\n", + "interface with `pysynphot`.
\n", "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": { + "scrolled": false + }, + "outputs": [], + "source": [ + "# Video showing how to fit a blackbody \n", + "HTML('')" + ] + }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ - "def blackbody_Fnu(lam, T, A):\n", - " \"\"\" Blackbody as a function of wavelength (um) and temperature (K).\n", - " Function returns the Planck function in f_nu units\n", - " # [Y Jy] = 1.0E+23 * [X erg/cm^2/s/Hz] = 10E+26 [X Watts/m^2/Hz]\n", - " \"\"\"\n", - " from scipy.constants import h, k, c\n", - " lam = 1e-6 * lam # convert to metres\n", - " bb_nu = 2*h*c / (lam**3 * (np.exp(h*c / (lam*k*T)) - 1)) # units of W/m^2/Hz/Steradian ; f_nu units\n", - " return A * bb_nu" + "spectra = specviz.get_spectra()\n", + " \n", + "a = checkKey(spectra, \"BB1\")\n", + "if a is True:\n", + " # Extract Blackbody fit from Specviz\n", + " blackbody = spectra[\"BB1\"]\n", + "if a is False:\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=False)\n", + " blackbody = Spectrum1D.read(fn)" ] }, { @@ -529,12 +629,11 @@ "metadata": {}, "outputs": [], "source": [ - "# Only want to fit to a small wavelength range at the start of the spectra\n", - "phot_fit_region = [5.0, 8.0]# Microns\n", - "\n", - "# Trim the specrum to the region showing a stellar photosphere\n", - "sub_region_phot = SpectralRegion([(phot_fit_region[0], phot_fit_region[1])] * u.micron)\n", - "sub_spectrum_phot = extract_region(spec, sub_region_phot)" + "# Delete any existing output in current directory\n", + "if os.path.exists(\"blackbody.fits\"):\n", + " os.remove(\"blackbody.fits\")\n", + "else:\n", + " print(\"The blackbody.fits file does not exist\")" ] }, { @@ -543,32 +642,40 @@ "metadata": {}, "outputs": [], "source": [ - "# fit BB to the data\n", - "def phot_fn(wa, T1, A):\n", - " return blackbody_Fnu(wa, T1, A) \n", - "\n", - "\n", - "popt, pcov = curve_fit(phot_fn, sub_spectrum_phot.spectral_axis.value,\n", - " sub_spectrum_phot.flux.value, p0=(3000, 10000),\n", - " sigma=sub_spectrum_phot.uncertainty.quantity)\n", - "\n", - "# Get the best fitting parameter value and their 1 sigma errors\n", - "best_t1, best_a1 = popt\n", - "sigma_t1, sigma_a1 = np.sqrt(np.diag(pcov))\n", - "\n", - "ybest = blackbody_Fnu(spec.spectral_axis.value, best_t1, best_a1)\n", - "\n", - "print('Parameters of best-fitting model:')\n", - "print('T1 = %.2f +/- %.2f' % (best_t1, sigma_t1))\n", - "\n", - "degrees_of_freedom = len(sub_spectrum_phot.spectral_axis.value) - 2\n", - "\n", - "resid = (sub_spectrum_phot.flux.value - phot_fn(sub_spectrum_phot.spectral_axis.value, *popt)) \\\n", - " / sub_spectrum_phot.uncertainty.quantity\n", - "\n", - "chisq = np.dot(resid, resid)\n", - "\n", - "print('nchi2 %.2f' % (chisq.value / degrees_of_freedom))" + "# Save if you so desire. Keep commented otherwise.\n", + "# blackbody.write('blackbody.fits')" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Rename blackbody.flux as ybest\n", + "ybest = blackbody.flux" + ] + }, + { + "cell_type": "raw", + "metadata": {}, + "source": [ + "def blackbody_Fnu(lam, T, A):\n", + " \"\"\" Blackbody as a function of wavelength (um) and temperature (K).\n", + " Function returns the Planck function in f_nu units\n", + " # [Y Jy] = 1.0E+23 * [X erg/cm^2/s/Hz] = 10E+26 [X Watts/m^2/Hz]\n", + " \"\"\"\n", + " from scipy.constants import h, k, c\n", + " lam = 1e-6 * lam # convert to metres\n", + " bb_nu = 2*h*c / (lam**3 * (np.exp(h*c / (lam*k*T)) - 1)) # units of W/m^2/Hz/Steradian ; f_nu units\n", + " return A * bb_nu" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "#### Developer Note: At this point, the 12 extracted 1D spectra need to get spliced together with a specialty function written for MRS. " ] }, { @@ -591,8 +698,8 @@ "plt.close()\n", "\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, color='purple', label='Dust spectra')\n", + "plt.figure(figsize=(8, 4))\n", + "plt.plot(spec.spectral_axis, spec.flux.value - 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", @@ -648,7 +755,7 @@ "outputs": [], "source": [ "# Subtract the continuum and plot in a new instance of specviz\n", - "bbsub_spectra = spec - ybest # continuum subtracted spectra - Dust only" + "bbsub_spectra = spec - ybest.value # continuum subtracted spectra - Dust only" ] }, { @@ -657,62 +764,132 @@ "metadata": {}, "outputs": [], "source": [ - "# Fit a spline to the 10 micron feature to isolate it.\n", - "\n", - "bbsub_spectra = spec - ybest # continuum subtracted spectra - Dust only\n", - "\n", + "specviz = Specviz()\n", + "specviz.app" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "specviz.load_spectrum(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('')" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "spectra = specviz.get_spectra()\n", + " \n", + "a = checkKey(spectra, \"PolyFit\")\n", + "if a is True:\n", + " # Extract Blackbody fit from Specviz\n", + " blackbody = spectra[\"PolyFit\"]\n", + "if a is False:\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=False)\n", + " poly = Spectrum1D.read(fn)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Delete any existing output in current directory\n", + "if os.path.exists(\"poly.fits\"):\n", + " os.remove(\"poly.fits\")\n", + "else:\n", + " print(\"The poly.fits file does not exist\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Save if you so desire. Keep commented otherwise.\n", + "# poly.write('poly.fits')" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": { + "scrolled": false + }, + "outputs": [], + "source": [ "# Fit a local continuum between the flux densities at: 8.0 - 8.1 & 14.9 - 15.0 microns\n", "# (i.e. excluding the line itself)\n", "\n", - "sw_region = 8.0 #lam0\n", - "sw_line = 8.1 #lam1\n", - "lw_line = 14.9 #lam2\n", - "lw_region = 15.0 #lam3\n", + "sw_region = 8.0 # lam0\n", + "sw_line = 8.1 # lam1\n", + "lw_line = 14.9 # lam2\n", + "lw_region = 15.0 # lam3\n", "\n", "# Zoom in on the line complex & extract\n", "line_reg_10 = SpectralRegion([(sw_region*u.um, lw_region*u.um)])\n", "line_spec = extract_region(bbsub_spectra, line_reg_10)\n", + "polysub = extract_region(poly, line_reg_10)\n", + "line_y_continuum = polysub.flux\n", "\n", - "# Fit a local continuum - exclude the actual dust feature when doing the fit\n", - "\n", - "lgl_fit = fit_generic_continuum(line_spec, \n", - " exclude_regions = SpectralRegion([(sw_line*u.um, \n", - " lw_line*u.um)])) \n", - "\n", - "# Determine Y values of the line continuum\n", - "line_y_continuum = lgl_fit(line_spec.spectral_axis) \n", - "\n", - "#-----------------------------------------------------------------\n", + "# -----------------------------------------------------------------\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/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", "\n", - "#-----------------------------------------------------------------\n", + "# -----------------------------------------------------------------\n", "# Plot the dust feature & continuum fit to the region\n", "\n", - "plt.figure(figsize = (8, 4))\n", + "plt.figure(figsize=(8, 4))\n", "\n", "plt.plot(line_spec.spectral_axis, line_spec.flux.value,\n", - " label = 'Dust spectra 10 micron region')\n", + " label='Dust spectra 10 micron region')\n", "\n", - "plt.plot(line_spec.spectral_axis, line_y_continuum, label = 'Local continuum')\n", + "plt.plot(line_spec.spectral_axis, line_y_continuum, label='Local continuum')\n", "\n", "plt.xlabel('Wavelength (microns)')\n", "plt.ylabel(\"Flux ({:latex})\".format(spec.flux.unit))\n", "plt.title(\"10$\\mu$m feature plus local continuum\")\n", - "plt.legend(frameon = False, fontsize = 'medium')\n", + "plt.legend(frameon=False, fontsize='medium')\n", "plt.tight_layout()\n", "plt.show()\n", "plt.close()\n", "\n", - "#-----------------------------------------------------------------\n", + "# -----------------------------------------------------------------\n", "# Plot the continuum subtracted 10 micron feature\n", "\n", - "plt.figure(figsize = (8,4))\n", + "plt.figure(figsize=(8, 4))\n", "\n", "plt.plot(line_spec.spectral_axis, line_spec_consub.flux, color='green',\n", - " label = 'continuum subtracted')\n", + " label='continuum subtracted')\n", "\n", "plt.xlabel('Wavelength (microns)')\n", "plt.ylabel(\"Flux ({:latex})\".format(spec.flux.unit))\n", @@ -728,7 +905,49 @@ "metadata": {}, "outputs": [], "source": [ - "# Calculate the Line flux; Line Centroid; Equivalent width\n", + "# Load 10 um feature back into specviz and calculate the Line flux; Line Centroid; Equivalent width\n", + "specviz = Specviz()\n", + "specviz.app" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "specviz.load_spectrum(line_spec_consub, data_label='Continuum Subtraction')\n", + "specviz.load_spectrum(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('')" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "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", @@ -766,7 +985,7 @@ "outputs": [], "source": [ "# Plot the optical depth of the 10 micron region vs wavelength\n", - "plt.figure(figsize=(10,6))\n", + "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", @@ -809,9 +1028,10 @@ }, "source": [ "## About this notebook\n", - "**Author:** Olivia Jones, Project Scientist, UK ATC.\n", - "**Updated On:** 2020-08-11\n", - "**Later Updated On:** 2021-09-06 by B. Sargent, STScI Scientist, Space Telescope Science Institute" + "**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)
" ] }, { @@ -845,7 +1065,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.8.12" + "version": "3.8.10" } }, "nbformat": 4, diff --git a/03_JDatNotebooks/MRS_Mstar_analysis/JWST_Mstar_dataAnalysis_runpipeline.ipynb b/03_JDatNotebooks/MRS_Mstar_analysis/JWST_Mstar_dataAnalysis_runpipeline.ipynb index e6950db..d19b1e7 100644 --- a/03_JDatNotebooks/MRS_Mstar_analysis/JWST_Mstar_dataAnalysis_runpipeline.ipynb +++ b/03_JDatNotebooks/MRS_Mstar_analysis/JWST_Mstar_dataAnalysis_runpipeline.ipynb @@ -36,32 +36,6 @@ "This is the first notebook, which will process the data and automatically detect and extract spectra for the point source. The workflow will use `photutils` to automatically detect sources in the cube to extract the spectrum of the point source. Then it will read in the spectra generated at Stage 3 of the JWST pipeline.\n" ] }, - { - "cell_type": "markdown", - "metadata": {}, - "source": [ - "# TEST" - ] - }, - { - "cell_type": "code", - "execution_count": 1, - "metadata": {}, - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "env: CRDS_PATH=/tmp/crds_cache\n", - "env: CRDS_SERVER_URL=https://jwst-crds.stsci.edu\n" - ] - } - ], - "source": [ - "%env CRDS_PATH=$HOME/crds_cache\n", - "%env CRDS_SERVER_URL=https://jwst-crds.stsci.edu" - ] - }, { "cell_type": "markdown", "metadata": { diff --git a/03_JDatNotebooks/MRS_Mstar_analysis/requirements.txt b/03_JDatNotebooks/MRS_Mstar_analysis/requirements.txt index 9934ab3..c77fe36 100644 --- a/03_JDatNotebooks/MRS_Mstar_analysis/requirements.txt +++ b/03_JDatNotebooks/MRS_Mstar_analysis/requirements.txt @@ -1,7 +1,6 @@ scipy >= 1.7.1 astropy >= 4.3.1 matplotlib == 3.4.3 -aplpy >= 2.0.3 git+https://github.com/radio-astro-tools/radio-beam.git spectral-cube >= 0.5.0 photutils == 1.1.0 From 035d24ad6452ab602f53f1358ce57d59c2ebc7fa Mon Sep 17 00:00:00 2001 From: Kyle Conroy Date: Tue, 21 Dec 2021 11:58:33 -0500 Subject: [PATCH 2/2] typo fix --- .../MRS_Mstar_analysis/JWST_Mstar_dataAnalysis_analysis.ipynb | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/03_JDatNotebooks/MRS_Mstar_analysis/JWST_Mstar_dataAnalysis_analysis.ipynb b/03_JDatNotebooks/MRS_Mstar_analysis/JWST_Mstar_dataAnalysis_analysis.ipynb index ef12dbe..4f24595 100644 --- a/03_JDatNotebooks/MRS_Mstar_analysis/JWST_Mstar_dataAnalysis_analysis.ipynb +++ b/03_JDatNotebooks/MRS_Mstar_analysis/JWST_Mstar_dataAnalysis_analysis.ipynb @@ -9,7 +9,7 @@ }, "source": [ "\n", - "# MIRI MRS Spectroscopy of a Late M Star 2 - Data Anlysis and Visualization" + "# MIRI MRS Spectroscopy of a Late M Star 2 - Data Analysis and Visualization" ] }, {