diff --git a/notebooks/MIRI/MIRI_MRS_background_subtraction/MIRI_MRS_bkg_subtraction_with_AstroBkgInterp.ipynb b/notebooks/MIRI/MIRI_MRS_background_subtraction/MIRI_MRS_bkg_subtraction_with_AstroBkgInterp.ipynb new file mode 100644 index 000000000..adf6da7dc --- /dev/null +++ b/notebooks/MIRI/MIRI_MRS_background_subtraction/MIRI_MRS_bkg_subtraction_with_AstroBkgInterp.ipynb @@ -0,0 +1,835 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "a4143a8e-03a4-4d02-a883-fc617db11be9", + "metadata": {}, + "source": [ + "# MIRI MRS Background Subtraction with AstroBkgInterp\n", + "\n", + "**Use case:** Background estimation and subtraction on Level2 data using dedicated package `AstroBkgInterp` (https://github.com/brynickson/AstroBkgInterp). `AstroBkgInterp` is a Python-based tool that performs flexible and accurate background estimation and subtraction for both 2D images and 3D data cubes. It is designed to address limitations in traditional background methods, particularly in environments where the background is spatially structured, the source is embedded in high surface-brightness regions, or no dedicated background frames are available.
\n", + "**Data**: JWST/MIRI Medium Resolution Spectrograph (MRS) observations of the Type IIn Supernova 2005ip.
\n", + "**Tools**: AstroBkgInterp, numpy, astropy, matplotlib
\n", + "**Intrument**: MIRI.
\n", + "**Documentation**: This notebook is part of STScI’s larger [post-pipeline Data Analysis Tools Ecosystem](https://jwst-docs.stsci.edu/jwst-post-pipeline-data-analysis) and can be downloaded directly from the [JDAT Notebook Github directory](https://github.com/spacetelescope/jdat_notebooks).
\n", + "**Author**: Bryony Nickson (STScI)
\n", + "**Date Published**: 06/27/25
\n", + "**Last Updated**: 09/11/25
\n", + "\n", + "## Contents\n", + "1. [Introduction](#intro)
\n", + " 1.1. [Package Imports](#imports)
\n", + " 1.2. [Setup](#setup)
\n", + "2. [Source-Masking](#source-masking)
\n", + " 2.1. [Locate the source using `Photutils`](#loc)
\n", + " 2.2. [Set the aperture and annulus sizes](#set)
\n", + " 2.3. [Verify the source-masking](#verify)
\n", + "3. [Background Modelling](#bkg)
\n", + " 3.1. [Estimate and subtract the background](#run)
\n", + "4. [Final Spectrum Extraction](#results)
\n", + " 4.1. [Run `Extract1d` on the background subtracted data](#extract1d)
\n", + " 4.2. [Inspect the final spectrum](#plot)
\n", + " 4.3. [Compare with the default pipeline extraction](#compare)
\n" + ] + }, + { + "cell_type": "markdown", + "id": "20b06d3e-2b8e-4bd1-b78a-9c9b84c67533", + "metadata": {}, + "source": [ + "\n", + "## Introduction \n", + "Core-collapse supernovae (CCSNe) are critical probes of massive stellar evolution and chemical enrichment in galaxies. However, they are often embedded in crowded or high surface-brightness regions of their host galaxies, where background emission is spatially structured and varies on small angular scales. In such environments, accurate background estimation is essential for isolating the supernova signal and deriving reliable photometry and spectra. \n", + "\n", + "This notebook analyzes observations of the Type IIn Supernova 2005ip captured 17 years post-explosion by [Shahbandeh et al. 2024](https://doi.org/10.48550/arXiv.2410.09142). For 2005ip, the global thermal background is roughly 10 to 40 times larger than the Supernova (SN) itself, making local background estimation both critical and challenging. Even small variability in the background across the field of view (FOV) can have a potentially significant impact on the source spectrum. This is particularly true at the longer wavelengths, where thermal background tends to dominate, and the intrinsic SN signal is weakest.\n", + "\n", + "To address this, we will use `AstroBkgInterp` to perform accurate, spatially resolved background subtraction tailored to the complex environment of SN 2005ip. First, we will apply its source-masking procedure to exclude the SN flux and replace it with an interpolated estimate of the underlying background. Then, we will use `AstroBkgInterp` to fit a global background model to the masked data and subtract it from the original data. This approach will improve the fidelity of the extracted spectrum, particularly in the background-dominated mid-infrared regime, enabling the detection of faint dust emission features and yielding more robust constraints on the cold dust mass. " + ] + }, + { + "cell_type": "markdown", + "id": "20782ead-244d-411d-8c58-555d876e2a62", + "metadata": { + "tags": [] + }, + "source": [ + "\n", + "### Import Packages\n", + "\n", + "- `AstroBkgInterp` is our background estimation tool\n", + "- `numpy` for array processing and math\n", + "- `atropy.io` for accessing the data\n", + "- `astropy.time` for timing\n", + "- `astropy.stats` for calculating statistics on the data\n", + "- `matplotlib` for plotting images and spectra\n", + "- `photutils.detection` for finding sources in the data\n", + "- `jwst` for running the pipeline" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "af6e8895-d8fa-4c35-a6e1-e8010cac6d80", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "import copy\n", + "import os\n", + "\n", + "# Import the background subtraction tool \n", + "from AstroBkgInterp.AstroBkgInterp import AstroBkgInterp\n", + "\n", + "# Import astropy packages\n", + "from astropy.io import fits\n", + "from astropy.table import Table\n", + "from astropy.stats import sigma_clipped_stats\n", + "import astropy.units as u\n", + "\n", + "# Import packages for displaying images in notebook\n", + "from matplotlib import pyplot as plt\n", + "from matplotlib.patches import Circle\n", + "from mpl_toolkits.axes_grid1 import make_axes_locatable\n", + "\n", + "# For data handling\n", + "import numpy as np\n", + "\n", + "# To find stars in the MRS spectralcubes \n", + "from photutils.detection import DAOStarFinder\n", + "\n", + "# For running the pipeline\n", + "from jwst.extract_1d import Extract1dStep\n", + "from jwst import datamodels\n", + "os.environ['CRDS_SERVER_URL'] = 'https://jwst-crds.stsci.edu'" + ] + }, + { + "cell_type": "markdown", + "id": "36fc1a34-4512-4e58-9bdc-349ec0712e6d", + "metadata": {}, + "source": [ + "\n", + "### Setup\n", + "\n", + "#### Set path to Data" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "66057c1b-07bd-45f2-a2b5-0635862a76a4", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "data_url = \"https://data.science.stsci.edu/redirect/JWST/jwst-data_analysis_tools/MIRI_MRS_bkg_notebook/Level3_ch1-2-3-4-shortmediumlong_s3d.fits\"" + ] + }, + { + "cell_type": "markdown", + "id": "eb4682e7-13b8-45d1-9ecd-cf48d92da781", + "metadata": {}, + "source": [ + "#### Open and display the data " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "0df31923-4b87-4206-89c0-dd5b87f6e1e0", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "hdu = fits.open(data_url)\n", + "data = hdu[1].data\n", + "\n", + "# set all NaN values to 0\n", + "data[np.isnan(data)] = 0\n", + "\n", + "plt.imshow(data[9000], origin='lower')\n", + "plt.colorbar()" + ] + }, + { + "cell_type": "markdown", + "id": "1e4e1fc4-c4b6-49f2-a162-3541d76b05e8", + "metadata": {}, + "source": [ + "\n", + "## Source-Masking Procedure\n", + "\n", + "A critical feature of `AstroBkgInterp` is its ability to exclude flux from astrophysical sources during the background estimation process. It implements a flexible and configurable source-masking routine that enables users to define regions of interest (typically containing one or more sources) to be excluded from the final background model. \n", + "\n", + "The user defines the position of the target source, along with an aperture geometry and size; both circular and elliptical apertures are supported. These apertures define the region to be masked and replaced with an interpolated estimate of the underlying background. In polar coordinates centered on the source, the algorithm computes the median flux at each angle from the surrounding annulus and interpolates across the masked region. This directionally informed approach preserves local background structure while excluding contamination from the source itself. To increase robustness and minimize interpolation artifacts, the procedure is repeated in a grid of eight dithered positions surrounding the central source, and the resulting estimates are median-combined to produce the final background.\n", + "\n", + "For 3D data cubes, `AstroBkgInterp` supports wavelength-dependent masking, allowing for the mask to vary as a function of wavelength. This capability accommodates the variation in PSF size with wavelength, and spectral variations in source morphology, ensuring accurate background estimation across the full spectral range." + ] + }, + { + "cell_type": "markdown", + "id": "8466a853-ce49-4a6d-a144-61a65f7e970f", + "metadata": {}, + "source": [ + "\n", + "### 1) Locate the source using `Photutils` \n", + "\n", + "In order to utilize `AstroBkgInterp`'s source-masking algorithm, we first need to identify the position of the source. Rather than relying solely on the coordinates provided in the image header, we begin by identifying the brightest source in the data cube, which we assume to be the supernova (SN). To accomplish this, we use `DAOStarFinder`, an implementation of the DAOFIND algorithm [Stetson 1987](https://ui.adsabs.harvard.edu/abs/1987PASP...99..191S/abstract) designed for point-source detection in astronomical images. DAOFIND locates local maxima in the image that exceed a specified `threshold` (applied to a convolved image) and have a size and shape similar to the defined 2D Gaussian kernel. \n", + "\n", + "**Note**: both the detection threshold and maximum separation allowed between neighboring sources should be adjusted based on the specific science case (e.g., the number of stars in the field, the degree of crowding, and the expected brightness distribution)." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "00ca64d5-5237-4905-85be-ebd276988b5b", + "metadata": {}, + "outputs": [], + "source": [ + "cube = np.zeros((data.shape[1], data.shape[2]))\n", + "for a in range(data.shape[1]):\n", + " for b in range(data.shape[2]):\n", + " cube[a, b] = np.median(data[:, a, b])\n", + "\n", + "mean, median, std = sigma_clipped_stats(cube, sigma=3.0)\n", + "\n", + "# Get a list of sources using a dedicated source detection algorithm\n", + "# Find sources at least 3* background (typically)\n", + "\n", + "daofind = DAOStarFinder(fwhm=2.0, threshold=3.*std)\n", + "sources = daofind(cube-median) \n", + "print(\"\\n Number of sources in field:\", len(sources))\n", + "\n", + "# Positions in pixels\n", + "positions = Table([sources['xcentroid'], sources['ycentroid']])\n", + "\n", + "# Convert to RA & Dec (ICRS)\n", + "peakpixval = np.zeros(len(sources['xcentroid']))\n", + "\n", + "for count_s, _ in enumerate(sources):\n", + " peakpixval[count_s] = cube[int(np.round(sources['xcentroid'][count_s])), int(np.round(sources['ycentroid'][count_s]))]\n", + "\n", + "# Set the peak pixel positions as the source position. \n", + "src_x, src_y = sources['xcentroid'][np.argmax(peakpixval)], sources['ycentroid'][np.argmax(peakpixval)]\n", + "print(f'peak pixel x = {src_x}')\n", + "print(f'peak pixel y = {src_y}')" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "70b105e1-69a5-4685-a7cc-790a048e624b", + "metadata": {}, + "outputs": [], + "source": [ + "# Plot all of the sources\n", + "plt.imshow(data[9000], origin='lower')\n", + "plt.colorbar()\n", + "plt.scatter(sources['xcentroid'], sources['ycentroid'], c=\"red\", marker=\"+\", s=50)\n", + "plt.scatter(src_x, src_y, c=\"black\", marker=\"+\", s=50)" + ] + }, + { + "cell_type": "markdown", + "id": "5cb5f34e-80c3-463b-92fd-59bbf182d6b1", + "metadata": {}, + "source": [ + "\n", + "### 2) Set the aperture and annulus sizes\n", + "\n", + "The next step is to define the aperture and annulus radii for the source-masking procedure. The aperture defines the region to be masked (i.e., the area occupied by the source), while the annulus provides the surrounding region from which the background underneath the source is estimated.\n", + "\n", + "When selecting an aperture radius, the goal is to encompass as much of the source flux as possible, ensuring that the annulus lies beyond the majority of the source light. At the same time, the annulus should remain close enough to the source that the sampled background is representative of the local environment. The annulus must also be sufficiently wide to include enough pixels for a reliable statistical estimate of the underlying background.\n", + "\n", + "#### Manual aperture definition" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "0023f62a-b131-40c4-9517-72b88dadb67d", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "aper_rad = 8 \n", + "ann_width = 4\n", + "\n", + "plt.figure(figsize=(5, 5))\n", + "plt.imshow(data[8000], vmin=-50, vmax=175, origin='lower')\n", + "plt.colorbar()\n", + "\n", + "circ = Circle((src_x, src_y), radius=aper_rad, color='r', fill=False)\n", + "annin = Circle((src_x, src_y), radius=aper_rad+ann_width, color='r', fill=False)\n", + "plt.gca().add_patch(circ)\n", + "plt.gca().add_patch(annin)" + ] + }, + { + "cell_type": "markdown", + "id": "09d073ae-768e-4ce8-9e07-1e76be072b09", + "metadata": {}, + "source": [ + "#### Wavelength dependent masking\n", + "\n", + "For 3D data cubes, `AstroBkgInterp` also supports wavelength-dependent masking, allowing the aperture radius to vary as a function of wavelength. This capability accommodates the variation in PSF size with wavelength, and spectral variations in source morphology, ensuring accurate background estimation across the full spectral range. \n", + "\n", + "To define an aperture mask size that appropriately encompasses the flux, users can base the aperture radius on the PSF FWHM. To do this, they will need to configure two parameters:\n", + "\n", + "1. `fwhm`: the PSF FWHM array, i.e., an array matching the spectral dimension of the data cube, where each element corresponds to the PSF FWHM at that wavelength. This allows `AstroBkgInterp` to adjust the aperture size dynamically across the cube.\n", + "2. `fwhm_scale`: a constant scaling factor (default is 1.25) that multiplies the PSF FWHM to set the aperture radius. \n", + "\n", + "Together, these parameters will set the aperture radius to `fwhm_scale x fwhm` at each wavelength slice.\n", + "\n", + "Below we will define the FWHM array using WCS parameters and `np.linspace`, and the linear fit to the MIRI PSF found in Figure 2 https://jwst-docs.stsci.edu/jwst-mid-infrared-instrument/miri-operations/miri-dithering/miri-mrs-psf-and-dithering#gsc.tab=0." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "d943d49e-fdda-462b-8fdd-fe3871c6aff8", + "metadata": {}, + "outputs": [], + "source": [ + "# Extract wavelength array from WCS metadata table\n", + "WAVEfull = hdu[5].data['WAVELENGTH'][0]\n", + "\n", + "dat = datamodels.open(data_url)\n", + "\n", + "# Get spatial pixel scale along x-axis from the WCS metadata (degrees/pixel)\n", + "cdelt1 = dat.meta.wcsinfo.cdelt1*u.deg\n", + "\n", + "# Convert pixel scale from deg/pix to arcsec/pix\n", + "cdelt1 = cdelt1.to(u.arcsec).value\n", + "\n", + "# Compute angular PSF FWHM (θ) in arcsec as a function of wavelength.\n", + "# This uses an empirical model: θ(λ) = 0.033 * λ + 0.106, with λ in microns.\n", + "θ = 0.033 * (WAVEfull) + 0.106 \n", + "\n", + "# Convert FWHM from arcsec to pixels by dividing by the pixel scale\n", + "fwhm = (θ / cdelt1) \n", + "\n", + "# FWHM array of PSF size in pixels as a function of wavelength\n", + "fwhm\n", + "\n", + "# Set the FWHM scale\n", + "fwhm_scale = 1.5 # This will produce an aperture radius 1.25 times the local PSF FWHM." + ] + }, + { + "cell_type": "markdown", + "id": "a9dc9b3c-944e-46d0-8bb2-ed79189bfb6e", + "metadata": {}, + "source": [ + "Below, we compare the two aperture masking approaches across different slices of the cube: the top row shows the fixed-radius method that applies the same aperture size at all wavelengths, and the second row shows the wavelength-dependent method that dynamically adjusts the aperture size at each spectral slice based on the `fwhm` and `fwhm_scale` we set above. " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "b7b4ce38-7142-4ed5-af62-5b270678cb49", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "fig, ax = plt.subplots(2, 3, figsize=(8, 6), sharex=True)\n", + "\n", + "circ = Circle((src_x, src_y), radius=aper_rad, color='r', fill=False)\n", + "\n", + "fig.suptitle('Set vs Wavelength Dependent Aperture Masking')\n", + "\n", + "ax[0, 0].set_title('Slice 1000')\n", + "ax[0, 1].set_title('Slice 5000')\n", + "ax[0, 2].set_title('Slice 8000')\n", + "ax[0, 0].set_ylabel('Set aperture radius', fontsize=12)\n", + "ax[1, 0].set_ylabel('Wavelength-dependent radius', fontsize=12)\n", + "\n", + "ax[0, 0].imshow(data[1000], vmin=-50, vmax=175, origin='lower')\n", + "circ = Circle((src_x, src_y), radius=aper_rad, color='r', fill=False)\n", + "annin = Circle((src_x, src_y), radius=aper_rad+ann_width, color='r', fill=False)\n", + "ax[0, 0].add_patch(circ)\n", + "ax[0, 0].add_patch(annin)\n", + "\n", + "ax[0, 1].imshow(data[5000], vmin=-50, vmax=175, origin='lower')\n", + "circ = Circle((src_x, src_y), radius=aper_rad, color='r', fill=False)\n", + "annin = Circle((src_x, src_y), radius=aper_rad+ann_width, color='r', fill=False)\n", + "ax[0, 1].add_patch(circ)\n", + "ax[0, 1].add_patch(annin)\n", + "\n", + "ax[0, 2].imshow(data[8000], vmin=-50, vmax=175, origin='lower')\n", + "circ = Circle((src_x, src_y), radius=aper_rad, color='r', fill=False)\n", + "annin = Circle((src_x, src_y), radius=aper_rad+ann_width, color='r', fill=False)\n", + "ax[0, 2].add_patch(circ)\n", + "ax[0, 2].add_patch(annin)\n", + "\n", + "ax[1, 0].imshow(data[1000], vmin=-50, vmax=175, origin='lower')\n", + "aper_rad = fwhm_scale*fwhm[1000]\n", + "circ = Circle((src_x, src_y), radius=aper_rad, color='r', fill=False)\n", + "annin = Circle((src_x, src_y), radius=aper_rad+ann_width, color='r', fill=False)\n", + "ax[1, 0].add_patch(circ)\n", + "ax[1, 0].add_patch(annin)\n", + "\n", + "ax[1, 1].imshow(data[5000], vmin=-50, vmax=175, origin='lower')\n", + "aper_rad = fwhm_scale*fwhm[5000]\n", + "circ = Circle((src_x, src_y), radius=aper_rad, color='r', fill=False)\n", + "annin = Circle((src_x, src_y), radius=aper_rad+ann_width, color='r', fill=False)\n", + "ax[1, 1].add_patch(circ)\n", + "ax[1, 1].add_patch(annin)\n", + "\n", + "im = ax[1, 2].imshow(data[8000], vmin=-50, vmax=175, origin='lower')\n", + "aper_rad = fwhm_scale*fwhm[8000]\n", + "circ = Circle((src_x, src_y), radius=aper_rad, color='r', fill=False)\n", + "annin = Circle((src_x, src_y), radius=aper_rad+ann_width, color='r', fill=False)\n", + "ax[1, 2].add_patch(circ)\n", + "ax[1, 2].add_patch(annin)\n", + "\n", + "plt.tight_layout()" + ] + }, + { + "cell_type": "markdown", + "id": "0a190815-62ca-400c-8498-d92299e60227", + "metadata": {}, + "source": [ + "The above plot illustrates how applying a fixed aperture across the entire MIRI MRS data cube can result in both over-masking at some wavelengths and under-masking at others, ultimately compromising the quality of the background subtraction.\n", + "\n", + "At shorter wavelengths (where the PSF is narrower), the fixed aperture may be too large, including excess background and potentially nearby sources. Conversely, at longer wavelengths (where the PSF is broader), the same aperture may be too small, missing significant portions of the source flux. This leaves residual source flux in the background model, contaminating the background estimate and under-subtracting the source. \n", + "\n", + "Using wavelength-dependent masking helps address these issues by ensuring the aperture appropriately scales with the PSF across the spectral range, leading to cleaner and more accurate background estimates.\n" + ] + }, + { + "cell_type": "markdown", + "id": "23b56cbb-a903-4788-a9f4-14fadf2205e6", + "metadata": {}, + "source": [ + "\n", + "### 3) Verify the source-masking \n", + "\n", + "To verify the source masking and manually refine the aperture and annulus radii if needed, `AstroBkgInterp` can be run with `bkg_mode = None`. In this mode, the source-masking routine executes as usual, but the tool does not construct a spatial background model across the entire frame. The result is effectively a \"source replacement\" only, where the background is interpolated locally beneath each masked region without estimating or subtracting a broader background across the image." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "402154a5-095b-4f26-a6ac-b42b3116abb9", + "metadata": {}, + "outputs": [], + "source": [ + "bi = AstroBkgInterp()\n", + "\n", + "# Source position\n", + "bi.src_y = src_y\n", + "bi.src_x = src_x\n", + "\n", + "# Source masking params\n", + "bi.aper_rad = aper_rad\n", + "bi.ann_width = ann_width\n", + "\n", + "bi.fwhm = fwhm\n", + "bi.fwhm_scale = 1.5\n", + "\n", + "# Background params\n", + "bi.bkg_mode = 'None' \n", + "\n", + "# Multiprocessing params\n", + "bi.pool_size = 12 \n", + "\n", + "diff, bkg, mask = bi.run(data)" + ] + }, + { + "cell_type": "markdown", + "id": "a4408eef-0ad8-465d-9510-17f5582d6d66", + "metadata": {}, + "source": [ + "Inspect the source masked data." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "a5630814-3fb4-4a34-9b17-52cc385bfb2e", + "metadata": { + "scrolled": true + }, + "outputs": [], + "source": [ + "fig, ax = plt.subplots(1, 2, figsize=(6, 5), sharex=True)\n", + "\n", + "ax[0].set_title('Raw Data')\n", + "ax[1].set_title('Source Masked Data')\n", + "\n", + "a = ax[0].imshow(data[8000], vmin=-50, vmax=175, origin='lower', cmap='viridis')\n", + "b = ax[1].imshow(bkg[8000], vmin=-50, vmax=175, origin='lower', cmap='viridis')\n", + "\n", + "divider = make_axes_locatable(ax[1])\n", + "cax = divider.append_axes('right', size='5%', pad=0.05)\n", + "fig.colorbar(b, cax=cax, orientation='vertical')\n", + "\n", + "plt.tight_layout()\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "id": "18ff94d4-b4b5-4bd9-ba07-f90734b7b0e4", + "metadata": {}, + "source": [ + "Once we are happy with the source masking, we can move on to the background modelling stage." + ] + }, + { + "cell_type": "markdown", + "id": "b5e80b5a-cf82-4ecb-a780-5b021511b14f", + "metadata": {}, + "source": [ + "\n", + "## Background Modelling\n", + "\n", + "In addition to localized source masking, `AstroBkgInterp` provides robust background modelling capabilities designed to handle spatially complex astrophysical sources. After masking flux in the user-defined source regions, the tool constructs a background model using one of several configurable fitting methods. These include a fast, one-dimensional row/column median estimator (`bkg_mode = simple`) and a more sophisticated two-dimensional polynomial surface fitting approach (`bkg_mode = polynomial`). The background model is generated from the source-masked image and is applied on a per-frame basis for 2D images, and on a per-slice basis for 3D data cubes. " + ] + }, + { + "cell_type": "markdown", + "id": "de7a1699-6dc7-4a4f-ae24-0190f4bc4af5", + "metadata": {}, + "source": [ + "\n", + "### 4) Estimate and subtract the background\n", + "\n", + "For the purpose of this notebook, we will set the `\"polynomial\"` background mode. This mode constructs a background model using a 2D polynomial surface fitting approach implemented in a sliding-window framework. The image is divided into a series of overlapping subregions, and a low-order polynomial is independently fit to each region. These fits are then median combined at each pixel location, effectively mitigating the influence of outliers, cosmic rays, and localized noise features. The result is a smooth and robust background model that preserves large-scale background structure while expressing small-scale fluctuations. \n", + "\n", + "We set the following parameters to control the behavior and resolution of the polynomial fitting process:\n", + "\n", + "- `k`: Specifies the degree of the 2D polynomial used for each local fit. Higher values allow for greater flexibility in modeling complex background structure but may increase the risk of overfitting. We adopt the default of `k = 3`, which provides a good balance between smoothness and adaptability.\n", + "\n", + "- `bin_size`: Defines the size of the subregions (windows) over which the local polynomial fitting is performed. Smaller windows capture finer spatial variations but increase computational cost, while larger windows produce smoother fits at the expense of resolving small-scale structure. For this analysis, we use `bin size = 5`.\n", + "\n", + "- `cube_resolution`: Controls the spatial sampling and density of polynomial fits across the image. It determines both the number of local fits and the step size between successive fitting regions. We use `cube_resolution = medium`, which offers a balanced trade-off between sampling and computation time, and is recommended for most use cases." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "c24aa9d5-5d8c-495d-b882-127d6c415478", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "bi = AstroBkgInterp()\n", + "\n", + "# Source position\n", + "bi.src_y = src_y\n", + "bi.src_x = src_x\n", + "\n", + "# Source masking params\n", + "bi.aper_rad = aper_rad\n", + "bi.ann_width = ann_width\n", + "\n", + "bi.fwhm = fwhm\n", + "bi.fwhm_scale = 1.5\n", + "\n", + "# Background params\n", + "bi.bkg_mode = 'polynomial' \n", + "bi.k = 3 \n", + "bi.bin_size = 9 \n", + "bi.cube_resolution = 'medium' \n", + "\n", + "# Multiprocessing params\n", + "bi.pool_size = 12 \n", + "\n", + "diff, bkg, mask = bi.run(data)" + ] + }, + { + "cell_type": "markdown", + "id": "fb7cdba9-ee4a-41b9-9443-5851ea2f83c0", + "metadata": {}, + "source": [ + "Save out new background subtracted data." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "8fbdf100-e1c3-4d87-90b9-c3a7a39f5317", + "metadata": {}, + "outputs": [], + "source": [ + "newdata = np.array([s for s in diff])\n", + "newhdu = copy.deepcopy(hdu)\n", + "newhdu[1].data = newdata\n", + "newhdu.writeto('newdata_med_res.fits', overwrite=True)" + ] + }, + { + "cell_type": "markdown", + "id": "f8398d86-7d64-4eeb-8a58-f2b39653c053", + "metadata": {}, + "source": [ + "Now lets plot out all the products returned by `AstroBkgInterp`." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "d0e4d36b-a5ac-452a-8954-7e956e7010b1", + "metadata": {}, + "outputs": [], + "source": [ + "fig, ax = plt.subplots(2, 2, figsize=(10, 12))\n", + "\n", + "ax[0, 0].set_title('Raw Data')\n", + "ax[0, 1].set_title('Source-Masked data')\n", + "ax[1, 0].set_title('Background Fit')\n", + "ax[1, 1].set_title('Residual')\n", + "\n", + "a = ax[0, 0].imshow(data[8000], vmin=-10, vmax=180, origin='lower')\n", + "b = ax[0, 1].imshow(mask[8000][0], vmin=-10, vmax=180, origin='lower')\n", + "c = ax[1, 0].imshow(bkg[8000], vmin=-10, vmax=180, origin='lower')\n", + "d = ax[1, 1].imshow(diff[8000], vmin=-10, vmax=180, origin='lower')\n", + "\n", + "divider = make_axes_locatable(ax[1, 1])\n", + "cax = divider.append_axes('right', size='5%', pad=0.05)\n", + "fig.colorbar(d, cax=cax, orientation='vertical')\n", + "\n", + "divider = make_axes_locatable(ax[1, 1])\n", + "cax = divider.append_axes('right', size='5%', pad=0.05)\n", + "fig.colorbar(d, cax=cax, orientation='vertical')\n", + "\n", + "divider = make_axes_locatable(ax[1, 1])\n", + "cax = divider.append_axes('right', size='5%', pad=0.05)\n", + "fig.colorbar(d, cax=cax, orientation='vertical')\n", + "\n", + "divider = make_axes_locatable(ax[1, 1])\n", + "cax = divider.append_axes('right', size='5%', pad=0.05)\n", + "fig.colorbar(d, cax=cax, orientation='vertical')\n", + "\n", + "plt.tight_layout()\n", + "\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "id": "343f6306-7508-4aae-bda4-9f65b1a24a80", + "metadata": {}, + "source": [ + "\n", + "## Final Spectrum Extraction\n", + "\n", + "Finally, we extract the spectrum of the target source at the supernova position using the `Extract1dStep` from the JWST calibration pipeline. Instead of relying solely on the header coordinates, we use the refined source position identified earlier with `DAOSourceFinder`. The resulting 1D spectrum represents the final science product, with local background emission minimized through subtraction of our custom background model generated by `AstroBkgInterp`. We omit the default background subtraction step provided by the pipeline, as our tailored approach is considered sufficient to account for the background at the source location accurately.\n", + "\n", + "\n", + "### 5) Run `Extract1d` on the background-subtracted data" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "3643667e-1195-48ef-88b6-fb25e79a1fc4", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "step = Extract1dStep()\n", + "\n", + "cube = datamodels.open('newdata_med_res.fits')\n", + "\n", + "result = step.call(cube, \n", + " subtract_background=False, \n", + " center_xy=[src_x, src_y],\n", + " ifu_rfcorr=True)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "5c48def2-cd58-43e1-bacf-a88e96a3ba1b", + "metadata": {}, + "outputs": [], + "source": [ + "result.to_fits('newdata_med_res_spec2.fits', overwrite=True)\n", + "res_pipe = fits.open('newdata_med_res_spec2.fits')" + ] + }, + { + "cell_type": "markdown", + "id": "74281b6a-5220-4794-97eb-16f7b62a6f85", + "metadata": {}, + "source": [ + "\n", + "### 6) Plot the final spectrum" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "4ece7891-5264-4e80-8ae2-1f15c0fa82eb", + "metadata": {}, + "outputs": [], + "source": [ + "plt.figure(figsize=(15, 9))\n", + "plt.tick_params(size=7, width=2, direction='inout', labelsize=12)\n", + "\n", + "spec = res_pipe[1].data\n", + "WAVE = spec['WAVELENGTH']\n", + "FLUX = spec['FLUX']\n", + "FLUX_mjy = (FLUX*u.Jy).to(u.mJy)\n", + "\n", + "plt.plot(WAVE, FLUX_mjy, lw=0.5, label='2D Interp Bkg')\n", + "\n", + "plt.title('AstroBkgInterp Background Subtracted Spectrum', fontsize=20)\n", + "\n", + "plt.ylim(-3, 10)\n", + "plt.xlim(4.8, 28)\n", + "\n", + "plt.xlabel(r'$\\mu m$', fontsize=15)\n", + "plt.ylabel('Flux (mJy)', fontsize=15)\n", + "\n", + "plt.xscale('linear')\n", + "plt.yscale('linear')\n", + "\n", + "plt.tight_layout()\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "id": "f7eee995-7596-4206-ba77-349e5c50af8c", + "metadata": {}, + "source": [ + "\n", + "### 7) Compare with the default JWST pipeline spectrum \n", + "\n", + "Here, we will run the original data through the JWST pipeline (with the background subtraction method used by the pipeline turned on) and compare it with the AstroBkgInterp-subtracted spectra. " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "419c47e0-3cbe-4dfe-a010-e446a57a2530", + "metadata": {}, + "outputs": [], + "source": [ + "# Original pipeline data\n", + "origdata = datamodels.open('Level3_ch1-2-3-4-shortmediumlong_s3d.fits')\n", + "\n", + "result = step.call(origdata, \n", + " subtract_background=True, \n", + " center_xy=[src_x, src_y],\n", + " ifu_rfcorr=True)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "7d7390db-dc75-477c-9a3b-94e10345ced2", + "metadata": {}, + "outputs": [], + "source": [ + "result.to_fits('origdata_spec2.fits', overwrite=True)\n", + "orig_res = fits.open('origdata_spec2.fits')" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "a38b5be9-4ee7-4859-873b-be548be9b89c", + "metadata": {}, + "outputs": [], + "source": [ + "fig, ax = plt.subplots(2, 1, figsize=(15, 9), sharex=True)\n", + "ax[0].tick_params(size=7, width=2, direction='inout', labelsize=12)\n", + "ax[1].tick_params(size=7, width=2, direction='inout', labelsize=12)\n", + "\n", + "origspec = orig_res[1].data\n", + "origWAVE = origspec['WAVELENGTH']\n", + "origFLUX = origspec['FLUX']\n", + "origFLUX_mjy = (origFLUX*u.Jy).to(u.mJy)\n", + "\n", + "ax[0].plot(origWAVE, origFLUX, lw=0.5, c='c', label='Pipeline')\n", + "ax[0].plot(WAVE, FLUX, lw=0.5, c='m', label='ABI')\n", + "ax[1].plot(origWAVE, origFLUX-FLUX, lw=0.5, c='k', label='Difference')\n", + "\n", + "ax[0].legend()\n", + "ax[1].legend()\n", + "\n", + "ax[0].set_title('Background Subtraction comparison', fontsize=20)\n", + "\n", + "ax[0].set_ylim(1e-4, 5e-2)\n", + "ax[1].set_ylim(-0.012, 0.012)\n", + "\n", + "ax[0].set_xlim(4.8, 28)\n", + "ax[1].set_xlim(4.8, 28)\n", + "\n", + "plt.xlabel(r'$\\mu m$', fontsize=15)\n", + "ax[0].set_ylabel('Flux (Jy)', fontsize=15)\n", + "ax[1].set_ylabel('Flux (Jy)', fontsize=15)\n", + "\n", + "ax[1].axhline(0, ls='--', c='r')\n", + "\n", + "ax[0].set_xscale('log')\n", + "ax[0].set_yscale('log')\n", + "ax[1].set_yscale('symlog')\n", + "\n", + "ax[0].set_xticks([5, 7.5, 10, 15, 20, 25])\n", + "ax[0].set_xticklabels([5, 7.5, 10, 15, 20, 25])\n", + "\n", + "ax[1].set_xticks([5, 7.5, 10, 15, 20, 25])\n", + "ax[1].set_xticklabels([5, 7.5, 10, 15, 20, 25])\n", + "\n", + "ax[1].set_yticks([-1e-2, -5e-3, -1e-3, 1e-3, 5e-3, 1e-2])\n", + "ax[1].set_yticklabels([-1e-2, -5e-3, -1e-3, 1e-3, 5e-3, 1e-2])\n", + "\n", + "plt.tight_layout()\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "id": "cc4de85c-b7c0-4030-b317-c42d494257b5", + "metadata": {}, + "source": [ + "A comparison between the default JWST pipeline output (teal) and the output after applying background subtraction with `AstroBkgInterp` (pink) highlights the improvement in background removal at the source position. This enhancement is particularly pronounced at longer wavelengths, where thermal background emission from the observatory dominates. Accurate subtraction in this regime is essential for isolating intrinsic dust emission and for placing meaningful constraints on cold dust mass and the underlying dust models." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "f034b6e4-ad5f-42b4-a496-61c71a238281", + "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.11.11" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/notebooks/MIRI/MIRI_MRS_background_subtraction/requirements.txt b/notebooks/MIRI/MIRI_MRS_background_subtraction/requirements.txt new file mode 100644 index 000000000..4a99effeb --- /dev/null +++ b/notebooks/MIRI/MIRI_MRS_background_subtraction/requirements.txt @@ -0,0 +1,6 @@ +jwst >= 1.18.1 +astropy >= 7.1.0 +photutils >= 2.2.0 +numpy >= 2.2.4 +matplotlib >= 3.10.1 +git+https://github.com/brynickson/AstroBkgInterp.git