diff --git a/notebooks/MIRI_LRS_spectral_extraction/miri_lrs_advanced_extraction_part1.ipynb b/notebooks/MIRI_LRS_spectral_extraction/miri_lrs_advanced_extraction_part1.ipynb new file mode 100644 index 000000000..ad83f47bf --- /dev/null +++ b/notebooks/MIRI_LRS_spectral_extraction/miri_lrs_advanced_extraction_part1.ipynb @@ -0,0 +1,791 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "20ce4694", + "metadata": {}, + "source": [ + "# MIRI LRS Slit Spectroscopy: Spectral Extraction using the JWST Pipeline\n", + "\n", + "July 2023\n", + "\n", + "**Use case:** Spectral extraction of slit spectra with the JWST calibration pipeline.
\n", + "**Data:** Publicly available science data
\n", + "**Tools:** jwst, matplotlib, astropy.
\n", + "**Cross-intrument:** NIRSpec, MIRI.
\n", + "**Documentation:** This notebook is part of a STScI's larger [post-pipeline Data Analysis Tools Ecosystem](https://jwst-docs.stsci.edu/jwst-post-pipeline-data-analysis) 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", + "\n", + "### Introduction: Spectral extraction in the JWST calibration pipeline\n", + "\n", + "The JWST calibration pipeline performs spectrac extraction for all spectroscopic data using basic default assumptions that are tuned to produce accurately calibrated spectra for the majority of science cases. This default method is a simple fixed-width boxcar extraction, where the spectrum is summed over a number of pixels along the cross-dispersion axis, over the valid wavelength range. An aperture correction is applied at each pixel along the spectrum to account for flux lost from the finite-width aperture. \n", + "\n", + "The ``extract_1d`` step uses the following inputs for its algorithm:\n", + "- the spectral extraction reference file: this is a json-formatted file, available as a reference file from the [JWST CRDS system](https://jwst-crds.stsci.edu)\n", + "- the bounding box: the ``assign_wcs`` step attaches a bounding box definition to the data, which defines the region over which a valid calibration is available. We will demonstrate below how to visualize this region. \n", + "\n", + "However the ``extract_1d`` step has the capability to perform more complex spectral extractions, requiring some manual editing of parameters and re-running of the pipeline step. \n", + "\n", + "\n", + "### Aims\n", + "\n", + "This notebook will demonstrate how to re-run the spectral extraction step with different settings to illustrate the capabilities of the JWST calibration pipeline. \n", + "\n", + "\n", + "### Assumptions\n", + "\n", + "We will demonstrate the spectral extraction methods on resampled, calibrated spectral images. The basic demo and two examples run on Level 3 data, in which the nod exposures have been combined into a single spectral image. Two examples will use the Level 2b data - one of the nodded exposures. \n", + "\n", + "\n", + "### Test data\n", + "\n", + "The data used in this notebook is an observation of the Type Ia supernova SN2021aefx, observed by Jha et al in PID 2072 (Obs 1). These data were taken with zero exclusive access period, and published in [Kwok et al 2023](https://ui.adsabs.harvard.edu/abs/2023ApJ...944L...3K/abstract). You can retrieve the data from [this Box folder](https://stsci.box.com/s/i2xi18jziu1iawpkom0z2r94kvf9n9kb), and we recommend you place the files in the ``data/`` folder of this repository, or change the directory settings in the notebook prior to running. \n", + "\n", + "You can of course use your own data instead of the demo data. \n", + "\n", + "\n", + "### JWST pipeline version and CRDS context\n", + "\n", + "This notebook was written using the calibration pipeline version 1.10.2. We set the CRDS context explicitly to 1089 to match the current latest version in MAST. If you use different pipeline versions or CRDS context, please read the relevant release notes ([here for pipeline](https://github.com/spacetelescope/jwst), [here for CRDS](https://jwst-crds.stsci.edu)) for possibly relevant changes.\n", + "\n", + "### Contents\n", + "\n", + "1. [The Level 3 data products](#l3data)\n", + "2. [The spectral extraction reference file](#x1dref)\n", + "3. [Example 1: Changing the aperture width](#ex1)\n", + "4. [Example 2: Changing the aperture location](#ex2)\n", + "5. [Example 3: Extraction with background subtraction](#ex3)\n", + "6. [Example 4: Tapered column extraction](#ex4)\n", + "\n" + ] + }, + { + "cell_type": "markdown", + "id": "14ff9543", + "metadata": {}, + "source": [ + "## Import Packages" + ] + }, + { + "cell_type": "markdown", + "id": "e698ce3a-fdaf-4d3b-9109-b980794e94aa", + "metadata": {}, + "source": [ + "- `astropy.io` fits for accessing FITS files\n", + "- `os` for managing system paths\n", + "- `matplotlib` for plotting data\n", + "- `urllib` for downloading data\n", + "- `tarfile` for unpacking data\n", + "- `numpy` for basic array manipulation\n", + "- `jwst` for running JWST pipeline and handling data products\n", + "- `json` for working with json files\n", + "- `crds` for working with JWST reference files" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "3a130dd2", + "metadata": {}, + "outputs": [], + "source": [ + "# Set CRDS variables first\n", + "import os\n", + "\n", + "os.environ['CRDS_CONTEXT'] = 'jwst_1089.pmap'\n", + "os.environ['CRDS_PATH'] = os.environ['HOME']+'/crds_cache'\n", + "os.environ['CRDS_SERVER_URL'] = 'https://jwst-crds.stsci.edu'\n", + "print(f'CRDS cache location: {os.environ[\"CRDS_PATH\"]}')" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "08ddf5f7", + "metadata": {}, + "outputs": [], + "source": [ + "%matplotlib inline\n", + "import urllib.request\n", + "import tarfile\n", + "\n", + "import numpy as np\n", + "import matplotlib.pyplot as plt\n", + "\n", + "import astropy.io.fits as fits\n", + "import astropy.units as u\n", + "from astropy.modeling import models, fitting\n", + "\n", + "import jwst\n", + "from jwst import datamodels\n", + "from jwst.extract_1d import Extract1dStep\n", + "\n", + "from matplotlib.patches import Rectangle\n", + "\n", + "import json\n", + "import crds\n", + "\n", + "print(f'Using JWST calibration pipeline version {jwst.__version__}')" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "305103d5", + "metadata": {}, + "outputs": [], + "source": [ + "data_tar_url = 'https://data.science.stsci.edu/redirect/JWST/jwst-data_analysis_tools/MIRI_LRS_notebook/data.tar.gz'\n", + "\n", + "# Download and unpack data if needed\n", + "if not os.path.exists(\"data.tar.gz\"):\n", + " print(\"Downloading Data\")\n", + " urllib.request.urlretrieve(data_tar_url, 'data.tar.gz')\n", + "if not os.path.exists(\"data/\"):\n", + " print(\"Unpacking Data\")\n", + " with tarfile.open('./data.tar.gz', \"r:gz\") as tar:\n", + " tar.extractall()" + ] + }, + { + "cell_type": "markdown", + "id": "611086f4", + "metadata": {}, + "source": [ + "## 1. The Level 3 Data Products \n", + "\n", + "\n", + "Let's start by plotting the main default Level 3 output products:\n", + "* the ``s2d`` file: this is the 2D image built from the co-added resampled individual nod exposures. \n", + "* the ``x1d`` file: this is the 1-D extracted spectrum, extracted from the Level 3 ``s2d`` file. \n", + "\n", + "The ``s2d`` image shows a bright central trace, flanked by two negative traces. These result from the combination of the nod exposures, each of which also contains a positive and negative trace due to being mutually subtracted for background subtraction. \n", + "\n", + "We restrict the short-wavelength end of the x-axis to 5 micron, as our calibration is very poor below this wavelength. The Level 3 spectrum is extracted from the resampled, dither-combined, calibrated exposure. " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "a8012bfa", + "metadata": {}, + "outputs": [], + "source": [ + "l3_s2d_file = 'data/jw02072-o001_t010_miri_p750l_s2d_1089.fits'\n", + "l3_s2d = datamodels.open(l3_s2d_file)\n", + "fig, ax = plt.subplots(figsize=[2, 8])\n", + "im2d = ax.imshow(l3_s2d.data, origin='lower', aspect='auto', cmap='gist_gray')\n", + "ax.set_xlabel('column')\n", + "ax.set_ylabel('row')\n", + "ax.set_title('SN2021aefx - Level 3 resampled 2D spectral image')\n", + "fig.colorbar(im2d)\n", + "fig.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "c51f421b", + "metadata": {}, + "outputs": [], + "source": [ + "l3_file = 'data/jw02072-o001_t010_miri_p750l_x1d_1089.fits'\n", + "l3_spec = datamodels.open(l3_file)\n", + "\n", + "fig2, ax2 = plt.subplots(figsize=[12, 4])\n", + "ax2.plot(l3_spec.spec[0].spec_table['WAVELENGTH'], l3_spec.spec[0].spec_table['FLUX'])\n", + "ax2.set_xlabel('wavelength (um)')\n", + "ax2.set_ylabel('flux (Jy)')\n", + "ax2.set_title('SN2021aefx - Level 3 spectrum in MAST (pmap 1089)')\n", + "ax2.set_xlim(5., 14.)\n", + "fig2.show()" + ] + }, + { + "cell_type": "markdown", + "id": "3ae110b4", + "metadata": {}, + "source": [ + "## The spectral extraction reference file \n", + "\n", + "The reference file that tells the ``extract_1d`` algorithm what parameters to use is a text file using the `json` format that is available in [CRDS](https://jwst-crds.stsci.edu). The second reference file used in the extraction is the aperture correction; this corrects for the flux lost as a function of wavelength for the extraction aperture size used. You can use the datamodel attributes of the ``x1d`` file to check which extraction reference file was called by the algorithm. \n", + "\n", + "We show below how to examine the file programmatically to see what aperture was used to produce the default Level 3 spectrum shown above. **Note: this json file can easily be opened and edited with a simple text editor**. \n", + "\n", + "Full documentation of the ``extract_1d`` reference file is available [here](https://jwst-pipeline.readthedocs.io/en/latest/jwst/extract_1d/reference_files.html). We recommend you read this page and any links therein carefully to understand how the parameters in the file are applied to the data. " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "f0574ee0-84a0-4fa8-ae54-d7b6ca34a7a7", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "print(f'Spectral extraction reference file used: {l3_spec.meta.ref_file.extract1d.name}')" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "95f20a0a-0f24-4d4b-8480-2c37574ad6e8", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "file_path = 'data/jw02072-o001_t010_miri_p750l_x1d_1089.fits'\n", + "with fits.open(file_path) as hdul:\n", + " header = hdul[0].header\n", + " json_ref_default = crds.getreferences(header)['extract1d']\n", + "\n", + " with open(json_ref_default) as json_ref:\n", + " x1dref_default = json.load(json_ref)\n", + " print('Settings for SLIT data: {}'.format(x1dref_default['apertures'][0]))\n", + " print(' ')\n", + " print('Settings for SLITLESS data: {}'.format(x1dref_default['apertures'][1]))" + ] + }, + { + "cell_type": "markdown", + "id": "ecc11d6b", + "metadata": {}, + "source": [ + "Let's look at what's in this file. \n", + "\n", + "* **id**: identification label, in this case specifying the exposure type the parameters will be applied to.\n", + "* **region_type**: optional, if included must be set to 'target'\n", + "* **disp_axis**: defines the direction of dispersion (1 for x-axis, 2 for y-axis). **For MIRI LRS this should always be set to 2**. \n", + "* **xstart** (int): first pixel in the horizontal direction (x-axis; 0-indexed) \n", + "* **xstop** (int): last pixel in the horizontal direction (x-axis; 0-indexed; limit is **inclusive**)\n", + "* **bkg_order**: \n", + "* **use_source_posn** (True/False): if True, this will use the target coordinates to locate the target in the field, and offset the extraction aperture to this location. **We recommend this is set to False**. \n", + "* **bkg_order**: the polynomial order to be used for background fitting. if the accompanying parameter **bkg_coeff** is not provided, no background fitting will be performed. **For MIRI LRS slit data, default background subtraction is achieved in the Spec2Pipeline, by mutually subtracting nod expsosures**.\n", + "\n", + "As for MIRI LRS the dispersion is in the vertical direction (i.e. `disp_axis` = 2), the extraction aperture width is specified with the coordinates `xstart` and `xstop`. If no coordinates `ystart` and `ystop` are provided, the spectrum will be extracted over the full height of the ``s2d`` cutout region. We can illustrate the default extraction parameters on the Level 3 ``s2d`` file. \n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "703f59cd", + "metadata": {}, + "outputs": [], + "source": [ + "xstart = x1dref_default['apertures'][0]['xstart']\n", + "xstop = x1dref_default['apertures'][0]['xstop']\n", + "ap_height = np.shape(l3_s2d.data)[0]\n", + "ap_width = xstop - xstart + 1\n", + "x1d_rect = Rectangle(xy=(xstart, 0), width=ap_width, height=ap_height, angle=0., edgecolor='red',\n", + " facecolor='None', ls='-', lw=1.5)\n", + "\n", + "fig, ax = plt.subplots(figsize=[2, 8])\n", + "im2d = ax.imshow(l3_s2d.data, origin='lower', aspect='auto', cmap='gist_gray')\n", + "ax.add_patch(x1d_rect)\n", + "ax.set_xlabel('column')\n", + "ax.set_ylabel('row')\n", + "ax.set_title('SN2021aefx - Level 3 resampled 2D spectral image')\n", + "fig.colorbar(im2d)\n", + "fig.show()" + ] + }, + { + "cell_type": "markdown", + "id": "1fd784d9", + "metadata": {}, + "source": [ + "## Example 1: Changing the extraction width \n", + "\n", + "In this example, we demonstrate how to change the extraction width from the default. Instead of 8 pixels, we'll extract 12, keeping the aperture centred on the trace. \n", + "\n", + "We will modify the values in the json files in python in this notebook, but the file can also simply be edited in a text editor. " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "06dc8eb5", + "metadata": {}, + "outputs": [], + "source": [ + "xstart2 = xstart - 2\n", + "xstop2 = xstop + 2\n", + "print('New xstart, xstop values = {0},{1}'.format(xstart2, xstop2))\n", + "\n", + "with open(json_ref_default) as json_ref:\n", + " x1dref_default = json.load(json_ref)\n", + " x1dref_ex1 = x1dref_default.copy()\n", + " x1dref_ex1['apertures'][0]['xstart'] = xstart2\n", + " x1dref_ex1['apertures'][0]['xstop'] = xstop2\n", + "\n", + "with open('x1d_reffile_example1.json', 'w') as jsrefout:\n", + " json.dump(x1dref_ex1, jsrefout, indent=4)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "5bc85413", + "metadata": {}, + "outputs": [], + "source": [ + "ap_width2 = xstop2 - xstart2 + 1\n", + "x1d_rect1 = Rectangle(xy=(xstart, 0), width=ap_width, height=ap_height, angle=0., edgecolor='red',\n", + " facecolor='None', ls='-', lw=1, label='8-px aperture (default)')\n", + "\n", + "x1d_rect2 = Rectangle(xy=(xstart2, 0), width=ap_width2, height=ap_height, angle=0., edgecolor='cyan',\n", + " facecolor='None', ls='-', lw=1, label='12-px aperture')\n", + "\n", + "fig4, ax4 = plt.subplots(figsize=[2, 8])\n", + "im2d = ax4.imshow(l3_s2d.data, origin='lower', aspect='auto', cmap='gist_gray')\n", + "# ax4.add_collection(aps_collection)\n", + "ax4.add_patch(x1d_rect1)\n", + "ax4.add_patch(x1d_rect2)\n", + "\n", + "ax4.set_xlabel('column')\n", + "ax4.set_ylabel('row')\n", + "ax4.set_title('Example 1: Default vs modified extraction aperture')\n", + "ax4.legend(loc=3)\n", + "fig.colorbar(im2d)\n", + "fig.show()" + ] + }, + { + "cell_type": "markdown", + "id": "87efcb9f", + "metadata": {}, + "source": [ + "Next we run the spectral extraction step, using this modified reference file. Note: when a step is run individually the file name suffix is different from when we run the Spec3Pipeline in its entirety. The extracted spectrum will now have ``extract1dstep.fits`` in the filename. The custom parameters we pass to the step call:\n", + "\n", + "* ``output_file``: we provide a custom output filename for this example (including an output filename renders the ``save_results`` parameter obsolete)\n", + "* ``override_extract1d``: here we provide the name of the custom reference file we created above\n", + "\n", + "We will plot the output against the default extracted product. We expect the spectra to be almost identical; differences can be apparent at the longer wavelengths as our path loss correction is less well calibrated in this low SNR region. " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "7304f758", + "metadata": {}, + "outputs": [], + "source": [ + "sp3_ex1 = Extract1dStep.call(l3_s2d, output_dir='data/', \n", + " output_file='lrs_slit_extract_example1', override_extract1d='x1d_reffile_example1.json')" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "91199fd1", + "metadata": {}, + "outputs": [], + "source": [ + "print(sp3_ex1)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "91ebfc64", + "metadata": {}, + "outputs": [], + "source": [ + "fig5, ax5 = plt.subplots(figsize=[12, 4])\n", + "ax5.plot(l3_spec.spec[0].spec_table['WAVELENGTH'], l3_spec.spec[0].spec_table['FLUX'], label='8-px aperture')\n", + "ax5.plot(sp3_ex1.spec[0].spec_table['WAVELENGTH'], sp3_ex1.spec[0].spec_table['FLUX'], label='12-px aperture')\n", + "ax5.set_xlabel('wavelength (um)')\n", + "ax5.set_ylabel('flux (Jy)')\n", + "ax5.set_title('Example 1: Difference aperture sizes')\n", + "ax5.set_xlim(5., 14.)\n", + "ax5.legend()\n", + "fig5.show()" + ] + }, + { + "cell_type": "markdown", + "id": "f28a4f8e", + "metadata": {}, + "source": [ + "## Example 2: Changing aperture location\n", + "\n", + "In this example we will demonstrate spectral extraction at a different location in the slit. A good use case for this is to extract a spectrum from one of the nodded exposures, prior to combination of the nods in the Spec3Pipeline. We will take the ``s2d`` output from the Spec2Pipeline, and extract the spectrum. In the nod 1 exposure we see the spectrum peak is located in column 13 (0-indexed), and we extract a default 8-px fixed-width aperture. " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "8a8cf793", + "metadata": {}, + "outputs": [], + "source": [ + "l2_s2d_file = 'data/jw02072001001_06101_00001_mirimage_s2d.fits'\n", + "l2_s2d = datamodels.open(l2_s2d_file)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "55c81453", + "metadata": {}, + "outputs": [], + "source": [ + "xstart3 = 9\n", + "xstop3 = 17\n", + "\n", + "with open(json_ref_default) as json_ref:\n", + " x1dref_default = json.load(json_ref)\n", + " x1dref_ex2 = x1dref_default.copy()\n", + " x1dref_ex2['apertures'][0]['xstart'] = xstart3\n", + " x1dref_ex2['apertures'][0]['xstop'] = xstop3\n", + "\n", + "with open('x1d_reffile_example2.json', 'w') as jsrefout:\n", + " json.dump(x1dref_ex2, jsrefout, indent=4)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "fe340506", + "metadata": {}, + "outputs": [], + "source": [ + "ap_width3 = xstop3 - xstart3 + 1\n", + "x1d_rect3 = Rectangle(xy=(xstart3, 0), width=ap_width3, height=ap_height, angle=0., edgecolor='red',\n", + " facecolor='None', ls='-', lw=1, label='8-px aperture at nod 1 location')\n", + "\n", + "fig6, ax6 = plt.subplots(figsize=[2, 8])\n", + "im2d = ax6.imshow(l2_s2d.data, origin='lower', aspect='auto', cmap='gist_gray')\n", + "ax6.add_patch(x1d_rect3)\n", + "ax6.set_xlabel('column')\n", + "ax6.set_ylabel('row')\n", + "ax6.set_title('Example 2: Different aperture location')\n", + "ax6.legend(loc=3)\n", + "fig6.colorbar(im2d)\n", + "fig6.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "3b0b287b", + "metadata": {}, + "outputs": [], + "source": [ + "sp2_ex2 = Extract1dStep.call(l2_s2d_file, output_dir='data/', output_file='lrs_slit_extract_example2',\n", + " override_extract1d='x1d_reffile_example2.json')" + ] + }, + { + "cell_type": "markdown", + "id": "32eaa24e", + "metadata": {}, + "source": [ + "Let's again plot the output against the default extracted product. We expect this 1-nod spectrum to be noisier but not significantly different from the combined product. The spectrum may have more bad pixels that manifest as spikes or dips in the spectrum. " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "ce8eccfb", + "metadata": {}, + "outputs": [], + "source": [ + "fig7, ax7 = plt.subplots(figsize=[12, 4])\n", + "ax7.plot(l3_spec.spec[0].spec_table['WAVELENGTH'], l3_spec.spec[0].spec_table['FLUX'], label='default location (nods combined)')\n", + "ax7.plot(sp2_ex2.spec[0].spec_table['WAVELENGTH'], sp2_ex2.spec[0].spec_table['FLUX'], label='nod 1 location (single nod)')\n", + "ax7.set_xlabel('wavelength (um)')\n", + "ax7.set_ylabel('flux (Jy)')\n", + "ax7.set_title('Example 2: Different aperture locations')\n", + "ax7.set_xlim(5., 14.)\n", + "ax7.legend()\n", + "fig7.show()" + ] + }, + { + "cell_type": "markdown", + "id": "39492ab3", + "metadata": {}, + "source": [ + "## Example 3: Extraction with background subtraction\n", + "\n", + "For LRS slit observations, the default background subtraction strategy is performed in the ``background`` step in the Spec2Pipeline; the 2 nodded exposures are mutually subtracted, resulting in each returning a 2D spectral image with a positive and a negative trace, and the background subtracted. \n", + "\n", + "For non-standard cases or slitless LRS data it is however possible to subtract a background as part of the spectral extraction in ``extract_1d``. In the ``extract_1d`` reference file we can pass specific parameters for the background:\n", + "* bkg_coeff (list or list of floats): the regions to be used as background. **This is the main parameter required for background subtraction**\n", + "* bkg_fit (string): the type or method of the background computation. (e.g. None, 'poly', 'mean' or 'median')\n", + "* bkg_order (int): the order of polynomial to fit to background regions. if bkg_fit is not set to 'poly', this parameter will be ignored. \n", + "* smoothing_length (odd int; optional): the width of the boxcar filter that will be used to smooth the background signal in the dispersion direction. This can provide a better quality in case of noisy data. \n", + "\n", + "The 'poly' option for the ``bkg_fit`` parameter will take the value of all pixels in the background region on a given row, and fit a polynomial of order ``bkg_order`` to them. This option can be useful in cases where a gradient is present in the background. \n", + "\n", + "The data we're using here already has the background subtracted so we expect the impact of this to be minimal, but we provide a demonstration using the nod 1, level 2b spectral image. In this example we will calculate the background from 2 4-column windows, setting the ``bkg_fit`` to 'median'. \n", + "\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "b4ea99ea", + "metadata": {}, + "outputs": [], + "source": [ + "rows = [140, 200, 325]\n", + "fig8, ax8 = plt.subplots(figsize=[8, 4])\n", + "ncols = np.shape(l2_s2d.data)[1]\n", + "pltx = np.arange(ncols)\n", + "for rr in rows:\n", + " label = 'row {}'.format(rr)\n", + " ax8.plot(pltx, l2_s2d.data[rr, :], label=label)\n", + "ax8.axvline(x=1, ymin=0, ymax=1, ls='--', lw=1., color='coral', label='background regions')\n", + "ax8.axvline(x=5, ymin=0, ymax=1, ls='--', lw=1., color='coral')\n", + "ax8.axvline(x=39, ymin=0, ymax=1, ls='--', lw=1., color='coral')\n", + "ax8.axvline(x=43, ymin=0, ymax=1, ls='--', lw=1., color='coral')\n", + "ax8.legend()\n", + "fig8.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "1238d6c9", + "metadata": {}, + "outputs": [], + "source": [ + "with open(json_ref_default) as json_ref:\n", + " x1dref_default = json.load(json_ref)\n", + " x1dref_ex3 = x1dref_default.copy()\n", + " x1dref_ex3['apertures'][0]['xstart'] = xstart3\n", + " x1dref_ex3['apertures'][0]['xstop'] = xstop3\n", + " x1dref_ex3['apertures'][0]['bkg_coeff'] = [[0.5], [4.5], [38.5], [43.5]]\n", + " x1dref_ex3['apertures'][0]['bkg_fit'] = 'median'\n", + "\n", + "with open('x1d_reffile_example3.json', 'w') as jsrefout:\n", + " json.dump(x1dref_ex3, jsrefout, indent=4)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "ac0d2746", + "metadata": {}, + "outputs": [], + "source": [ + "sp2_ex3 = Extract1dStep.call(l2_s2d_file, output_dir='data/', output_file='lrs_slit_extract_example3',\n", + " override_extract1d='x1d_reffile_example3.json')" + ] + }, + { + "cell_type": "markdown", + "id": "45b69d9b", + "metadata": {}, + "source": [ + "When the ``extract_1d`` step performs a background subtraction, the background spectrum is part of the output product, so you can check what was subtracted. In the plot below we can see that, as expected, the background for this particular exposure is near-zero (apart from the noisy long-wavelength end), as the subtraction was already performed. " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "7210c9ac", + "metadata": {}, + "outputs": [], + "source": [ + "fig9, ax9 = plt.subplots(nrows=2, ncols=1, figsize=[12, 4])\n", + "# ax9.plot(l3_spec.spec[0].spec_table['WAVELENGTH'], l3_spec.spec[0].spec_table['FLUX'], label='default location (nods combined)')\n", + "ax9[0].plot(sp2_ex2.spec[0].spec_table['WAVELENGTH'], sp2_ex2.spec[0].spec_table['FLUX'], label='nod 1 spectrum - no bkg sub')\n", + "ax9[0].plot(sp2_ex3.spec[0].spec_table['WAVELENGTH'], sp2_ex3.spec[0].spec_table['FLUX'], label='nod 1 spectrum - with bkg sub')\n", + "ax9[1].plot(sp2_ex3.spec[0].spec_table['WAVELENGTH'], sp2_ex3.spec[0].spec_table['BACKGROUND'], label='background')\n", + "ax9[1].set_xlabel('wavelength (um)')\n", + "ax9[0].set_ylabel('flux (Jy)')\n", + "ax9[0].set_title('Example 3: Extraction with background subtraction')\n", + "ax9[0].set_xlim(5., 14.)\n", + "ax9[1].set_xlim(5., 14.)\n", + "ax9[0].legend()\n", + "ax9[1].legend()\n", + "fig9.show()" + ] + }, + { + "cell_type": "markdown", + "id": "e01786eb", + "metadata": {}, + "source": [ + "## Example 4: Tapered column extraction\n", + "\n", + "In this example we will use the JWST calibration pipeline to perform a spectral extraction in a tapered column aperture. The way to specify this in the extraction reference file is to use the ``src_soeff`` parameter instead of the simpler ``xstart``, ``xstop`` settings. The ``src_coeff`` parameter can take polynomial coefficients rather than fixed pixel values. In this example, we will define a tapered column aperture corresponding to 3 * the FWHM of the spatial profile. \n", + "\n", + "Polynomial definitions for the extraction aperture can be specified as a function of pixels or wavelength, which is defined in the ``independent_var`` parameter. \n", + "\n", + "We will use pre-measured FWHM values as a function of **wavelength** to fit a straight line to the FWHM($\\lambda$) profile, and set the extraction parameters according to this fit. The FWHM can of course also be measured directly from the data as well. " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "9e3c2433", + "metadata": {}, + "outputs": [], + "source": [ + "def calc_xap_fit():\n", + " # these are values measured from commissioning data. FWHM is in arcsec.\n", + " lam = [5.0, 7.5, 10.0, 12.0]\n", + " fwhm = [0.29, 0.3, 0.36, 0.42]\n", + "\n", + " # convert from arcsec to pixel using MIRI pixel scaling of 0.11 arcsec/px\n", + " fwhm_px = fwhm / (0.11*u.arcsec/u.pixel)\n", + "\n", + " # we want to extract 3 * fwhm, which means 1.5 * fwhm on either side of the trace\n", + " xap_pix = 1.5 * fwhm_px\n", + "\n", + " # now we want to fit a line to these points\n", + " line_init = models.Linear1D()\n", + " fit = fitting.LinearLSQFitter()\n", + "\n", + " fitted_line = fit(line_init, lam, xap_pix.value)\n", + " print(fitted_line)\n", + "\n", + " fig, ax = plt.subplots(figsize=[8, 4])\n", + " xplt = np.linspace(4.0, 14., num=50)\n", + " ax.plot(lam, xap_pix.value, 'rx', label='1.5 * FWHM(px)')\n", + " ax.plot(xplt, fitted_line(xplt), 'b-', label='best-fit line')\n", + " ax.set_xlabel('wavelength')\n", + " ax.set_ylabel('px')\n", + " ax.legend()\n", + "\n", + " return fitted_line" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "e21fcec5", + "metadata": {}, + "outputs": [], + "source": [ + "poly_pos = calc_xap_fit()\n", + "print(poly_pos.slope, poly_pos.intercept)" + ] + }, + { + "cell_type": "markdown", + "id": "f7dc8785", + "metadata": {}, + "source": [ + "The above polynomial defines the relationship between wavelength and the number of pixels to extract. To ensure that the extractio location is centred on the location of the spectrum, we add to the intercept value the central location of the trace, which is at column 30.5. \n", + "\n", + "In the next cell, we provide these parameters to the ``src_coeff`` parameter in the extraction reference file. **Note: The ``src_coeff`` parameter takes precedence over the ``xstart`` and ``xstop`` parameters if all 3 are present; for clarity we remove the latter from our reference file.**" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "a9c9d403", + "metadata": {}, + "outputs": [], + "source": [ + "trace_cen = 30.5\n", + "\n", + "with open(json_ref_default) as json_ref:\n", + " x1dref_default = json.load(json_ref)\n", + " x1dref_ex4 = x1dref_default.copy()\n", + " x1dref_ex4['apertures'][0]['xstart'] = None\n", + " x1dref_ex4['apertures'][0]['xstop'] = None\n", + " x1dref_ex4['apertures'][0]['independent_var'] = 'wavelength'\n", + " x1dref_ex4['apertures'][0]['src_coeff'] = [[-1*poly_pos.intercept.value + trace_cen, -1*poly_pos.slope.value], [poly_pos.intercept.value + trace_cen, poly_pos.slope.value]]\n", + "\n", + "with open('x1d_reffile_example4.json', 'w') as jsrefout:\n", + " json.dump(x1dref_ex4, jsrefout, indent=4)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "71789e83", + "metadata": {}, + "outputs": [], + "source": [ + "sp3_ex4 = Extract1dStep.call(l3_s2d, output_dir='data/', \n", + " output_file='lrs_slit_extract_example4', override_extract1d='x1d_reffile_example4.json')" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "9d1bc74c", + "metadata": {}, + "outputs": [], + "source": [ + "fig10, ax10 = plt.subplots(figsize=[12, 4])\n", + "ax10.plot(l3_spec.spec[0].spec_table['WAVELENGTH'], l3_spec.spec[0].spec_table['FLUX'], label='default fixed-width aperture')\n", + "ax10.plot(sp3_ex4.spec[0].spec_table['WAVELENGTH'], sp3_ex4.spec[0].spec_table['FLUX'], label='tapered column aperture')\n", + "ax10.set_xlabel('wavelength (um)')\n", + "ax10.set_ylabel('flux (Jy)')\n", + "ax10.set_title('Example 4: Tapered column vs. fixed-width extraction aperture')\n", + "ax10.set_xlim(5., 14.)\n", + "ax10.legend()\n", + "fig10.show()" + ] + }, + { + "cell_type": "markdown", + "id": "378f8393", + "metadata": {}, + "source": [ + "The output spectrum also contains a reference to the number of pixels used for the extraction as a function of wavelength. Let's visualize that too. \n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "78ca0c68", + "metadata": {}, + "outputs": [], + "source": [ + "fig11, ax11 = plt.subplots(figsize=[12, 4])\n", + "ax11.plot(l3_spec.spec[0].spec_table['WAVELENGTH'], l3_spec.spec[0].spec_table['NPIXELS'], label='default fixed-width aperture')\n", + "ax11.plot(sp3_ex4.spec[0].spec_table['WAVELENGTH'], sp3_ex4.spec[0].spec_table['NPIXELS'], label='tapered column aperture')\n", + "ax11.set_xlabel('wavelength (um)')\n", + "ax11.set_ylabel('number of pixels')\n", + "ax11.set_title('Example 4: Number of pixels extracted')\n", + "ax11.set_xlim(5., 14.)\n", + "ax11.legend()\n", + "fig11.show()" + ] + }, + { + "cell_type": "markdown", + "id": "f4501ee7", + "metadata": {}, + "source": [ + "## Summary\n", + "\n", + "We hope this notebook was useful in helping you understand the capabilities of the JWST calibration for spectral extraction. The above examples are not an exhaustive list of all the possibilities: different methods of source and background extraction can be combined for more complex extraction operations. \n", + "\n", + "**If you have any questions, comments, or requests for further demos of these capabilities, please contact the [JWST Helpdesk](http://jwsthelp.stsci.edu/).**" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3 (ipykernel)", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.9.17" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/notebooks/MIRI_LRS_spectral_extraction/requirements.txt b/notebooks/MIRI_LRS_spectral_extraction/requirements.txt index 32260cfff..0c89bc8d8 100644 --- a/notebooks/MIRI_LRS_spectral_extraction/requirements.txt +++ b/notebooks/MIRI_LRS_spectral_extraction/requirements.txt @@ -1,3 +1,5 @@ -matplotlib==3.2.1 -astropy==4.0.1 -git+https://github.com/spacetelescope/jwst@0.16.1 +jdaviz >= 3.6.0 +astropy >= 5.3.1 +jwst >= 1.11.3 +specreduce >= 1.3.0 + diff --git a/notebooks/MIRI_MRS_1987A/87a_part1_download.ipynb b/notebooks/MIRI_MRS_1987A/87a_part1_download.ipynb new file mode 100755 index 000000000..164c07ebe --- /dev/null +++ b/notebooks/MIRI_MRS_1987A/87a_part1_download.ipynb @@ -0,0 +1,767 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "298b8519-c695-4562-a96c-96208063c1c3", + "metadata": {}, + "source": [ + "# MIRI MRS IFU Spectroscopy Part 1: \n", + "# Downloading Data\n", + "\n", + "Aug 2023\n", + "\n", + "**Use case:** Reduce MRS Data With User Defined Master Background Step. This is particularly relevant if you did not obtain a Dedicated Background with your observations. While the pipeline will subtract a sky background derived from an annulus, the underlying background may be prohibitively complicated and the user may wish to measure their own background from elsewhere in the cube.
\n", + "**Data:** Publicly available science data for SN 1987A (Program 1232). For this notebook, we will follow the science workflow outlined by [Jones et al. 2023](https://ui.adsabs.harvard.edu/abs/2023arXiv230706692J/abstract).
\n", + "**Tools:** jwst, jdaviz, matplotlib, astropy.
\n", + "**Cross-intrument:** NIRSpec, MIRI.
\n", + "**Documentation:** This notebook is part of a STScI's larger [post-pipeline Data Analysis Tools Ecosystem](https://jwst-docs.stsci.edu/jwst-post-pipeline-data-analysis) 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", + "\n", + "### Introduction: Spectral extraction in the JWST calibration pipeline\n", + "\n", + "The JWST calibration pipeline performs spectrac extraction for all spectroscopic data using basic default assumptions that are tuned to produce accurately calibrated spectra for the majority of science cases. This default method is a simple fixed-width boxcar extraction, where the spectrum is summed over a number of pixels along the cross-dispersion axis, over the valid wavelength range. An aperture correction is applied at each pixel along the spectrum to account for flux lost from the finite-width aperture. \n", + "\n", + "The ``extract_1d`` step uses the following inputs for its algorithm:\n", + "- the spectral extraction reference file: this is a json-formatted file, available as a reference file from the [JWST CRDS system](https://jwst-crds.stsci.edu)\n", + "- the bounding box: the ``assign_wcs`` step attaches a bounding box definition to the data, which defines the region over which a valid calibration is available. We will demonstrate below how to visualize this region. \n", + "\n", + "However the ``extract_1d`` step has the capability to perform more complex spectral extractions, requiring some manual editing of parameters and re-running of the pipeline step. \n", + "\n", + "\n", + "### Aims\n", + "\n", + "This notebook will demonstrate how to re-run the spectral extraction step with different settings to illustrate the capabilities of the JWST calibration pipeline. \n", + "\n", + "\n", + "### Assumptions\n", + "\n", + "We will demonstrate the spectral extraction methods on resampled, calibrated spectral images. The basic demo and two examples run on Level 3 data, in which the nod exposures have been combined into a single spectral image. Two examples will use the Level 2b data - one of the nodded exposures. \n", + "\n", + "\n", + "### Test data\n", + "\n", + "The data used in this notebook is an observation of the Type Ia supernova SN2021aefx, observed by Jha et al in PID 2072 (Obs 1). These data were taken with zero exclusive access period, and published in [Kwok et al 2023](https://ui.adsabs.harvard.edu/abs/2023ApJ...944L...3K/abstract). You can retrieve the data from [this Box folder](https://stsci.box.com/s/i2xi18jziu1iawpkom0z2r94kvf9n9kb), and we recommend you place the files in the ``data/`` folder of this repository, or change the directory settings in the notebook prior to running. \n", + "\n", + "You can of course use your own data instead of the demo data. \n", + "\n", + "\n", + "### JWST pipeline version and CRDS context\n", + "\n", + "This notebook was written using the calibration pipeline version 1.10.2. We set the CRDS context explicitly to 1089 to match the current latest version in MAST. If you use different pipeline versions or CRDS context, please read the relevant release notes ([here for pipeline](https://github.com/spacetelescope/jwst), [here for CRDS](https://jwst-crds.stsci.edu)) for possibly relevant changes.\n", + "\n", + "### Contents\n", + "\n", + "1. [The Level 3 data products](#l3data)\n", + "2. [The spectral extraction reference file](#x1dref)\n", + "3. [Example 1: Changing the aperture width](#ex1)\n", + "4. [Example 2: Changing the aperture location](#ex2)\n", + "5. [Example 3: Extraction with background subtraction](#ex3)\n", + "6. [Example 4: Tapered column extraction](#ex4)" + ] + }, + { + "cell_type": "markdown", + "id": "96ddb4af", + "metadata": {}, + "source": [ + "## Import Packages" + ] + }, + { + "cell_type": "markdown", + "id": "06416a11", + "metadata": {}, + "source": [ + "- `astropy.io` fits for accessing FITS files\n", + "- `os` for managing system paths\n", + "- `matplotlib` for plotting data\n", + "- `urllib` for downloading data\n", + "- `tarfile` for unpacking data\n", + "- `numpy` for basic array manipulation\n", + "- `jwst` for running JWST pipeline and handling data products\n", + "- `json` for working with json files\n", + "- `crds` for working with JWST reference files" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "36b96e73", + "metadata": {}, + "outputs": [], + "source": [ + "# Set CRDS variables first\n", + "import os\n", + "\n", + "os.environ['CRDS_CONTEXT'] = 'jwst_1089.pmap'\n", + "os.environ['CRDS_PATH'] = os.environ['HOME']+'/crds_cache'\n", + "os.environ['CRDS_SERVER_URL'] = 'https://jwst-crds.stsci.edu'\n", + "print(f'CRDS cache location: {os.environ[\"CRDS_PATH\"]}')" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "21efc012", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "import sys,os, pdb\n", + "# Basic system utilities for interacting with files\n", + "import glob\n", + "import time\n", + "import shutil\n", + "import warnings\n", + "import zipfile\n", + "import urllib.request\n", + "\n", + "# Astropy utilities for opening FITS and ASCII files\n", + "from astropy.io import fits\n", + "from astropy.io import ascii\n", + "from astropy.utils.data import download_file\n", + "from regions import Regions\n", + "from astropy import units as u\n", + "\n", + "from astroquery.mast import Observations\n", + "\n", + "# Astropy utilities for making plots\n", + "from astropy.visualization import (LinearStretch, LogStretch, ImageNormalize, ZScaleInterval)\n", + "\n", + "# Numpy for doing calculations\n", + "import numpy as np\n", + "\n", + "# Matplotlib for making plots\n", + "import matplotlib.pyplot as plt\n", + "from matplotlib import rc\n", + "\n", + "# Import the base JWST package\n", + "import jwst\n", + "\n", + "# JWST pipelines (encompassing many steps)\n", + "from jwst.pipeline import Detector1Pipeline\n", + "from jwst.pipeline import Spec2Pipeline\n", + "from jwst.pipeline import Spec3Pipeline\n", + "\n", + "# JWST pipeline utilities\n", + "from jwst import datamodels # JWST datamodels\n", + "from jwst.associations import asn_from_list as afl # Tools for creating association files\n", + "from jwst.associations.lib.rules_level2_base import DMSLevel2bBase # Definition of a Lvl2 association file\n", + "from jwst.associations.lib.rules_level3_base import DMS_Level3_Base # Definition of a Lvl3 association file\n", + "\n", + "from stcal import dqflags # Utilities for working with the data quality (DQ) arrays\n", + "\n", + "import shutil\n", + "\n", + "# Import packages for multiprocessing. These won't be used on the online demo, but can be\n", + "# very useful for local data processing unless/until they get integrated natively into\n", + "# the cube building code. These need to be imported before anything else.\n", + "\n", + "import multiprocessing\n", + "#multiprocessing.set_start_method('fork')\n", + "from multiprocessing import Pool\n", + "import os\n", + "\n", + "# Set the maximum number of processes to spawn based on available cores\n", + "usage = 'all' # Either 'none' (single thread), 'quarter', 'half', or 'all' available cores\n", + "\n", + "from specutils import Spectrum1D\n", + "from matplotlib.pyplot import cm\n", + "\n", + "from jdaviz import Cubeviz\n", + "\n", + "#shutil.copytree('/astro/armin/data/mshahbandeh/aefx/input_dir/', '/astro/armin/data/mshahbandeh/aefx/input_dir_sc/')\n", + "#shutil.copytree('/astro/armin/data/mshahbandeh/aefx/input_dir/', '/astro/armin/data/mshahbandeh/aefx/input_dir_bkg/')" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "37350f4d-4dbc-445d-b743-f26e7898e73e", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "# Set parameters to be changed here.\n", + "# It should not be necessary to edit cells below this in general unless modifying pipeline processing steps.\n", + "\n", + "import sys,os, pdb\n", + "\n", + "# CRDS context (if overriding)\n", + "#%env CRDS_CONTEXT jwst_0771.pmap\n", + "\n", + "# Point to where the uncalibrated FITS files are from the science observation\n", + "input_dir = './mastDownload/1232/uncal/'\n", + "\n", + "# Point to where you want the output science results to go\n", + "output_dir = './output/87A/'\n", + "\n", + "# Point to where the uncalibrated FITS files are from the background observation\n", + "# If no background observation, leave this blank\n", + "input_bgdir = ' '\n", + "\n", + "# Point to where the output background observations should go\n", + "# If no background observation, leave this blank\n", + "output_bgdir = ' '\n", + "\n", + "# Whether or not to run a given pipeline stage\n", + "# Science and background are processed independently through det1+spec2, and jointly in spec3\n", + "\n", + "# Science processing\n", + "dodet1=True\n", + "dospec2=True\n", + "dospec3=True\n", + "\n", + "# Background processing\n", + "dodet1bg=True\n", + "dospec2bg=True\n", + "\n", + "# If there is no background folder, ensure we don't try to process it\n", + "if (input_bgdir == ''):\n", + " dodet1bg=False\n", + " dospec2bg=False" + ] + }, + { + "cell_type": "raw", + "id": "a1c5254b-6e38-4b43-82c5-44fa67a4f4b0", + "metadata": {}, + "source": [ + "## Point to where the uncalibrated FITS files are from the science observation\n", + "input_dir = '/Users/ofox/data/1860/mast/01860/obsnum03/'\n", + "#\n", + "## Point to where you want the output science results to go\n", + "output_dir = '/Users/ofox/data/1860/output/05ip_3/'" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "7c10e734", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "## Output subdirectories to keep science data products organized\n", + "## Note that the pipeline might complain about this as it is intended to work with everything in a single\n", + "## directory, but it nonetheless works fine for the examples given here.\n", + "det1_dir = os.path.join(output_dir, 'stage1/') # Detector1 pipeline outputs will go here\n", + "#spec2_dir = os.path.join(output_dir, 'stage2/') # Spec2 pipeline outputs will go here\n", + "spec2_dir = os.path.join(output_dir, 'stage2/') # Spec2 pipeline outputs will go here\n", + "spec2_bgdir = ' '\n", + "#spec3_dir = os.path.join(output_dir, 'stage3/') # Spec3 pipeline outputs will go here\n", + "spec3_dir = os.path.join(output_dir, 'stage3/') # Spec3 pipeline outputs will go here\n", + "\n", + "# We need to check that the desired output directories exist, and if not create them\n", + "if not os.path.exists(det1_dir):\n", + " os.makedirs(det1_dir)\n", + "if not os.path.exists(spec2_dir):\n", + " os.makedirs(spec2_dir)\n", + "if not os.path.exists(spec3_dir):\n", + " os.makedirs(spec3_dir)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "ba5341f5-08f7-409b-a764-c3afc160faa2", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "# Output subdirectories to keep background data products organized\n", + "det1_bgdir = os.path.join(output_bgdir, 'stage1/') # Detector1 pipeline outputs will go here\n", + "spec2_bgdir = os.path.join(output_bgdir, 'stage2/') # Spec2 pipeline outputs will go here\n", + "\n", + "# We need to check that the desired output directories exist, and if not create them\n", + "if (output_bgdir != ''):\n", + " if not os.path.exists(det1_bgdir):\n", + " os.makedirs(det1_bgdir)\n", + " if not os.path.exists(spec2_bgdir):\n", + " os.makedirs(spec2_bgdir)" + ] + }, + { + "cell_type": "markdown", + "id": "75a62ae9", + "metadata": {}, + "source": [ + "# 2. Download all MRS data from SN 1987A PID 1232 (Public)" + ] + }, + { + "cell_type": "markdown", + "id": "5bd8798f", + "metadata": {}, + "source": [ + "#### If you want to run the entire MRS pipeline from start to finish, you will need to download nearly 100 GB of data. The vast majority of these data are the Level0 raw ramp (uncal.fits) and Level1 ramp (rate.fits and rateints.fits) files. For our purposes, we encourage you to simply download the Level2 calibrated data (cal.fits), which totals only 3 GB." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "5b4e12dd", + "metadata": {}, + "outputs": [], + "source": [ + "# Let's get a list of all observations associated with this proposal\n", + "obs_list = Observations.query_criteria(proposal_id=1232)\n", + "\n", + "# We can chooose the columns we want to display in our table\n", + "disp_col = ['dataproduct_type','instrument_name','calib_level','obs_id',\n", + " 'target_name','filters','proposal_pi', 'obs_collection']\n", + "obs_list[disp_col].show_in_notebook()" + ] + }, + { + "cell_type": "raw", + "id": "e378fe92", + "metadata": {}, + "source": [ + "# Level 0 uncal.fits downloads. Can skip for this workflow.\n", + "\n", + "mask = (obs_list['instrument_name'] == 'MIRI/IFU')\n", + "data_products = Observations.get_product_list(obs_list[mask])\n", + "\n", + "filtered_prod = Observations.filter_products(data_products, calib_level=[1], productType=\"SCIENCE\")\n", + "\n", + "# Again, we choose columns of interest for convenience\n", + "disp_col = ['obsID','dataproduct_type','productFilename','size','calib_level']\n", + "filtered_prod.show_in_notebook(display_length=10)" + ] + }, + { + "cell_type": "raw", + "id": "ee113760", + "metadata": {}, + "source": [ + "total = sum(filtered_prod['size'])\n", + "print('{:.2f} GB'.format(total/10**9))" + ] + }, + { + "cell_type": "raw", + "id": "1e6988f7", + "metadata": {}, + "source": [ + "# Don't forget to login, if accessing non-public data! You can un-comment the line below:\n", + "# Observations.login()\n", + "\n", + "# You can download all of the products by removing the '[:5]' from the line below:\n", + "manifest = Observations.download_products(filtered_prod)\n", + "print(manifest['Status'])" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "e92c6d23", + "metadata": {}, + "outputs": [], + "source": [ + "# Level 2b cal.fits\n", + "\n", + "mask = (obs_list['instrument_name'] == 'MIRI/IFU')\n", + "data_products = Observations.get_product_list(obs_list[mask])\n", + "\n", + "filtered_prod = Observations.filter_products(data_products, calib_level=[2], productType=\"SCIENCE\", productSubGroupDescription=\"CAL\")\n", + "\n", + "# Again, we choose columns of interest for convenience\n", + "disp_col = ['obsID','dataproduct_type','productFilename','size','calib_level']\n", + "filtered_prod.show_in_notebook(display_length=10)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "d96ecb35", + "metadata": {}, + "outputs": [], + "source": [ + "total = sum(filtered_prod['size'])\n", + "print('{:.2f} GB'.format(total/10**9))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "6bebdd1e", + "metadata": {}, + "outputs": [], + "source": [ + "# Don't forget to login, if accessing non-public data! You can un-comment the line below:\n", + "# Observations.login()\n", + "\n", + "# You can download all of the products by removing the '[:5]' from the line below:\n", + "manifest = Observations.download_products(filtered_prod)\n", + "print(manifest['Status'])" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "8b44f15a", + "metadata": {}, + "outputs": [], + "source": [ + "# Check to see if the input directory exists. If not, create it. Move all _cal.fits files into that directory.\n", + "\n", + "if os.path.exists(input_dir):\n", + " print(input_dir+\" already exists\")\n", + "else:\n", + " print(\"Creating Directory \"+input_dir)\n", + " os.mkdir(input_dir)\n", + " \n", + "if os.path.exists(\"./output\"):\n", + " print(\"./output already exists\")\n", + "else:\n", + " print(\"Creating Directory ./output\")\n", + " os.mkdir(\"./output\")\n", + " \n", + "if os.path.exists(output_dir):\n", + " print(output_dir+\" already exists\")\n", + "else:\n", + " print(\"Creating Directory \"+output_dir)\n", + " os.mkdir(output_dir)\n", + " \n", + "if os.path.exists(det1_dir):\n", + " print(det1_dir+\" already exists\")\n", + "else:\n", + " print(\"Creating Directory \"+det1_dir)\n", + " os.mkdir(det1_dir)\n", + " \n", + "if os.path.exists(spec2_dir):\n", + " print(spec2_dir+\" already exists\")\n", + "else:\n", + " print(\"Creating Directory \"+spec2_dir)\n", + " os.mkdir(spec2_dir)\n", + " \n", + "if os.path.exists(spec3_dir):\n", + " print(spec3_dir+\" already exists\")\n", + "else:\n", + " print(\"Creating Directory \"+spec3_dir)\n", + " os.mkdir(spec3_dir)\n", + " \n", + "print(\"Moving All Uncal Files To Input Directory\")\n", + "for file in glob.glob('./mastDownload/JWST/*/*_uncal.fits'):\n", + " root = file.split('/')\n", + " print(root[-1])\n", + " if os.path.isfile(input_dir+'/'+root[-1]):\n", + " print('Deleting '+input_dir+'/'+root[-1])\n", + " os.remove(input_dir+'/'+root[-1])\n", + " print('Moving '+input_dir+'/'+root[-1])\n", + " shutil.move(file, input_dir)\n", + " \n", + "print(\"Moving All Cal Files To Input Directory\")\n", + "for file in glob.glob('./mastDownload/JWST/*/*_cal.fits'):\n", + " root = file.split('/')\n", + " print(root[-1])\n", + " if os.path.isfile(spec2_dir+'/'+root[-1]):\n", + " print('Deleting '+spec2_dir+'/'+root[-1])\n", + " os.remove(spec2_dir+'/'+root[-1])\n", + " print('Moving '+spec2_dir+'/'+root[-1])\n", + " shutil.move(file, spec2_dir)" + ] + }, + { + "cell_type": "markdown", + "id": "235e74e6-52cb-4c56-b150-2945f3a00230", + "metadata": { + "tags": [] + }, + "source": [ + "# 2. Detector1 Pipeline" + ] + }, + { + "cell_type": "markdown", + "id": "c4591e04", + "metadata": {}, + "source": [ + "#### Not necessary to run the Detector1 stage of the pipeline for this notebook. But here is sample code in case you do." + ] + }, + { + "cell_type": "raw", + "id": "4f9cbdd5", + "metadata": { + "tags": [] + }, + "source": [ + "# If you did want to run the entire pipeline, here are the steps.\n", + "\n", + "# First we'll define a function that will call the detector1 pipeline with our desired set of parameters\n", + "# We won't enumerate the individual steps\n", + "def rundet1(filename, outdir):\n", + " print(filename)\n", + " det1 = Detector1Pipeline() # Instantiate the pipeline\n", + " det1.output_dir = outdir # Specify where the output should go\n", + " \n", + " # The jump and ramp fitting steps can benefit from multi-core processing, but this is off by default\n", + " # Turn them on here if desired by choosing how many cores to use (quarter, half, or all)\n", + " det1.jump.maximum_cores='half'\n", + " det1.ramp_fit.maximum_cores='half'\n", + " \n", + " det1.jump.find_showers=True\n", + " det1.jump.only_use_ints = False\n", + " det1.minimum_sigclip_groups = 75\n", + " \n", + " det1.save_results = True # Save the final resulting _rate.fits files\n", + " det1(filename) # Run the pipeline on an input list of files" + ] + }, + { + "cell_type": "raw", + "id": "6f6a0ab7", + "metadata": { + "tags": [] + }, + "source": [ + "# Now let's look for input files of the form *uncal.fits from the science observation\n", + "sstring = input_dir + 'jw*mirifu*uncal.fits'\n", + "lvl1b_files = sorted(glob.glob(sstring))\n", + "print('Found ' + str(len(lvl1b_files)) + ' science input files to process')" + ] + }, + { + "cell_type": "raw", + "id": "33ecb948", + "metadata": { + "tags": [] + }, + "source": [ + "# Run the pipeline on these input files by a simple loop over our pipeline function\n", + "if dodet1:\n", + " for file in lvl1b_files:\n", + " rundet1(file, det1_dir)\n", + "else:\n", + " print('Skipping Detector1 processing')" + ] + }, + { + "cell_type": "markdown", + "id": "b96ea41b-0b77-4a76-9de2-5f237c551226", + "metadata": {}, + "source": [ + "# 3. Spec2 Pipeline\n" + ] + }, + { + "cell_type": "markdown", + "id": "6e0012b1", + "metadata": {}, + "source": [ + "#### Not necessary to run the Spec2 stage of the pipeline for this notebook. But here is sample code in case you do." + ] + }, + { + "cell_type": "raw", + "id": "032bf7cf", + "metadata": { + "tags": [] + }, + "source": [ + "# Define a function that will call the spec2 pipeline with our desired set of parameters\n", + "# We'll list the individual steps just to make it clear what's running\n", + "def runspec2(filename, outdir, nocubes=False):\n", + " spec2 = Spec2Pipeline()\n", + " spec2.output_dir = outdir\n", + " \n", + " spec2.save_results = True\n", + " spec2(filename)" + ] + }, + { + "cell_type": "raw", + "id": "14d23946", + "metadata": { + "tags": [] + }, + "source": [ + "# Look for uncalibrated science slope files from the Detector1 pipeline\n", + "sstring = det1_dir + 'jw*mirifu*rate.fits'\n", + "ratefiles = sorted(glob.glob(sstring))\n", + "ratefiles=np.array(ratefiles)\n", + "print('Found ' + str(len(ratefiles)) + ' input files to process')" + ] + }, + { + "cell_type": "raw", + "id": "3d6b2d1f", + "metadata": {}, + "source": [ + "if dospec2:\n", + " for file in ratefiles:\n", + " runspec2(file, spec2_dir, nocubes=True)\n", + "else:\n", + " print('Skipping Spec2 processing')" + ] + }, + { + "cell_type": "markdown", + "id": "0e9a59fe", + "metadata": {}, + "source": [ + "# 4. Spec3 Pipeline" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "cecc7e6f", + "metadata": {}, + "outputs": [], + "source": [ + "# Define a useful function to write out a Lvl3 association file from an input list\n", + "# Note that any background exposures have to be of type x1d.\n", + "def writel3asn(scifiles, bgfiles, asnfile, prodname):\n", + " # Define the basic association of science files\n", + " asn = afl.asn_from_list(scifiles, rule=DMS_Level3_Base, product_name=prodname)\n", + " \n", + " # Add background files to the association\n", + " nbg=len(bgfiles)\n", + " for ii in range(0,nbg):\n", + " asn['products'][0]['members'].append({'expname': bgfiles[ii], 'exptype': 'background'})\n", + " \n", + " # Write the association to a json file\n", + " _, serialized = asn.dump()\n", + " with open(asnfile, 'w') as outfile:\n", + " outfile.write(serialized)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "a65d12e6", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "# Find and sort all of the input files\n", + "\n", + "# Science Files need the cal.fits files\n", + "sstring = spec2_dir + 'jw*mirifu*cal.fits'\n", + "calfiles = np.array(sorted(glob.glob(sstring)))\n", + "\n", + "# Background Files need the x1d.fits files\n", + "sstring = spec2_bgdir + 'jw*mirifu*x1d.fits'\n", + "bgfiles = np.array(sorted(glob.glob(sstring)))\n", + "\n", + "print('Found ' + str(len(calfiles)) + ' science files to process')\n", + "print('Found ' + str(len(bgfiles)) + ' background files to process')" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "4d4e0e2e", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "# Make an association file that includes all of the different exposures\n", + "#asnfile=os.path.join(output_dir, 'spec2_l3asn.json')\n", + "asnfile='spec2_l3asn.json'\n", + "dospec3 = 1.\n", + "if dospec3:\n", + " writel3asn(calfiles, bgfiles, asnfile, 'Level3')" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "284cfcf2", + "metadata": {}, + "outputs": [], + "source": [ + "asnfile" + ] + }, + { + "cell_type": "markdown", + "id": "1965c445", + "metadata": {}, + "source": [ + "#### Running spec3 in 'multi' output mode to create a single cube with all sub-bands stitched together." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "ec50de18", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "# Define a function that will call the spec3 pipeline with our desired set of parameters\n", + "# This is designed to run on an association file\n", + "def runspec3(filename):\n", + " crds_config = Spec3Pipeline.get_config_from_reference(filename)\n", + " spec3 = Spec3Pipeline.from_config_section(crds_config)\n", + "\n", + " spec3.output_dir = spec3_dir\n", + " spec3.save_results = True\n", + "\n", + " # Cube building configuration options\n", + " spec3.cube_build.output_type = 'multi' # 'band', 'channel', or 'multi' type cube output\n", + "\n", + " # Overrides for whether or not certain steps should be skipped\n", + " spec3.master_background.skip = True\n", + " spec3.subtract_background = False\n", + " spec3.extract_1d.subtract_background=False\n", + " spec3(filename)\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "fd96a3a5", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "spec3 = 1.\n", + "if dospec3:\n", + " runspec3(asnfile)\n", + "else:\n", + " print('Skipping Spec3 processing')" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "403777ba", + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3 (ipykernel)", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.10.12" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/notebooks/MIRI_MRS_1987A/87a_part2_background_sub.ipynb b/notebooks/MIRI_MRS_1987A/87a_part2_background_sub.ipynb new file mode 100755 index 000000000..eb7408eb7 --- /dev/null +++ b/notebooks/MIRI_MRS_1987A/87a_part2_background_sub.ipynb @@ -0,0 +1,756 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "298b8519-c695-4562-a96c-96208063c1c3", + "metadata": {}, + "source": [ + "# MIRI MRS IFU Spectroscopy Part 2: \n", + "# Defining and Extracting a Background Spectrum\n", + "\n", + "Aug 2023\n", + "\n", + "**Use case:** Reduce MRS Data With User Defined Master Background Step. This is particularly relevant if you did not obtain a Dedicated Background with your observations. While the pipeline will subtract a sky background derived from an annulus, the underlying background may be prohibitively complicated and the user may wish to measure their own background from elsewhere in the cube.
\n", + "**Data:** Publicly available science data for SN 1987A (Program 1232). For this notebook, we will follow the science workflow outlined by [Jones et al. 2023](https://ui.adsabs.harvard.edu/abs/2023arXiv230706692J/abstract).
\n", + "**Tools:** jwst, jdaviz, matplotlib, astropy.
\n", + "**Cross-intrument:** NIRSpec, MIRI.
\n", + "**Documentation:** This notebook is part of a STScI's larger [post-pipeline Data Analysis Tools Ecosystem](https://jwst-docs.stsci.edu/jwst-post-pipeline-data-analysis) 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", + "\n", + "### Introduction: Spectral extraction in the JWST calibration pipeline\n", + "\n", + "The JWST calibration pipeline performs spectrac extraction for all spectroscopic data using basic default assumptions that are tuned to produce accurately calibrated spectra for the majority of science cases. This default method is a simple fixed-width boxcar extraction, where the spectrum is summed over a number of pixels along the cross-dispersion axis, over the valid wavelength range. An aperture correction is applied at each pixel along the spectrum to account for flux lost from the finite-width aperture. \n", + "\n", + "The ``extract_1d`` step uses the following inputs for its algorithm:\n", + "- the spectral extraction reference file: this is a json-formatted file, available as a reference file from the [JWST CRDS system](https://jwst-crds.stsci.edu)\n", + "- the bounding box: the ``assign_wcs`` step attaches a bounding box definition to the data, which defines the region over which a valid calibration is available. We will demonstrate below how to visualize this region. \n", + "\n", + "However the ``extract_1d`` step has the capability to perform more complex spectral extractions, requiring some manual editing of parameters and re-running of the pipeline step. \n", + "\n", + "\n", + "### Aims\n", + "\n", + "This notebook will demonstrate how to re-run the spectral extraction step with different settings to illustrate the capabilities of the JWST calibration pipeline. \n", + "\n", + "\n", + "### Assumptions\n", + "\n", + "We will demonstrate the spectral extraction methods on resampled, calibrated spectral images. The basic demo and two examples run on Level 3 data, in which the nod exposures have been combined into a single spectral image. Two examples will use the Level 2b data - one of the nodded exposures. \n", + "\n", + "\n", + "### Test data\n", + "\n", + "The data used in this notebook is an observation of the Type Ia supernova SN2021aefx, observed by Jha et al in PID 2072 (Obs 1). These data were taken with zero exclusive access period, and published in [Kwok et al 2023](https://ui.adsabs.harvard.edu/abs/2023ApJ...944L...3K/abstract). You can retrieve the data from [this Box folder](https://stsci.box.com/s/i2xi18jziu1iawpkom0z2r94kvf9n9kb), and we recommend you place the files in the ``data/`` folder of this repository, or change the directory settings in the notebook prior to running. \n", + "\n", + "You can of course use your own data instead of the demo data. \n", + "\n", + "\n", + "### JWST pipeline version and CRDS context\n", + "\n", + "This notebook was written using the calibration pipeline version 1.10.2. We set the CRDS context explicitly to 1089 to match the current latest version in MAST. If you use different pipeline versions or CRDS context, please read the relevant release notes ([here for pipeline](https://github.com/spacetelescope/jwst), [here for CRDS](https://jwst-crds.stsci.edu)) for possibly relevant changes.\n", + "\n", + "### Contents\n", + "\n", + "1. [The Level 3 data products](#l3data)\n", + "2. [The spectral extraction reference file](#x1dref)\n", + "3. [Example 1: Changing the aperture width](#ex1)\n", + "4. [Example 2: Changing the aperture location](#ex2)\n", + "5. [Example 3: Extraction with background subtraction](#ex3)\n", + "6. [Example 4: Tapered column extraction](#ex4)" + ] + }, + { + "cell_type": "markdown", + "id": "e58c7be6", + "metadata": {}, + "source": [ + "## Import Packages" + ] + }, + { + "cell_type": "markdown", + "id": "5b8c9d77", + "metadata": {}, + "source": [ + "- `astropy.io` fits for accessing FITS files\n", + "- `os` for managing system paths\n", + "- `matplotlib` for plotting data\n", + "- `urllib` for downloading data\n", + "- `tarfile` for unpacking data\n", + "- `numpy` for basic array manipulation\n", + "- `jwst` for running JWST pipeline and handling data products\n", + "- `json` for working with json files\n", + "- `crds` for working with JWST reference files" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "0d1f42bc", + "metadata": {}, + "outputs": [], + "source": [ + "# Set CRDS variables first\n", + "import os\n", + "\n", + "os.environ['CRDS_CONTEXT'] = 'jwst_1089.pmap'\n", + "os.environ['CRDS_PATH'] = os.environ['HOME']+'/crds_cache'\n", + "os.environ['CRDS_SERVER_URL'] = 'https://jwst-crds.stsci.edu'\n", + "print(f'CRDS cache location: {os.environ[\"CRDS_PATH\"]}')" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "21efc012", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "import sys,os, pdb\n", + "# Basic system utilities for interacting with files\n", + "import glob\n", + "import time\n", + "import shutil\n", + "import warnings\n", + "import zipfile\n", + "import urllib.request\n", + "import requests\n", + "\n", + "# Astropy utilities for opening FITS and ASCII files\n", + "from astropy.io import fits\n", + "from astropy.io import ascii\n", + "from astropy.utils.data import download_file\n", + "from regions import Regions\n", + "from astropy import units as u\n", + "\n", + "from astroquery.mast import Observations\n", + "\n", + "# Astropy utilities for making plots\n", + "from astropy.visualization import (LinearStretch, LogStretch, ImageNormalize, ZScaleInterval)\n", + "\n", + "# Numpy for doing calculations\n", + "import numpy as np\n", + "\n", + "# Matplotlib for making plots\n", + "import matplotlib.pyplot as plt\n", + "from matplotlib import rc\n", + "\n", + "# Import the base JWST package\n", + "import jwst\n", + "\n", + "# JWST pipelines (encompassing many steps)\n", + "from jwst.pipeline import Detector1Pipeline\n", + "from jwst.pipeline import Spec2Pipeline\n", + "from jwst.pipeline import Spec3Pipeline\n", + "\n", + "# JWST pipeline utilities\n", + "from jwst import datamodels # JWST datamodels\n", + "from jwst.associations import asn_from_list as afl # Tools for creating association files\n", + "from jwst.associations.lib.rules_level2_base import DMSLevel2bBase # Definition of a Lvl2 association file\n", + "from jwst.associations.lib.rules_level3_base import DMS_Level3_Base # Definition of a Lvl3 association file\n", + "from jwst.datamodels import SpecModel, MultiSpecModel, IFUCubeModel\n", + "\n", + "\n", + "from stcal import dqflags # Utilities for working with the data quality (DQ) arrays\n", + "\n", + "import shutil\n", + "\n", + "# Import packages for multiprocessing. These won't be used on the online demo, but can be\n", + "# very useful for local data processing unless/until they get integrated natively into\n", + "# the cube building code. These need to be imported before anything else.\n", + "\n", + "import multiprocessing\n", + "#multiprocessing.set_start_method('fork')\n", + "from multiprocessing import Pool\n", + "import os\n", + "\n", + "# Set the maximum number of processes to spawn based on available cores\n", + "usage = 'all' # Either 'none' (single thread), 'quarter', 'half', or 'all' available cores\n", + "\n", + "from specutils import Spectrum1D\n", + "from matplotlib.pyplot import cm\n", + "\n", + "from jdaviz import Cubeviz\n", + "\n", + "# Display the video\n", + "from IPython.display import HTML, YouTubeVideo\n", + "\n", + "#shutil.copytree('/astro/armin/data/mshahbandeh/aefx/input_dir/', '/astro/armin/data/mshahbandeh/aefx/input_dir_sc/')\n", + "#shutil.copytree('/astro/armin/data/mshahbandeh/aefx/input_dir/', '/astro/armin/data/mshahbandeh/aefx/input_dir_bkg/')" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "37350f4d-4dbc-445d-b743-f26e7898e73e", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "# Set parameters to be changed here.\n", + "# It should not be necessary to edit cells below this in general unless modifying pipeline processing steps.\n", + "\n", + "import sys,os, pdb\n", + "\n", + "# CRDS context (if overriding)\n", + "#%env CRDS_CONTEXT jwst_0771.pmap\n", + "\n", + "# Point to where the uncalibrated FITS files are from the science observation\n", + "input_dir = './mastDownload/1232/uncal/'\n", + "\n", + "# Point to where you want the output science results to go\n", + "output_dir = './output/87A/'\n", + "\n", + "# Point to where the uncalibrated FITS files are from the background observation\n", + "# If no background observation, leave this blank\n", + "input_bgdir = ' '\n", + "\n", + "# Point to where the output background observations should go\n", + "# If no background observation, leave this blank\n", + "output_bgdir = ' '\n", + "\n", + "# Whether or not to run a given pipeline stage\n", + "# Science and background are processed independently through det1+spec2, and jointly in spec3\n", + "\n", + "# Science processing\n", + "dodet1=True\n", + "dospec2=True\n", + "dospec3=True\n", + "\n", + "# Background processing\n", + "dodet1bg=True\n", + "dospec2bg=True\n", + "\n", + "# If there is no background folder, ensure we don't try to process it\n", + "if (input_bgdir == ''):\n", + " dodet1bg=False\n", + " dospec2bg=False" + ] + }, + { + "cell_type": "raw", + "id": "a1c5254b-6e38-4b43-82c5-44fa67a4f4b0", + "metadata": {}, + "source": [ + "## Point to where the uncalibrated FITS files are from the science observation\n", + "input_dir = '/Users/ofox/data/1860/mast/01860/obsnum03/'\n", + "#\n", + "## Point to where you want the output science results to go\n", + "output_dir = '/Users/ofox/data/1860/output/05ip_3/'" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "7c10e734", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "## Output subdirectories to keep science data products organized\n", + "## Note that the pipeline might complain about this as it is intended to work with everything in a single\n", + "## directory, but it nonetheless works fine for the examples given here.\n", + "det1_dir = os.path.join(output_dir, 'stage1/') # Detector1 pipeline outputs will go here\n", + "#spec2_dir = os.path.join(output_dir, 'stage2/') # Spec2 pipeline outputs will go here\n", + "spec2_dir = os.path.join(output_dir, 'stage2/') # Spec2 pipeline outputs will go here\n", + "spec2_bgdir = ' '\n", + "#spec3_dir = os.path.join(output_dir, 'stage3/') # Spec3 pipeline outputs will go here\n", + "spec3_dir = os.path.join(output_dir, 'stage3/') # Spec3 pipeline outputs will go here\n", + "\n", + "# We need to check that the desired output directories exist, and if not create them\n", + "if not os.path.exists(det1_dir):\n", + " os.makedirs(det1_dir)\n", + "if not os.path.exists(spec2_dir):\n", + " os.makedirs(spec2_dir)\n", + "if not os.path.exists(spec3_dir):\n", + " os.makedirs(spec3_dir)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "ba5341f5-08f7-409b-a764-c3afc160faa2", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "# Output subdirectories to keep background data products organized\n", + "det1_bgdir = os.path.join(output_bgdir, 'stage1/') # Detector1 pipeline outputs will go here\n", + "spec2_bgdir = os.path.join(output_bgdir, 'stage2/') # Spec2 pipeline outputs will go here\n", + "\n", + "# We need to check that the desired output directories exist, and if not create them\n", + "if (output_bgdir != ''):\n", + " if not os.path.exists(det1_bgdir):\n", + " os.makedirs(det1_bgdir)\n", + " if not os.path.exists(spec2_bgdir):\n", + " os.makedirs(spec2_bgdir)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "86205cd2", + "metadata": {}, + "outputs": [], + "source": [ + "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)" + ] + }, + { + "cell_type": "markdown", + "id": "7def969a", + "metadata": {}, + "source": [ + "# 2. Use Cubeviz to make mask" + ] + }, + { + "cell_type": "markdown", + "id": "9e4eef1a", + "metadata": {}, + "source": [ + "#### This step show how to interactively define a region to be used for extracting a background. If you skip this step, you can continue to run the notebook further in Step 3." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "b7aa348f", + "metadata": {}, + "outputs": [], + "source": [ + "# Video showing how to define an annulus background around SN 1987A using the cells below\n", + "\n", + "HTML('')" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "e89d0586", + "metadata": {}, + "outputs": [], + "source": [ + "cubefile = \"/astro/armin/data/ofox/1232/output/87A/stage3/Level3_ch1-2-3-4-shortmediumlong_s3d.fits\"" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "a7224f0b", + "metadata": { + "scrolled": false + }, + "outputs": [], + "source": [ + "from jdaviz import Cubeviz\n", + "cubeviz = Cubeviz()\n", + "cubeviz.load_data(cubefile, data_label='SN1987A')\n", + "cubeviz.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "6b9dc7eb", + "metadata": {}, + "outputs": [], + "source": [ + "# Get the collapsed cube, ideally from Cubeviz, but otherwise download from pre-defined files.\n", + "\n", + "collapse_cube = cubeviz.app.get_data_from_viewer(\"uncert-viewer\") # AGN Center Model Cube\n", + "if checkKey(collapse_cube,\"collapsed\") is True:\n", + " collapse_cube = cubeviz.app.get_data_from_viewer(\"uncert-viewer\",\"collapsed\") # AGN Center Model Cube\n", + " collapse_cube.write(spec3_dir+\"collapsed_cube.fits\",overwrite='True')\n", + "else:\n", + " print(\"No Collapsed Cube in Cubeviz.\")\n", + " if os.path.isfile(spec3_dir+'/'+\"collapsed_cube.fits\"):\n", + " print('File exists. Deleting '+spec3_dir+'/'+\"collapsed_cube.fits\")\n", + " os.remove(spec3_dir+'/'+\"collapsed_cube.fits\")\n", + " print(\"Downloading to \"+spec3_dir)\n", + " url = 'https://data.science.stsci.edu/redirect/JWST/jwst-data_analysis_tools/MRS_1987A/collapsed_cube.fits'\n", + " urllib.request.urlretrieve(url, './collapsed_cube.fits')\n", + " shutil.move(\"collapsed_cube.fits\",spec3_dir)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "b145607c", + "metadata": {}, + "outputs": [], + "source": [ + "# Get the ellipse region, ideally from Cubeviz, but otherwise download from pre-defined files.\n", + "\n", + "regions = cubeviz.get_interactive_regions()\n", + "if checkKey(regions,\"Subset 1\") is True:\n", + " regions['Subset 1'].write('my_elipse.reg', overwrite=True)\n", + "else:\n", + " print(\"No Background Region From Cubeviz.\")\n", + " if os.path.isfile(\"./my_elipse.reg\"):\n", + " print('File exists. Deleting ./my_elipse.reg')\n", + " os.remove(\"./my_elipse.reg\")\n", + " print(\"Downloading...\")\n", + " fname = \"./my_elipse.reg\"\n", + " url = 'https://data.science.stsci.edu/redirect/JWST/jwst-data_analysis_tools/MRS_1987A/my_elipse.txt'\n", + " urllib.request.urlretrieve(url, fname)\n", + " #fn = x = download_file(url, cache=False)\n", + " #reg = Regions.read(fn, format='ds9')[0]\n" + ] + }, + { + "cell_type": "markdown", + "id": "24590618", + "metadata": {}, + "source": [ + "# 3. Apply the mask to the weights extension of the data cube" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "dabb87c3", + "metadata": {}, + "outputs": [], + "source": [ + "# Read in data cube as a JWST data model\n", + "spec_model_cube = IFUCubeModel()\n", + "spec_model_cube.read(cubefile)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "a7334184", + "metadata": {}, + "outputs": [], + "source": [ + "# Print the source and aperture type being used in the header of the file (this can be Extended or Point). \n", + "# For SN 1987A, we have an EXTENDED source\n", + "\n", + "spec_model_cube.find_fits_keyword('SRCTYPE')\n", + "spec_model_cube.find_fits_keyword('SRCTYAPT')\n", + "\n", + "print(spec_model_cube.meta.target.source_type)\n", + "print(spec_model_cube.meta.target.source_type_apt)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "730eb1d2", + "metadata": {}, + "outputs": [], + "source": [ + "# If necessary, you can change your cube header as necessary. We don't need to change anything in this case.\n", + "# But you might want to if you have a point source, yet want to extract a user specified background spectrum.\n", + "# The pipeline extracts EXTENDED and POINT sources differently.\n", + "\n", + "spec_model_cube.meta.target.source_type = 'EXTENDED'\n", + "spec_model_cube.meta.target.source_type_apt = 'EXTENDED'" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "21adcb34", + "metadata": {}, + "outputs": [], + "source": [ + "# Read in previously extracted region\n", + "\n", + "reg = Regions.read('./my_elipse.reg', format='ds9')[0]\n", + "print(reg)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "4ce2d007", + "metadata": {}, + "outputs": [], + "source": [ + "# Create a weight map using the region mask\n", + "\n", + "tmp_wgts = spec_model_cube.weightmap[:]\n", + "mask = reg.to_mask('exact')\n", + "x1 = int(reg.center.x-reg.width/2.)\n", + "x2 = x1+mask.shape[1]\n", + "y1 = int(reg.center.y-reg.height/2.)\n", + "y2 = y1+mask.shape[0]\n", + "\n", + "### Note above, the region shape is slightly different than the mask shape that gets generated. \n", + "### This hack gets all the arrays to be the same size.\n", + "\n", + "# Start by setting all pixels to 1.\n", + "mask2d = tmp_wgts[13,:,:]\n", + "mask2d[mask2d>0] = 1.\n", + "\n", + "# Because we want an inverse array, we can't just use the mask, we have to subtract the mask (which is 1's) from the original mask2d above (make sense?) \n", + "mask2d[y1:y2,x1:x2] = mask2d[y1:y2,x1:x2]-mask\n", + "\n", + "# Take into account weird rounding errors\n", + "mask2d[mask2d<0.1] = 0." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "c51239ca", + "metadata": {}, + "outputs": [], + "source": [ + "# Visualize the 2D Mask\n", + "\n", + "from astropy.nddata import CCDData\n", + "from astropy.visualization import simple_norm\n", + "ccd = mask2d\n", + "norm = simple_norm(ccd, 'sqrt', min_cut=0, max_cut=0.5) \n", + "color = 'rgbmkrgbmk'\n", + "\n", + "xceni = [36, 44]\n", + "yceni = [66, 58]\n", + "\n", + "fig = plt.figure(figsize=(10, 10))\n", + "ax = fig.add_subplot()\n", + "counter = 0\n", + "plt.title(\"Masked Data\")\n", + "plt.imshow(ccd, norm=norm, origin=\"lower\") \n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "d800d33b", + "metadata": {}, + "outputs": [], + "source": [ + "# Create a 3D weightmap from the 2D map. Mask all NAN values, too.\n", + "\n", + "mask3d = np.broadcast_to(mask2d, spec_model_cube.weightmap.shape)\n", + "mask3d.flags.writeable = True\n", + "mask3d[np.isnan(spec_model_cube.data)] = 0\n", + "#mask_sci_cube = np.ma.masked_array(spec_model_cube.weightmap, mask=mask3d.astype(bool))\n", + "tmp_wgt_cube = np.swapaxes(mask3d,0,1)\n", + "tmp_wgt_cube = np.swapaxes(tmp_wgt_cube,1,2)\n", + "plotcube = Spectrum1D(tmp_wgt_cube*u.dimensionless_unscaled)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "e7b9a1bf", + "metadata": {}, + "outputs": [], + "source": [ + "# Visualize the 3D Cube in Cubeviz\n", + "\n", + "cubeviz2 = Cubeviz()\n", + "cubeviz2.load_data(plotcube, data_label='SN1987A MASK')\n", + "cubeviz2.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "03177ab1", + "metadata": {}, + "outputs": [], + "source": [ + "# Define the weightmap in the original cube as the new 3D Mask \n", + "# This will tell the pipeline which spaxels to use for extraction\n", + "\n", + "spec_model_cube.weightmap = mask3d\n", + "spec_model_cube.save(spec3_dir+'87A_skycube.fits',overwrite=True)" + ] + }, + { + "cell_type": "markdown", + "id": "6b381d59", + "metadata": {}, + "source": [ + "# 4. Extract Background Spectrum using Extract1dStep" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "17914f4a", + "metadata": {}, + "outputs": [], + "source": [ + "# Set 1D extraction parameters\n", + "\n", + "def runex(filename, outdir, outputfile):\n", + " ex1d = jwst.extract_1d.Extract1dStep()\n", + " ex1d.output_dir = outdir\n", + " ex1d.save_results = True\n", + " ex1d.subtract_background = False\n", + " ex1d.output_file = outputfile\n", + " ex1d(filename)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "63c46fc4", + "metadata": { + "scrolled": true + }, + "outputs": [], + "source": [ + "# We will extract a 1D spectrum from the cube created above with a weightmap defined by the region mask\n", + "# This extraction will create an average of all the spaxels in each frame that are not masked\n", + "# We will use this to be our master background to subtract from the entire cube\n", + "\n", + "cubefile_p1 = spec3_dir+'87A_skycube.fits'\n", + "outputfile = spec3_dir+'87A_bg'\n", + "runex(cubefile_p1,spec3_dir,outputfile=outputfile)" + ] + }, + { + "cell_type": "markdown", + "id": "f7f50b83", + "metadata": {}, + "source": [ + "### A single background spectrum is necessary in this workflow. But in some workflows, you may wish to work with more than one background region. This is further illustrated in Notebook 3 of this series on SN 1987A." + ] + }, + { + "cell_type": "markdown", + "id": "b36656ad", + "metadata": {}, + "source": [ + "# 4. Rerun Stage 3 With Master Background Turned On" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "41c7b505", + "metadata": {}, + "outputs": [], + "source": [ + "# Define a function that will call the spec3 pipeline with our desired set of parameters\n", + "# This is designed to run on an association file\n", + "def runspec3(filename):\n", + " \n", + " crds_config = Spec3Pipeline.get_config_from_reference(filename)\n", + " spec3 = Spec3Pipeline.from_config_section(crds_config)\n", + " spec3.output_dir = spec3_dir\n", + " spec3.save_results = True\n", + " spec3.cube_build.output_file = '87A_bg_sub' # Custom output name\n", + " spec3.cube_build.output_type = 'multi' # 'band', 'channel', or 'multi' type cube output\n", + " spec3.outlier_detection.threshold_percent = 98.5 # optimized threshold number\n", + " spec3.master_background.user_background=spec3_dir+'87A_bg_extract1dstep.fits' # Master Background Extracted Above\n", + " spec3.master_background.force_subtract=True\n", + " \n", + " spec3(filename)\n" + ] + }, + { + "cell_type": "markdown", + "id": "f450bfca", + "metadata": {}, + "source": [ + "### Developer Note: Right now, this association file can only be created manually. " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "bc38c7a3", + "metadata": {}, + "outputs": [], + "source": [ + "# Define the association file used for the background subtraction.\n", + "# Download the background subtract json file.\n", + "\n", + "asnfile_bg_sub = 'spec2_l3asn_bg_sub.json'\n", + "if os.path.isfile(asnfile_bg_sub):\n", + " print('File exists. Deleting '+asnfile_bg_sub)\n", + " os.remove(asnfile_bg_sub)\n", + "print(\"Downloading \"+asnfile_bg_sub)\n", + "url = 'https://data.science.stsci.edu/redirect/JWST/jwst-data_analysis_tools/MRS_1987A/spec2_l3asn_bg_sub.json'\n", + "urllib.request.urlretrieve(url, asnfile_bg_sub)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "da5d3a36", + "metadata": { + "scrolled": true + }, + "outputs": [], + "source": [ + "spec3 = 1.\n", + "if dospec3:\n", + " runspec3(asnfile_bg_sub)\n", + "else:\n", + " print('Skipping Spec3 processing')" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "ad227282", + "metadata": {}, + "outputs": [], + "source": [ + "datacube = spec3_dir+'87A_bg_sub_ch1-2-3-4-shortmediumlong_s3d.fits'\n", + "datacube" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "ad317338", + "metadata": {}, + "outputs": [], + "source": [ + "# Visualize Background Subtracted Cube\n", + "# Note, this is not a perfect subtraction. As we will see in the next notebook, we oversubtract the background.\n", + "# This is likely caused by poorly chose spaxels. A more careful selection of background regions should result in a flatter background post-subtraction.\n", + "\n", + "cubeviz3 = Cubeviz()\n", + "cubeviz3.load_data(datacube, data_label='SN1987A MASK')\n", + "cubeviz3.show()" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3 (ipykernel)", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.10.12" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/notebooks/MIRI_MRS_1987A/87a_part3_extract.ipynb b/notebooks/MIRI_MRS_1987A/87a_part3_extract.ipynb new file mode 100755 index 000000000..ee4524f73 --- /dev/null +++ b/notebooks/MIRI_MRS_1987A/87a_part3_extract.ipynb @@ -0,0 +1,983 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "298b8519-c695-4562-a96c-96208063c1c3", + "metadata": {}, + "source": [ + "# MIRI MRS IFU Spectroscopy Part 3: \n", + "# Extracting a Spectrum of the SN 1987A Annulus\n", + "\n", + "Aug 2023\n", + "\n", + "**Use case:** Reduce MRS Data With User Defined Master Background Step. This is particularly relevant if you did not obtain a Dedicated Background with your observations. While the pipeline will subtract a sky background derived from an annulus, the underlying background may be prohibitively complicated and the user may wish to measure their own background from elsewhere in the cube.
\n", + "**Data:** Publicly available science data for SN 1987A (Program 1232). For this notebook, we will follow the science workflow outlined by [Jones et al. 2023](https://ui.adsabs.harvard.edu/abs/2023arXiv230706692J/abstract).
\n", + "**Tools:** jwst, jdaviz, matplotlib, astropy.
\n", + "**Cross-intrument:** NIRSpec, MIRI.
\n", + "**Documentation:** This notebook is part of a STScI's larger [post-pipeline Data Analysis Tools Ecosystem](https://jwst-docs.stsci.edu/jwst-post-pipeline-data-analysis) 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", + "\n", + "### Introduction: Spectral extraction in the JWST calibration pipeline\n", + "\n", + "The JWST calibration pipeline performs spectrac extraction for all spectroscopic data using basic default assumptions that are tuned to produce accurately calibrated spectra for the majority of science cases. This default method is a simple fixed-width boxcar extraction, where the spectrum is summed over a number of pixels along the cross-dispersion axis, over the valid wavelength range. An aperture correction is applied at each pixel along the spectrum to account for flux lost from the finite-width aperture. \n", + "\n", + "The ``extract_1d`` step uses the following inputs for its algorithm:\n", + "- the spectral extraction reference file: this is a json-formatted file, available as a reference file from the [JWST CRDS system](https://jwst-crds.stsci.edu)\n", + "- the bounding box: the ``assign_wcs`` step attaches a bounding box definition to the data, which defines the region over which a valid calibration is available. We will demonstrate below how to visualize this region. \n", + "\n", + "However the ``extract_1d`` step has the capability to perform more complex spectral extractions, requiring some manual editing of parameters and re-running of the pipeline step. \n", + "\n", + "\n", + "### Aims\n", + "\n", + "This notebook will demonstrate how to re-run the spectral extraction step with different settings to illustrate the capabilities of the JWST calibration pipeline. \n", + "\n", + "\n", + "### Assumptions\n", + "\n", + "We will demonstrate the spectral extraction methods on resampled, calibrated spectral images. The basic demo and two examples run on Level 3 data, in which the nod exposures have been combined into a single spectral image. Two examples will use the Level 2b data - one of the nodded exposures. \n", + "\n", + "\n", + "### Test data\n", + "\n", + "The data used in this notebook is an observation of the Type Ia supernova SN2021aefx, observed by Jha et al in PID 2072 (Obs 1). These data were taken with zero exclusive access period, and published in [Kwok et al 2023](https://ui.adsabs.harvard.edu/abs/2023ApJ...944L...3K/abstract). You can retrieve the data from [this Box folder](https://stsci.box.com/s/i2xi18jziu1iawpkom0z2r94kvf9n9kb), and we recommend you place the files in the ``data/`` folder of this repository, or change the directory settings in the notebook prior to running. \n", + "\n", + "You can of course use your own data instead of the demo data. \n", + "\n", + "\n", + "### JWST pipeline version and CRDS context\n", + "\n", + "This notebook was written using the calibration pipeline version 1.10.2. We set the CRDS context explicitly to 1089 to match the current latest version in MAST. If you use different pipeline versions or CRDS context, please read the relevant release notes ([here for pipeline](https://github.com/spacetelescope/jwst), [here for CRDS](https://jwst-crds.stsci.edu)) for possibly relevant changes.\n", + "\n", + "### Contents\n", + "\n", + "1. [The Level 3 data products](#l3data)\n", + "2. [The spectral extraction reference file](#x1dref)\n", + "3. [Example 1: Changing the aperture width](#ex1)\n", + "4. [Example 2: Changing the aperture location](#ex2)\n", + "5. [Example 3: Extraction with background subtraction](#ex3)\n", + "6. [Example 4: Tapered column extraction](#ex4)" + ] + }, + { + "cell_type": "markdown", + "id": "2e9ff09d", + "metadata": {}, + "source": [ + "## Import Packages" + ] + }, + { + "cell_type": "markdown", + "id": "a3122a08", + "metadata": {}, + "source": [ + "- `astropy.io` fits for accessing FITS files\n", + "- `os` for managing system paths\n", + "- `matplotlib` for plotting data\n", + "- `urllib` for downloading data\n", + "- `tarfile` for unpacking data\n", + "- `numpy` for basic array manipulation\n", + "- `jwst` for running JWST pipeline and handling data products\n", + "- `json` for working with json files\n", + "- `crds` for working with JWST reference files" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "0cb67bbb", + "metadata": {}, + "outputs": [], + "source": [ + "# Set CRDS variables first\n", + "import os\n", + "\n", + "os.environ['CRDS_CONTEXT'] = 'jwst_1089.pmap'\n", + "os.environ['CRDS_PATH'] = os.environ['HOME']+'/crds_cache'\n", + "os.environ['CRDS_SERVER_URL'] = 'https://jwst-crds.stsci.edu'\n", + "print(f'CRDS cache location: {os.environ[\"CRDS_PATH\"]}')" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "21efc012", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "import sys,os, pdb\n", + "# Basic system utilities for interacting with files\n", + "import glob\n", + "import time\n", + "import shutil\n", + "import warnings\n", + "import zipfile\n", + "import urllib.request\n", + "import requests\n", + "\n", + "# Astropy utilities for opening FITS and ASCII files\n", + "from astropy.io import fits\n", + "from astropy.io import ascii\n", + "from astropy.utils.data import download_file\n", + "from regions import Regions\n", + "from astropy import units as u\n", + "\n", + "from astroquery.mast import Observations\n", + "\n", + "# Astropy utilities for making plots\n", + "from astropy.visualization import (LinearStretch, LogStretch, ImageNormalize, ZScaleInterval)\n", + "\n", + "# Numpy for doing calculations\n", + "import numpy as np\n", + "\n", + "# Matplotlib for making plots\n", + "import matplotlib.pyplot as plt\n", + "from matplotlib import rc\n", + "\n", + "# Import the base JWST package\n", + "import jwst\n", + "\n", + "# JWST pipelines (encompassing many steps)\n", + "from jwst.pipeline import Detector1Pipeline\n", + "from jwst.pipeline import Spec2Pipeline\n", + "from jwst.pipeline import Spec3Pipeline\n", + "\n", + "# JWST pipeline utilities\n", + "from jwst import datamodels # JWST datamodels\n", + "from jwst.associations import asn_from_list as afl # Tools for creating association files\n", + "from jwst.associations.lib.rules_level2_base import DMSLevel2bBase # Definition of a Lvl2 association file\n", + "from jwst.associations.lib.rules_level3_base import DMS_Level3_Base # Definition of a Lvl3 association file\n", + "from jwst.datamodels import SpecModel, MultiSpecModel, IFUCubeModel\n", + "\n", + "\n", + "from stcal import dqflags # Utilities for working with the data quality (DQ) arrays\n", + "\n", + "import shutil\n", + "\n", + "# Import packages for multiprocessing. These won't be used on the online demo, but can be\n", + "# very useful for local data processing unless/until they get integrated natively into\n", + "# the cube building code. These need to be imported before anything else.\n", + "\n", + "import multiprocessing\n", + "#multiprocessing.set_start_method('fork')\n", + "from multiprocessing import Pool\n", + "import os\n", + "\n", + "# Set the maximum number of processes to spawn based on available cores\n", + "usage = 'all' # Either 'none' (single thread), 'quarter', 'half', or 'all' available cores\n", + "\n", + "from specutils import Spectrum1D\n", + "from matplotlib.pyplot import cm\n", + "\n", + "from jdaviz import Cubeviz\n", + "\n", + "# Display the video\n", + "from IPython.display import HTML, YouTubeVideo\n", + "\n", + "#shutil.copytree('/astro/armin/data/mshahbandeh/aefx/input_dir/', '/astro/armin/data/mshahbandeh/aefx/input_dir_sc/')\n", + "#shutil.copytree('/astro/armin/data/mshahbandeh/aefx/input_dir/', '/astro/armin/data/mshahbandeh/aefx/input_dir_bkg/')" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "37350f4d-4dbc-445d-b743-f26e7898e73e", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "# Set parameters to be changed here.\n", + "# It should not be necessary to edit cells below this in general unless modifying pipeline processing steps.\n", + "\n", + "import sys,os, pdb\n", + "\n", + "# CRDS context (if overriding)\n", + "#%env CRDS_CONTEXT jwst_0771.pmap\n", + "\n", + "# Point to where the uncalibrated FITS files are from the science observation\n", + "input_dir = './mastDownload/1232/uncal/'\n", + "\n", + "# Point to where you want the output science results to go\n", + "output_dir = './output/87A/'\n", + "\n", + "# Point to where the uncalibrated FITS files are from the background observation\n", + "# If no background observation, leave this blank\n", + "input_bgdir = ' '\n", + "\n", + "# Point to where the output background observations should go\n", + "# If no background observation, leave this blank\n", + "output_bgdir = ' '\n", + "\n", + "# Whether or not to run a given pipeline stage\n", + "# Science and background are processed independently through det1+spec2, and jointly in spec3\n", + "\n", + "# Science processing\n", + "dodet1=True\n", + "dospec2=True\n", + "dospec3=True\n", + "\n", + "# Background processing\n", + "dodet1bg=True\n", + "dospec2bg=True\n", + "\n", + "# If there is no background folder, ensure we don't try to process it\n", + "if (input_bgdir == ''):\n", + " dodet1bg=False\n", + " dospec2bg=False" + ] + }, + { + "cell_type": "raw", + "id": "a1c5254b-6e38-4b43-82c5-44fa67a4f4b0", + "metadata": {}, + "source": [ + "## Point to where the uncalibrated FITS files are from the science observation\n", + "input_dir = '/Users/ofox/data/1860/mast/01860/obsnum03/'\n", + "#\n", + "## Point to where you want the output science results to go\n", + "output_dir = '/Users/ofox/data/1860/output/05ip_3/'" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "7c10e734", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "## Output subdirectories to keep science data products organized\n", + "## Note that the pipeline might complain about this as it is intended to work with everything in a single\n", + "## directory, but it nonetheless works fine for the examples given here.\n", + "det1_dir = os.path.join(output_dir, 'stage1/') # Detector1 pipeline outputs will go here\n", + "#spec2_dir = os.path.join(output_dir, 'stage2/') # Spec2 pipeline outputs will go here\n", + "spec2_dir = os.path.join(output_dir, 'stage2/') # Spec2 pipeline outputs will go here\n", + "spec2_bgdir = ' '\n", + "#spec3_dir = os.path.join(output_dir, 'stage3/') # Spec3 pipeline outputs will go here\n", + "spec3_dir = os.path.join(output_dir, 'stage3/') # Spec3 pipeline outputs will go here\n", + "\n", + "# We need to check that the desired output directories exist, and if not create them\n", + "if not os.path.exists(det1_dir):\n", + " os.makedirs(det1_dir)\n", + "if not os.path.exists(spec2_dir):\n", + " os.makedirs(spec2_dir)\n", + "if not os.path.exists(spec3_dir):\n", + " os.makedirs(spec3_dir)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "ba5341f5-08f7-409b-a764-c3afc160faa2", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "# Output subdirectories to keep background data products organized\n", + "det1_bgdir = os.path.join(output_bgdir, 'stage1/') # Detector1 pipeline outputs will go here\n", + "spec2_bgdir = os.path.join(output_bgdir, 'stage2/') # Spec2 pipeline outputs will go here\n", + "\n", + "# We need to check that the desired output directories exist, and if not create them\n", + "if (output_bgdir != ''):\n", + " if not os.path.exists(det1_bgdir):\n", + " os.makedirs(det1_bgdir)\n", + " if not os.path.exists(spec2_bgdir):\n", + " os.makedirs(spec2_bgdir)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "abd7f5dd", + "metadata": {}, + "outputs": [], + "source": [ + "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)" + ] + }, + { + "cell_type": "markdown", + "id": "a490cbb6", + "metadata": {}, + "source": [ + "# 2. Extract Same Background Region Post-Background Subtraction for Comparison" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "bbecca47", + "metadata": {}, + "outputs": [], + "source": [ + "# Read in data cube as a JWST data model\n", + "\n", + "cubefile = spec3_dir+'87A_bg_sub_ch1-2-3-4-shortmediumlong_s3d.fits'\n", + "\n", + "spec_model_cube = IFUCubeModel()\n", + "spec_model_cube.read(cubefile)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "9004d455", + "metadata": {}, + "outputs": [], + "source": [ + "# Print the source and aperture type being used in the header of the file (this can be Extended or Point). \n", + "# For SN 1987A, we have an EXTENDED source\n", + "\n", + "spec_model_cube.find_fits_keyword('SRCTYPE')\n", + "spec_model_cube.find_fits_keyword('SRCTYAPT')\n", + "\n", + "print(spec_model_cube.meta.target.source_type)\n", + "print(spec_model_cube.meta.target.source_type_apt)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "1d24b002", + "metadata": {}, + "outputs": [], + "source": [ + "# If necessary, you can change your cube header as necessary. We don't need to change anything in this case.\n", + "# But you might want to if you have a point source, yet want to extract a user specified background spectrum.\n", + "# The pipeline extracts EXTENDED and POINT sources differently.\n", + "\n", + "spec_model_cube.meta.target.source_type = 'EXTENDED'\n", + "spec_model_cube.meta.target.source_type_apt = 'EXTENDED'" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "6a440cc5", + "metadata": {}, + "outputs": [], + "source": [ + "# Read in previously extracted region (From Notebook 2)\n", + "\n", + "reg = Regions.read('my_elipse.reg', format='ds9')[0]\n", + "print(reg)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "2779d8ef", + "metadata": {}, + "outputs": [], + "source": [ + "# Create a weight map using the region mask\n", + "\n", + "tmp_wgts = spec_model_cube.weightmap[:]\n", + "mask = reg.to_mask('exact')\n", + "x1 = int(reg.center.x-reg.width/2.)\n", + "x2 = x1+mask.shape[1]\n", + "y1 = int(reg.center.y-reg.height/2.)\n", + "y2 = y1+mask.shape[0]\n", + "\n", + "### Note above, the region shape is slightly different than the mask shape that gets generated. \n", + "### This hack gets all the arrays to be the same size.\n", + "\n", + "# Start by setting all pixels to 1.\n", + "mask2d = tmp_wgts[13,:,:]\n", + "mask2d[mask2d>0] = 1.\n", + "\n", + "# Because we want an inverse array, we can't just use the mask, we have to subtract the mask (which is 1's) from the original mask2d above (make sense?) \n", + "mask2d[y1:y2,x1:x2] = mask2d[y1:y2,x1:x2]-mask\n", + "\n", + "# Take into account weird rounding errors\n", + "mask2d[mask2d<0.1] = 0." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "f2f5ff36", + "metadata": {}, + "outputs": [], + "source": [ + "# Visualize the 2D Mask (Same As Before)\n", + "\n", + "from astropy.nddata import CCDData\n", + "from astropy.visualization import simple_norm\n", + "ccd = mask2d\n", + "norm = simple_norm(ccd, 'sqrt', min_cut=0, max_cut=0.5) \n", + "color = 'rgbmkrgbmk'\n", + "\n", + "xceni = [36, 44]\n", + "yceni = [66, 58]\n", + "\n", + "fig = plt.figure(figsize=(10, 10))\n", + "ax = fig.add_subplot()\n", + "counter = 0\n", + "plt.title(\"Masked Data\")\n", + "plt.imshow(ccd, norm=norm, origin=\"lower\") \n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "24282902", + "metadata": {}, + "outputs": [], + "source": [ + "# Create a 3D weightmap from the 2D map. Mask all NAN values, too.\n", + "\n", + "mask3d = np.broadcast_to(mask2d, spec_model_cube.weightmap.shape)\n", + "mask3d.flags.writeable = True\n", + "mask3d[np.isnan(spec_model_cube.data)] = 0\n", + "#mask_sci_cube = np.ma.masked_array(spec_model_cube.weightmap, mask=mask3d.astype(bool))\n", + "tmp_wgt_cube = np.swapaxes(mask3d,0,1)\n", + "tmp_wgt_cube = np.swapaxes(tmp_wgt_cube,1,2)\n", + "plotcube = Spectrum1D(tmp_wgt_cube*u.dimensionless_unscaled)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "1d3996a0", + "metadata": {}, + "outputs": [], + "source": [ + "# Visualize the 3D Cube in Cubeviz\n", + "\n", + "cubeviz = Cubeviz()\n", + "cubeviz.load_data(plotcube, data_label='SN1987A MASK')\n", + "cubeviz.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "27fe1bac", + "metadata": {}, + "outputs": [], + "source": [ + "# Define the weightmap in the background subtracted cube as the new 3D Mask \n", + "# This will tell the pipeline which spaxels to use for extraction\n", + "\n", + "spec_model_cube.weightmap = mask3d\n", + "spec_model_cube.save(spec3_dir+'87A_bg_sub_skymask.fits',overwrite=True)" + ] + }, + { + "cell_type": "markdown", + "id": "967b38cc", + "metadata": {}, + "source": [ + "# 3. Extract Background Spectrum (again) using Extract1dStep" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "34a08066", + "metadata": {}, + "outputs": [], + "source": [ + "# Set 1D extraction parameters\n", + "\n", + "def runex(filename, outdir, outputfile):\n", + " ex1d = jwst.extract_1d.Extract1dStep()\n", + " ex1d.output_dir = outdir\n", + " ex1d.save_results = True\n", + " ex1d.subtract_background = False\n", + " ex1d.output_file = outputfile\n", + " ex1d(filename)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "565396df", + "metadata": { + "scrolled": true + }, + "outputs": [], + "source": [ + "# We will extract a 1D spectrum from the cube created above with a weightmap defined by the region mask\n", + "# This extraction will create an average of all the spaxels in each frame that are not masked\n", + "# We will use this to be our master background to subtract from the entire cube\n", + "\n", + "cubefile_p1 = spec3_dir+'87A_bg_sub_skymask.fits'\n", + "\n", + "outputfile = spec3_dir+'87A_bg_sub_skymask'\n", + "runex(cubefile_p1,spec3_dir,outputfile=outputfile)" + ] + }, + { + "cell_type": "markdown", + "id": "80466d2f", + "metadata": {}, + "source": [ + "# 4. Create Annulus to Define Extraction Region for SN Equitorial Ring" + ] + }, + { + "cell_type": "markdown", + "id": "b8519a0e", + "metadata": {}, + "source": [ + "#### This step show how to interactively define a region to be used for extracting a equatorial ring. If you skip this step, you can continue to run the notebook further in Step 5." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "357667d3", + "metadata": {}, + "outputs": [], + "source": [ + "# Video showing how to define an annulus background around SN 1987A using the cells below\n", + "\n", + "HTML('')" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "3fcd3de3", + "metadata": {}, + "outputs": [], + "source": [ + "# Open Background Subtracted Cube Again\n", + "\n", + "plotcube = spec3_dir+'87A_bg_sub_ch1-2-3-4-shortmediumlong_s3d.fits'\n", + "cubeviz2 = Cubeviz()\n", + "cubeviz2.load_data(plotcube, data_label='SN1987A BG SUB')\n", + "cubeviz2.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "42207155", + "metadata": {}, + "outputs": [], + "source": [ + "# Save Inner and Outer Annuli\n", + "# Get the ellipse region, ideally from Cubeviz, but otherwise download from pre-defined files.\n", + "\n", + "regions = cubeviz2.get_interactive_regions()\n", + "\n", + "if checkKey(regions,\"Subset 1\") is True:\n", + " regions['Subset 1'].write('87a_extract_outer.reg', overwrite=True)\n", + "else:\n", + " print(\"No Outer Annulus Region From Cubeviz.\")\n", + " if os.path.isfile(\"./87a_extract_outer.reg\"):\n", + " print('File exists. Deleting ./87a_extract_outer.reg')\n", + " os.remove(\"./87a_extract_outer.reg\")\n", + " print(\"Downloading...87a_extract_outer.txt\")\n", + " fname = \"./87a_extract_outer.reg\"\n", + " url = 'https://data.science.stsci.edu/redirect/JWST/jwst-data_analysis_tools/MRS_1987A/87a_extract_outer.txt'\n", + " urllib.request.urlretrieve(url, fname)\n", + " \n", + "if checkKey(regions,\"Subset 2\") is True:\n", + " regions['Subset 2'].write('87a_extract_inner.reg', overwrite=True)\n", + "else:\n", + " print(\"No Inner Annulus Region From Cubeviz.\")\n", + " if os.path.isfile(\"./87a_extract_inner.reg\"):\n", + " print('File exists. Deleting ./87a_extract_inner.reg')\n", + " os.remove(\"./87a_extract_inner.reg\")\n", + " print(\"Downloading...87a_extract_inner.txt\")\n", + " fname = \"./87a_extract_inner.reg\"\n", + " url = 'https://data.science.stsci.edu/redirect/JWST/jwst-data_analysis_tools/MRS_1987A/87a_extract_inner.txt'\n", + " urllib.request.urlretrieve(url, fname)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "fdffa541", + "metadata": {}, + "outputs": [], + "source": [ + "# Read in data cube as a JWST data model\n", + "cubefile = spec3_dir+'87A_bg_sub_ch1-2-3-4-shortmediumlong_s3d.fits'\n", + "\n", + "spec_model_cube = IFUCubeModel()\n", + "spec_model_cube.read(cubefile)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "3f968269", + "metadata": {}, + "outputs": [], + "source": [ + "# Print the source and aperture type being used in the header of the file (this can be Extended or Point). \n", + "# For SN 1987A, we have an EXTENDED source\n", + "\n", + "spec_model_cube.find_fits_keyword('SRCTYPE')\n", + "spec_model_cube.find_fits_keyword('SRCTYAPT')\n", + "\n", + "print(spec_model_cube.meta.target.source_type)\n", + "print(spec_model_cube.meta.target.source_type_apt)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "65565988", + "metadata": {}, + "outputs": [], + "source": [ + "# If necessary, you can change your cube header as necessary. We don't need to change anything in this case.\n", + "# But you might want to if you have a point source, yet want to extract a user specified background spectrum.\n", + "# The pipeline extracts EXTENDED and POINT sources differently.\n", + "\n", + "spec_model_cube.meta.target.source_type = 'EXTENDED'\n", + "spec_model_cube.meta.target.source_type_apt = 'EXTENDED'" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "9ad67489", + "metadata": {}, + "outputs": [], + "source": [ + "# Read in previously extracted regions (this time the annuli)\n", + "\n", + "reg_outer = Regions.read('87a_extract_outer.reg', format='ds9')[0]\n", + "reg_inner = Regions.read('87a_extract_inner.reg', format='ds9')[0]\n", + "print(reg_outer)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "448827ac", + "metadata": {}, + "outputs": [], + "source": [ + "# Create a weight map using the region mask for the outer annulus\n", + "\n", + "tmp_wgts = spec_model_cube.weightmap[:]\n", + "mask = reg_outer.to_mask('exact')\n", + "x1 = int(reg_outer.center.x-reg_outer.width/2.)\n", + "x2 = x1+mask.shape[1]\n", + "y1 = int(reg_outer.center.y-reg_outer.height/2.)\n", + "y2 = y1+mask.shape[0]\n", + "\n", + "### Note above, the region shape is slightly different than the mask shape that gets generated. That's why I had to put in this hack about setting the size.\n", + "\n", + "# Start by setting all pixels to 0.\n", + "mask2d = tmp_wgts[13,:,:]*0.\n", + "\n", + "# Because we want an inverse array, we can't just use the mask, we have to subtract the mask (which is 1's) from the original mask2d above (make sense?) \n", + "mask2d[y1:y2,x1:x2] = mask2d[y1:y2,x1:x2]+mask\n", + "\n", + "# Take into account weird rounding errors\n", + "mask2d[mask2d>0.1] = 1." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "2bd3c5de", + "metadata": {}, + "outputs": [], + "source": [ + "# Create a weight map using the region mask for the inner annulus (and subtract it)\n", + "\n", + "tmp_wgts = spec_model_cube.weightmap[:]\n", + "mask = reg_inner.to_mask('exact')\n", + "x1 = int(reg_inner.center.x-reg_inner.width/2.)\n", + "x2 = x1+mask.shape[1]\n", + "y1 = int(reg_inner.center.y-reg_inner.height/2.)\n", + "y2 = y1+mask.shape[0]\n", + "\n", + "### Note above, the region shape is slightly different than the mask shape that gets generated. That's why I had to put in this hack about setting the size.\n", + "\n", + "# mask2d already defined in cell above\n", + "\n", + "# Because we want an inverse array, we can't just use the mask, we have to subtract the mask (which is 1's) from the original mask2d above (make sense?) \n", + "mask2d[y1:y2,x1:x2] = mask2d[y1:y2,x1:x2]-mask\n", + "\n", + "# Take into account weird rounding errors\n", + "mask2d[mask2d<0.1] = 0.\n", + "npix = len(mask2d[mask2d==1])\n", + "npix" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "c0fa7cb5", + "metadata": {}, + "outputs": [], + "source": [ + "# Visualize the 2D Mask\n", + "\n", + "from astropy.nddata import CCDData\n", + "from astropy.visualization import simple_norm\n", + "ccd = mask2d\n", + "norm = simple_norm(ccd, 'sqrt', min_cut=0, max_cut=0.5) \n", + "color = 'rgbmkrgbmk'\n", + "\n", + "xceni = [36, 44]\n", + "yceni = [66, 58]\n", + "\n", + "fig = plt.figure(figsize=(10, 10))\n", + "ax = fig.add_subplot()\n", + "counter = 0\n", + "plt.title(\"Masked Data\")\n", + "plt.imshow(ccd, norm=norm, origin=\"lower\") \n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "e2a26371", + "metadata": {}, + "outputs": [], + "source": [ + "# Create a 3D weightmap from the 2D map. Mask all NAN values, too.\n", + "\n", + "mask3d = np.broadcast_to(mask2d, spec_model_cube.weightmap.shape)\n", + "mask3d.flags.writeable = True\n", + "mask3d[np.isnan(spec_model_cube.data)] = 0\n", + "#mask_sci_cube = np.ma.masked_array(spec_model_cube.weightmap, mask=mask3d.astype(bool))\n", + "tmp_wgt_cube = np.swapaxes(mask3d,0,1)\n", + "tmp_wgt_cube = np.swapaxes(tmp_wgt_cube,1,2)\n", + "plotcube = Spectrum1D(tmp_wgt_cube*u.dimensionless_unscaled)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "badea9cf", + "metadata": {}, + "outputs": [], + "source": [ + "# Visualize the 3D Mask in Cubeviz\n", + "\n", + "cubeviz3 = Cubeviz()\n", + "cubeviz3.load_data(plotcube, data_label='SN1987A MASK')\n", + "cubeviz3.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "18cd4752", + "metadata": {}, + "outputs": [], + "source": [ + "# Define the weightmap in the background subtracted cube as the new 3D Mask \n", + "# This will tell the pipeline which spaxels to use for extraction\n", + "# Reassign the weightmap to the original cube and save it\n", + "spec_model_cube.weightmap = mask3d\n", + "spec_model_cube.save(spec3_dir+'87A_bg_sub_annulus_extract.fits',overwrite=True)" + ] + }, + { + "cell_type": "markdown", + "id": "914bc219", + "metadata": {}, + "source": [ + "# 5. Extract SN Equitorial Ring Spectrum using Extract1dStep" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "ab919dd7", + "metadata": {}, + "outputs": [], + "source": [ + "# We will consider this our master background extraction\n", + "\n", + "cubefile_p1 = spec3_dir+'87A_bg_sub_annulus_extract.fits'\n", + "\n", + "outputfile = spec3_dir+'87A_bg_sub_annulus_extract.fits'\n", + "runex(cubefile_p1,spec3_dir,outputfile=outputfile)" + ] + }, + { + "cell_type": "markdown", + "id": "7e8b61e9", + "metadata": {}, + "source": [ + "# 6. Plot Final Results" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "6fdf7de1", + "metadata": {}, + "outputs": [], + "source": [ + "# First Plot Background in the Background Subtracted Cube\n", + "\n", + "filelist = glob.glob(spec3_dir+'87A_bg_sub_skymask_extract1dstep.fits')\n", + "filelist" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "06e13c9d", + "metadata": {}, + "outputs": [], + "source": [ + "# MRS Header Keyword value for the pixel area in sterdians \n", + "\n", + "ref_cube = spec3_dir+'87A_bg_sub_annulus_extract.fits'\n", + "hdu = fits.open(ref_cube)\n", + "hdr = hdu[1].header\n", + "pixar_sr = hdr['PIXAR_SR'] #should be 3.97217570860291E-13\n", + "print(pixar_sr)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "56cc11ec", + "metadata": {}, + "outputs": [], + "source": [ + "# Make Plot\n", + "\n", + "from specutils.manipulation import box_smooth, gaussian_smooth, trapezoid_smooth\n", + "\n", + "fig, ax = plt.subplots(figsize=(6, 6))\n", + "counter = 0\n", + "\n", + "#color = cm.rainbow(np.linspace(0, 1, len(filelist)))\n", + "color = 'rgbmkrgbmk'\n", + "ylim_low = -1.e2\n", + "ylim_high = 5.e-2\n", + "for filename in filelist:\n", + " jpipe_x1d = Spectrum1D.read(filename)\n", + " jpipe_x1d = gaussian_smooth(jpipe_x1d,stddev=50)\n", + " ax.plot(jpipe_x1d.spectral_axis, jpipe_x1d.flux, color = color[counter], label = filename)\n", + " ax.set_title(\"Manual Background Extractions AFTER Master Background Subtraction\")\n", + " ax.set_xlabel(r\"wavelength [$\\mu$m]\")\n", + " ax.set_ylabel(\"Flux Density [MJy/Str]\")\n", + " #ax.set_yscale(\"log\")\n", + " ax.set_ylim(ylim_low, ylim_high)\n", + " #ax.set_xlim(xlim_low, xlim_high)\n", + " ax.legend(bbox_to_anchor=(1.1, 1.05))\n", + " counter = counter+1\n", + "\n", + "plt.savefig(spec3_dir+'87A_background.png')\n", + "\n", + "jpipe_x1d.flux[0:20]" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "c14eb1f2", + "metadata": {}, + "outputs": [], + "source": [ + "# Next Plot the SN 1987A Equitorial Ring from the Background Subtracted Cube\n", + "\n", + "filelist = glob.glob(spec3_dir+'87A_bg_sub_annulus_extract_extract1dstep.fits')\n", + "filelist" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "f0c700b6", + "metadata": {}, + "outputs": [], + "source": [ + "# MRS Header Keyword value for the pixel area in sterdians \n", + "\n", + "ref_cube = spec3_dir+'87A_bg_sub_annulus_extract.fits'\n", + "hdu = fits.open(ref_cube)\n", + "hdr = hdu[1].header\n", + "pixar_sr = hdr['PIXAR_SR'] #should be 3.97217570860291E-13\n", + "print(pixar_sr)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "13b57670", + "metadata": {}, + "outputs": [], + "source": [ + "# Make Plot\n", + "\n", + "# Remember, the pipeline treats an EXTENDED source by averaging all the spaxels. \n", + "# But we want the total number flux, so we need to multiply by the number of spaxels in the annulus\n", + "# We calculated this value above when we were creating our mask (npix)\n", + "\n", + "from specutils.manipulation import box_smooth, gaussian_smooth, trapezoid_smooth\n", + "\n", + "fig, ax = plt.subplots(figsize=(6, 6))\n", + "counter = 0\n", + "\n", + "#color = cm.rainbow(np.linspace(0, 1, len(filelist)))\n", + "color = 'rgbmkrgbmk'\n", + "ylim_low = 1.e-4\n", + "ylim_high = 1.5e-1\n", + "for filename in filelist:\n", + " jpipe_x1d = Spectrum1D.read(filename)\n", + " #jpipe_x1d = gaussian_smooth(jpipe_x1d,stddev=50)\n", + " flux_convert = jpipe_x1d.flux*1.e6*pixar_sr*npix # Convert from MJy/str to Jy\n", + " ax.plot(jpipe_x1d.spectral_axis, flux_convert, color = color[counter], label = filename)\n", + " ax.set_title(\"SN 1987A Ring AFTER Master Background Subtraction\")\n", + " ax.set_xlabel(r\"wavelength [$\\mu$m]\")\n", + " ax.set_ylabel(\"Flux Density [Jy]\")\n", + " #ax.set_yscale(\"log\")\n", + " ax.set_ylim(ylim_low, ylim_high)\n", + " #ax.set_xlim(xlim_low, xlim_high)\n", + " ax.legend(bbox_to_anchor=(1.1, 1.05))\n", + " counter = counter+1\n", + "\n", + "plt.savefig(spec3_dir+'87A_ring_bg_sub.png')\n", + "\n", + "jpipe_x1d.flux[0:20]" + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3 (ipykernel)", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.10.12" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/notebooks/MIRI_MRS_1987A/requirements.txt b/notebooks/MIRI_MRS_1987A/requirements.txt new file mode 100644 index 000000000..0e0b08691 --- /dev/null +++ b/notebooks/MIRI_MRS_1987A/requirements.txt @@ -0,0 +1,4 @@ +git+https://github.com/spacetelescope/jdaviz.git +git+https://github.com/spacetelescope/jwst.git +astropy >= 5.3.1 +astroquery