diff --git a/notebooks/psf_photometry/MIRI/miri_phot_comparison.ipynb b/notebooks/psf_photometry/MIRI/miri_phot_comparison.ipynb new file mode 100644 index 000000000..eb64326d1 --- /dev/null +++ b/notebooks/psf_photometry/MIRI/miri_phot_comparison.ipynb @@ -0,0 +1,141 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "5e6982fe-42fe-42d7-984b-226ae08a01d0", + "metadata": {}, + "source": [ + "# MIRI Photometry Comparison of Photutils vs. Space_Phot\n", + "\n", + "**Author**: Ori Fox
\n", + "\n", + "**Submitted**: August, 2023
\n", + "**Updated**: November, 2023
\n", + "\n", + "**Use case**: A comparison of the photometry produced by the two other notebooks in this folder.
\n", + "**Data**: MIRI Data PID 1028 (Calibration Program; Single Star Visit 006 A5V dwarf 2MASSJ17430448+6655015) and MIRI Data PID 1171 (LMC; Multiple Stars).
\n", + "**Tools**: photutils, space_phot drizzlepac, jupyter
\n", + "**Cross-Instrument**: NIRCam, MIRI.
\n", + "**Documentation**: This notebook is part of a STScI's larger post-pipeline Data Analysis Tools Ecosystem and can be downloaded directly from the JDAT Notebook Github directory.
\n", + "**Pipeline Version**: JWST Pipeline
\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "18ff6a53-9e2d-4bef-aa2b-e2abd0758ec7", + "metadata": {}, + "outputs": [], + "source": [ + "import pandas as pd\n", + "import astropy.units as u\n", + "from astropy.coordinates import SkyCoord\n", + "from pandas import DataFrame" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "99887406-cc2a-4e62-8b8e-38de616e643d", + "metadata": {}, + "outputs": [], + "source": [ + "# Read space_phot data\n", + "\n", + "sphot = pd.read_csv('miri_photometry_space_phot_lvl2.txt') \n", + "sphot" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "672f0fde-cd6e-4e85-9518-a07b53f5f581", + "metadata": {}, + "outputs": [], + "source": [ + "# Read photutils data\n", + "\n", + "phot = pd.read_csv('miri_photometry_photutils.txt') \n", + "phot" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "1bf4436f-47ad-4f59-9992-1ede269724cd", + "metadata": {}, + "outputs": [], + "source": [ + "sphot['dec'].to_numpy()*u.deg" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "4a88bf03-50a0-4c9b-bc1d-1b9d93032ed0", + "metadata": {}, + "outputs": [], + "source": [ + "#Find closest neighbors\n", + "\n", + "c1 = SkyCoord(sphot['ra'].to_numpy()*u.deg, sphot['dec'].to_numpy()*u.deg, frame='icrs')\n", + "c2 = SkyCoord(phot['RA'].to_numpy()*u.deg, phot['DEC'].to_numpy()*u.deg, frame='icrs')\n", + "idx, d2d, d3d = c2.match_to_catalog_3d(c1)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "b32e3c67-397a-4515-94d6-afd8da0b850c", + "metadata": {}, + "outputs": [], + "source": [ + "# Calculate Delta Mag\n", + "\n", + "delta_mag = sphot['mag'][idx].to_numpy() - phot['Mag']" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "34f6ed58-d8dc-4af8-8f2a-d7498e118fd1", + "metadata": {}, + "outputs": [], + "source": [ + "# Cross match Catalogs and Look at Delta Mags\n", + "\n", + "df = DataFrame({\"photutils_skycoord\": c1[idx], \"spacephot_skycoord\": c2, \"skycoord_separation\": d3d, \"photutils_mag\": phot['Mag'], \"spacephot_mag\": sphot['mag'][idx].to_numpy(), \"delta_mag\": delta_mag})\n", + "df" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "8d0b5394-c11b-4fa0-922b-567225c7dc65", + "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/psf_photometry/MIRI/miri_photutils.ipynb b/notebooks/psf_photometry/MIRI/miri_photutils.ipynb new file mode 100644 index 000000000..54bf269ac --- /dev/null +++ b/notebooks/psf_photometry/MIRI/miri_photutils.ipynb @@ -0,0 +1,1066 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "b58d6954", + "metadata": {}, + "source": [ + "# MIRI PSF Photometry with Photutils\n", + "\n", + "**Author**: Ori Fox
\n", + "\n", + "**Submitted**: November, 2023
\n", + "**Updated**: November, 2023
\n", + "\n", + "**Use case**: PSF Photometry using [Photutils](https://photutils.readthedocs.io/en/stable/). The purpose here is to illustrate the workflow and runtime for using Photutils in a variety of use cases.\n", + "\n", + "Generally, PSF photometry for data from a space telescope is most accurately performed on pre-mosaiced data. The mosaic process changes the inherent PSF, blurring it both due to resampling and mixing PSFs at different detector positions and rotations. Additionally, accurate theoretical PSF models (e.g., from [WebbPSF](https://webbpsf.readthedocs.io/en/latest/)) are not available for mosaiced data. While an empirical PSF could be constructed (e.g., using Photutils [ePSFBuilder](https://photutils.readthedocs.io/en/latest/epsf.html)) for mosaiced data, the results will generally not be as accurate as performing PSF photometry on the pre-mosaiced data.\n", + "\n", + "**NOTE:** A companion notebook exists that illustrates how to use perform PSF photometry on both Level 2 and Level 3 data using a new software program called space_phot.
\n", + "**Data**: MIRI Data PID 1028 (Calibration Program; Single Star Visit 006 A5V dwarf 2MASSJ17430448+6655015) and MIRI Data PID 1171 (LMC; Multiple Stars).
\n", + "**Tools**: photutils, webbpsf, jwst
\n", + "**Cross-Instrument**: 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).
" + ] + }, + { + "cell_type": "markdown", + "id": "88c61bcf-1c4d-407a-b80c-aa13a01fd746", + "metadata": { + "tags": [] + }, + "source": [ + "## Table of contents\n", + "1. [Introduction](#intro)
\n", + " 1.1 [Python Imports](#imports)
\n", + " 1.2 [Set up WebbPSF and Synphot](#setup)
\n", + "2. [Download JWST MIRI Data](#data)
\n", + "3. [Single Bright Object](#bso)
\n", + " 3.1 [Single Level 2 File](#bso2)
\n", + " 3.2 [Generate empirical PSF grid for MIRI F770W using WebbPSF](#bso3)
\n", + " 3.3 [PSF Photometry](#bso4)
\n", + "4. [Faint/Upper Limit, Single Object](#fso)
\n", + " 4.1 [Multiple, Level2 Files](#fso2)
\n", + "5. [Stellar Field (LMC)](#lmc)
\n", + " 5.1 [Multiple Stars, Single Level 2 File](#lmc2)
\n", + " 5.2 [Generate empirical PSF grid for MIRI F560W using WebbPSF](#grid2)
\n", + " 5.3 [PSF Photometry](#lmc3)
" + ] + }, + { + "cell_type": "markdown", + "id": "4f572688", + "metadata": {}, + "source": [ + "# 1. Introduction " + ] + }, + { + "cell_type": "markdown", + "id": "95891849", + "metadata": {}, + "source": [ + "**GOALS**:
\n", + "\n", + "Perform PSF photometry on JWST MIRI images with the [Photutils PSF Photometry tools](https://photutils.readthedocs.io/en/latest/psf.html) using a grid of empirical PSF models from WebbPSF.\n", + "\n", + "\n", + "The notebook shows how to:
\n", + "\n", + "* generate a [grid of empirical PSF models](https://webbpsf.readthedocs.io/en/latest/psf_grids.html) from WebbPSF
\n", + "* perform PSF photometry on the image using the [PSFPhotometry class](https://photutils.readthedocs.io/en/latest/api/photutils.psf.PSFPhotometry.html#photutils.psf.PSFPhotometry)
\n", + "\n", + "**Data**:
\n", + "\n", + "MIRI Data PID 1028 (Calibration Program), F770W
\n", + "MIRI Data PID 1171 (LMC), F560W/F770W" + ] + }, + { + "cell_type": "markdown", + "id": "5e534877-5c31-4020-9263-4f234f19e1cd", + "metadata": {}, + "source": [ + "## 1.1 Python Imports " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "1ca81097-08d5-470c-942a-d8e7e8fd4479", + "metadata": {}, + "outputs": [], + "source": [ + "import glob\n", + "import os\n", + "import shutil\n", + "import tarfile\n", + "from urllib import request\n", + "\n", + "import astropy.units as u\n", + "import matplotlib.pyplot as plt\n", + "import numpy as np\n", + "import webbpsf\n", + "from astropy.coordinates import SkyCoord\n", + "from astropy.io import fits\n", + "from astropy.nddata import extract_array\n", + "from astropy.table import Table\n", + "from astropy.visualization import simple_norm\n", + "from astroquery.mast import Observations\n", + "from jwst.datamodels import ImageModel\n", + "from photutils.aperture import CircularAperture\n", + "from photutils.background import LocalBackground, MADStdBackgroundRMS, MMMBackground\n", + "from photutils.detection import DAOStarFinder\n", + "from photutils.psf import GriddedPSFModel, PSFPhotometry\n", + "\n", + "from pandas import DataFrame" + ] + }, + { + "cell_type": "markdown", + "id": "5b762602", + "metadata": {}, + "source": [ + "## 1.2 Download and Set up Required Data for WebbPSF and Synphot " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "8c50eace", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "# Set environmental variables\n", + "os.environ[\"WEBBPSF_PATH\"] = \"./webbpsf-data/webbpsf-data\"\n", + "os.environ[\"PYSYN_CDBS\"] = \"./grp/redcat/trds/\"\n", + "\n", + "# required webbpsf data\n", + "boxlink = 'https://stsci.box.com/shared/static/qxpiaxsjwo15ml6m4pkhtk36c9jgj70k.gz' \n", + "boxfile = './webbpsf-data/webbpsf-data-LATEST.tar.gz'\n", + "synphot_url = 'http://ssb.stsci.edu/trds/tarfiles/synphot5.tar.gz'\n", + "synphot_file = './synphot5.tar.gz'\n", + "\n", + "webbpsf_folder = './webbpsf-data'\n", + "synphot_folder = './grp'\n", + "\n", + "# Gather webbpsf files\n", + "if not os.path.exists(webbpsf_folder):\n", + " os.makedirs(webbpsf_folder)\n", + " request.urlretrieve(boxlink, boxfile)\n", + " gzf = tarfile.open(boxfile)\n", + " gzf.extractall(webbpsf_folder)\n", + "\n", + "# Gather synphot files\n", + "if not os.path.exists(synphot_folder):\n", + " os.makedirs(synphot_folder)\n", + " request.urlretrieve(synphot_url, synphot_file)\n", + " gzf = tarfile.open(synphot_file)\n", + " gzf.extractall('./')" + ] + }, + { + "cell_type": "markdown", + "id": "68f0b2d7-45a1-4511-858e-51425a50de00", + "metadata": {}, + "source": [ + "# 2. Download JWST MIRI Data " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "e6629878-a5d4-4e29-a56e-0f12271016a5", + "metadata": {}, + "outputs": [], + "source": [ + "# Download Proposal ID 1028 F770W data\n", + "\n", + "# Define source and destination directories\n", + "source_dir = 'mastDownload/JWST/'\n", + "destination_dir = 'mast/01028/'\n", + "\n", + "if os.path.isdir(destination_dir):\n", + " print(f'Data already downloaded to {os.path.abspath(destination_dir)}')\n", + "else:\n", + " # Query the MAST (Mikulski Archive for Space Telescopes) database for observations\n", + " # with proposal ID 1028 and the F770W filter\n", + " obs = Observations.query_criteria(proposal_id=1028, filters=['F770W'])\n", + " \n", + " # Get a list of products associated with the located observation\n", + " plist = Observations.get_product_list(obs)\n", + " \n", + " # Filter the product list to include only specific product subgroups\n", + " fplist = Observations.filter_products(plist, productSubGroupDescription=['CAL', 'I2D', 'ASN'])\n", + " \n", + " # Download the selected products from the MAST database\n", + " Observations.download_products(fplist)\n", + " \n", + " # Create the destination directory\n", + " os.makedirs(destination_dir)\n", + " \n", + " # Use glob to find all files matching the pattern\n", + " files_to_copy = glob.glob(os.path.join(source_dir, 'j*/jw01028*'))\n", + "\n", + " # Copy the matching files to the destination directory\n", + " for file_path in files_to_copy:\n", + " shutil.copy(file_path, destination_dir)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "b649b6f4-aa18-4760-a7d0-979b4e3caec2", + "metadata": {}, + "outputs": [], + "source": [ + "# Download Proposal ID 1171 F550W and F770W data\n", + "\n", + "# Define source and destination directories\n", + "source_dir = 'mastDownload/JWST/'\n", + "destination_dir = 'mast/01171/'\n", + "\n", + "if os.path.isdir(destination_dir):\n", + " print(f'Data already downloaded to {os.path.abspath(destination_dir)}')\n", + "else:\n", + " # Query the MAST (Mikulski Archive for Space Telescopes) database for observations\n", + " # with proposal ID 1171 and the F550W and F770W filters\n", + " obs = Observations.query_criteria(proposal_id=1171, filters=['F560W', 'F770W'])\n", + " \n", + " # Get a list of products associated with the located observation\n", + " plist = Observations.get_product_list(obs)\n", + " \n", + " # Filter the product list to include only specific product subgroups\n", + " fplist = Observations.filter_products(plist, productSubGroupDescription=['CAL', 'I2D', 'ASN'])\n", + " \n", + " # Download the selected products from the MAST database\n", + " Observations.download_products(fplist)\n", + " \n", + " # Create the destination directory\n", + " os.makedirs(destination_dir)\n", + " \n", + " # Use glob to find all files matching the pattern\n", + " files_to_copy = glob.glob(os.path.join(source_dir, 'j*/jw01171*'))\n", + " \n", + " # Copy the matching files to the destination directory\n", + " for file_path in files_to_copy:\n", + " shutil.copy(file_path, destination_dir)" + ] + }, + { + "cell_type": "markdown", + "id": "5611799d", + "metadata": {}, + "source": [ + "# 3. Single Bright Star " + ] + }, + { + "cell_type": "markdown", + "id": "6d052f0e-dcb4-4c2a-bc11-4b467dad07c2", + "metadata": {}, + "source": [ + "The purpose of this section is to illustrate how to perform PSF photometry on a single bright star. While aperture photometry is feasible in isolated cases, the user may find PSF photometry preferable in crowded fields or complicated backgrounds." + ] + }, + { + "cell_type": "markdown", + "id": "55c52f95", + "metadata": {}, + "source": [ + "## 3.1 Single Level 2 File " + ] + }, + { + "cell_type": "markdown", + "id": "058a14ad-be89-4d7e-934e-1e0a909319c8", + "metadata": {}, + "source": [ + "In this example, we fit a single, bright source in a single Level 2 images. For a collection of Level 2 images, we could fit each Level 2 image individually and then average the measured fluxes.\n", + "\n", + "Useful references:
\n", + "HST Documentation on PSF Photometry: https://www.stsci.edu/hst/instrumentation/wfc3/data-analysis/psf
\n", + "WFPC2 Stellar Photometry with HSTPHOT: https://ui.adsabs.harvard.edu/abs/2000PASP..112.1383D/abstract
\n", + "Photutils PSF Fitting Photometry: https://photutils.readthedocs.io/en/stable/psf.html" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "2af76a35-14ee-44b2-adb5-d92e66ddfafc", + "metadata": {}, + "outputs": [], + "source": [ + "# get the level 2 filenames\n", + "path = \"./mast/01028/\"\n", + "level2 = sorted(glob.glob(os.path.join(path, '*cal.fits')))\n", + "level2" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "ab0ac799-a832-40af-9cd4-b5b3934cfee4", + "metadata": {}, + "outputs": [], + "source": [ + "# display the first level-2 image\n", + "data = fits.getdata(level2[0])\n", + "norm = simple_norm(data, 'sqrt', percent=99)\n", + "\n", + "fig, ax = plt.subplots(figsize=(20, 12))\n", + "im = ax.imshow(data, origin='lower', norm=norm, cmap='gray')\n", + "clb = plt.colorbar(im, label='MJy/sr')\n", + "ax.set_xlabel('Pixels')\n", + "ax.set_ylabel('Pixels');" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "4ab947da-be0a-4cde-88eb-2d947adf2b81", + "metadata": {}, + "outputs": [], + "source": [ + "# Change all DQ flagged pixels to NaN.\n", + "# Here, we'll overwrite the original CAL file.\n", + "# Reference for JWST DQ Flag Definitions: https://jwst-pipeline.readthedocs.io/en/latest/jwst/references_general/references_general.html\n", + "# In this case, we choose all DQ > 10, but users are encouraged to choose their own values accordingly.\n", + "filename = level2[0]\n", + "with ImageModel(filename) as model:\n", + " model.data[model.dq >= 10] = np.nan\n", + " model.save(filename)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "87673dd8-1bc6-4663-bebc-8e2434974ceb", + "metadata": {}, + "outputs": [], + "source": [ + "# Re-display the image\n", + "data = fits.getdata(level2[0])\n", + "norm = simple_norm(data, 'sqrt', percent=99)\n", + "\n", + "fig, ax = plt.subplots(figsize=(20, 12))\n", + "im = ax.imshow(data, origin='lower', norm=norm, cmap='gray')\n", + "clb = plt.colorbar(im, label='MJy/sr')\n", + "ax.set_xlabel('Pixels')\n", + "ax.set_ylabel('Pixels');" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "2f138285-bfc5-4f20-854b-fe5d4530e99b", + "metadata": {}, + "outputs": [], + "source": [ + "# Zoom in to see the source. In this case, our source is from MIRI Program ID #1028, a Calibration Program.\n", + "# We are using Visit 006, which targets the A5V dwarf 2MASSJ17430448+6655015\n", + "# Reference Link: http://simbad.cds.unistra.fr/simbad/sim-basic?Ident=2MASSJ17430448%2B6655015&submit=SIMBAD+search\n", + "source_location = SkyCoord('17:43:04.4879', '+66:55:01.837', unit=(u.hourangle, u.deg))\n", + "with ImageModel(filename) as model:\n", + " x, y = model.meta.wcs.world_to_pixel(source_location)\n", + "\n", + "cutout = extract_array(data, (21, 21), (y, x))\n", + "\n", + "fig, ax = plt.subplots(figsize=(8, 8))\n", + "norm2 = simple_norm(cutout, 'log', percent=99)\n", + "im = ax.imshow(cutout, origin='lower', norm=norm2, cmap='gray')\n", + "clb = plt.colorbar(im, label='MJy/sr', shrink=0.8)\n", + "ax.set_title('PID1028, Obs006')\n", + "\n", + "ax.set_xlabel('Pixels')\n", + "ax.set_ylabel('Pixels');" + ] + }, + { + "cell_type": "markdown", + "id": "b24288b6-c6c7-433f-866f-126eea4a9ab3", + "metadata": { + "execution": { + "iopub.execute_input": "2024-02-10T00:33:17.900216Z", + "iopub.status.busy": "2024-02-10T00:33:17.899831Z", + "iopub.status.idle": "2024-02-10T00:33:17.903626Z", + "shell.execute_reply": "2024-02-10T00:33:17.902732Z", + "shell.execute_reply.started": "2024-02-10T00:33:17.900187Z" + } + }, + "source": [ + "## 3.2 Generate empirical PSF grid for MIRI F770W using WebbPSF " + ] + }, + { + "cell_type": "markdown", + "id": "f7472c4e-f418-4f39-b8b6-5d824c1c4c65", + "metadata": {}, + "source": [ + "Let's now use WebbPSF to generate an empirical grid of ePSF models for MIRI F770W.\n", + "The output will be a Photutils [GriddedPSFModel](https://photutils.readthedocs.io/en/latest/api/photutils.psf.GriddedPSFModel.html#photutils.psf.GriddedPSFModel) containing a 2x2 grid of detector-position-dependent empirical PSFs, each oversampled by a factor of 4. Note that we save the PSF grid to a FITS file (via `save=True`) called `miri_mirim_f770w_fovp101_samp4_npsf4.fits`. To save time in future runs, we load this FITS file directly into a `GriddedPSFModel` object:" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "19b97fc8-35b4-460f-b154-72bdd4c45f66", + "metadata": {}, + "outputs": [], + "source": [ + "psfgrid_filename = 'miri_mirim_f770w_fovp101_samp4_npsf4.fits'\n", + "\n", + "if not os.path.exists(psfgrid_filename):\n", + " miri = webbpsf.MIRI()\n", + " miri.filter = 'F770W'\n", + " psf_model = miri.psf_grid(num_psfs=4, all_detectors=True, verbose=True, save=True)\n", + "else:\n", + " psf_model = GriddedPSFModel.read(psfgrid_filename)\n", + "\n", + "psf_model" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "7d1640dc-6c64-4fe7-a208-ddbcdc529512", + "metadata": {}, + "outputs": [], + "source": [ + "# display the PSF grid\n", + "psf_model.plot_grid();" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "f7294300-6ba8-418b-a61c-4194af270de3", + "metadata": {}, + "outputs": [], + "source": [ + "# display the PSF grid deltas from the mean ePSF\n", + "psf_model.plot_grid(deltas=True);" + ] + }, + { + "cell_type": "markdown", + "id": "bcab882b-16f6-4aef-9f4f-b3af8b3b1b14", + "metadata": {}, + "source": [ + "## 3.3 PSF Photometry " + ] + }, + { + "cell_type": "markdown", + "id": "b34c6e74-1a52-4874-afc1-5a444662e966", + "metadata": {}, + "source": [ + "Now let's use our gridded PSF model to perform PSF photometry." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "f42d294f-a5db-4a2a-be4a-e40a26374207", + "metadata": {}, + "outputs": [], + "source": [ + "# load data and convert units from MJy/sr to uJy\n", + "with ImageModel(filename) as model:\n", + " unit = u.Unit(model.meta.bunit_data)\n", + " data = model.data << unit\n", + " error = model.err << unit\n", + "\n", + " # use pixel area map because of geometric distortion in level-2 data \n", + " pixel_area = model.area * model.meta.photometry.pixelarea_steradians * u.sr\n", + " data *= pixel_area\n", + " error *= pixel_area\n", + " \n", + " data = data.to(u.uJy)\n", + " error = error.to(u.uJy)\n", + "\n", + "data.unit, error.unit" + ] + }, + { + "cell_type": "markdown", + "id": "318ab078-ec51-4b71-9641-dd1fe566138a", + "metadata": {}, + "source": [ + "To perform photometry on a single source we can input a Table containing its (x, y) position." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "6efc06f3-445d-4139-a404-fff43ad15804", + "metadata": {}, + "outputs": [], + "source": [ + "init_params = Table()\n", + "init_params['x'] = [x]\n", + "init_params['y'] = [y]\n", + "init_params" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "61424bfb-7c76-489b-9ccf-efc9c324472d", + "metadata": {}, + "outputs": [], + "source": [ + "# we turn off the finder because we input the source position\n", + "fit_shape = 5\n", + "localbkg_estimator = LocalBackground(5, 10, bkg_estimator=MMMBackground())\n", + "psfphot = PSFPhotometry(psf_model, fit_shape, finder=None, aperture_radius=fit_shape, \n", + " localbkg_estimator=localbkg_estimator, progress_bar=True)\n", + "phot = psfphot(data, error=error, init_params=init_params)\n", + "phot" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "62945e0b-6eb7-487f-a7b6-8e3c718b3c26", + "metadata": {}, + "outputs": [], + "source": [ + "# convert fit flux from uJy to ABmag\n", + "flux = phot['flux_fit']\n", + "flux_err = phot['flux_err']\n", + "mag = phot['flux_fit'].to(u.ABmag)\n", + "magerr = 2.5 * np.log10(1.0 + (flux_err / flux))\n", + "magerr = magerr.value * u.ABmag\n", + "mag, magerr" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "6c8baf43-51bc-4ff5-97e9-10db55ab6128", + "metadata": {}, + "outputs": [], + "source": [ + "fig, ax = plt.subplots(ncols=3, figsize=(12, 4))\n", + "\n", + "shape = (21, 21)\n", + "cutout1 = extract_array(data.value, shape, (y, x))\n", + "norm = simple_norm(cutout1, 'log', percent=98)\n", + "im1 = ax[0].imshow(cutout1, origin='lower', norm=norm)\n", + "ax[0].set_title(r'Data ($\\mu$Jy)')\n", + "plt.colorbar(im1, shrink=0.7)\n", + "\n", + "model = psfphot.make_model_image(data.shape, shape)\n", + "cutout2 = extract_array(model, shape, (y, x))\n", + "im2 = ax[1].imshow(cutout2, origin='lower', norm=norm)\n", + "ax[1].set_title('Fit PSF Model')\n", + "plt.colorbar(im2, shrink=0.7)\n", + "\n", + "resid = psfphot.make_residual_image(data.value, shape)\n", + "cutout3 = extract_array(resid, shape, (y, x))\n", + "norm3 = simple_norm(cutout3, 'sqrt', percent=99)\n", + "im3 = ax[2].imshow(cutout3, origin='lower', norm=norm3)\n", + "ax[2].set_title('Residual')\n", + "plt.colorbar(im3, shrink=0.7)" + ] + }, + { + "cell_type": "markdown", + "id": "5b5f0ad5-b59e-4eff-8687-f6a2199d8bd9", + "metadata": {}, + "source": [ + "# 4. Faint/Upper Limit, Single Object " + ] + }, + { + "cell_type": "markdown", + "id": "1dc60da1-0f4d-4e5e-a109-f7352dfd0fdc", + "metadata": {}, + "source": [ + "The purpose of this section is to illustrate how to calculate an upper limit at a fixed (x, y) position using forced PSF photometry a blank part of the sky.\n", + "\n", + "We'll use the same data as Section 3." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "9311116e-630b-4a2e-8fc6-be75f6de61fd", + "metadata": {}, + "outputs": [], + "source": [ + "# load data and convert units from MJy/sr to uJy\n", + "with ImageModel(filename) as model:\n", + " unit = u.Unit(model.meta.bunit_data)\n", + " data = model.data << unit\n", + " error = model.err << unit\n", + " \n", + " pixel_area = pixel_area = model.meta.photometry.pixelarea_steradians * u.sr\n", + " data *= pixel_area\n", + " error *= pixel_area\n", + " \n", + " data = data.to(u.uJy)\n", + " error = error.to(u.uJy)\n", + "\n", + "source_location = SkyCoord('17:43:00.0332', '+66:54:42.677', unit=(u.hourangle, u.deg))\n", + "with ImageModel(filename) as model:\n", + " x, y = model.meta.wcs.world_to_pixel(source_location)\n", + "\n", + "cutout = extract_array(data.value, (21, 21), (y, x))\n", + "\n", + "fig, ax = plt.subplots()\n", + "norm = simple_norm(cutout, 'sqrt', percent=95)\n", + "im = ax.imshow(cutout, origin='lower', norm=norm, cmap='gray')\n", + "clb = plt.colorbar(im, label=r'$\\mu$Jy')\n", + "ax.set_title('PID1028, Obs006')\n", + "\n", + "ax.set_xlabel('Pixels')\n", + "ax.set_ylabel('Pixels');" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "63ac28fd-ad3b-479b-8561-ca5913925aba", + "metadata": {}, + "outputs": [], + "source": [ + "# to perform forced photometry, we set the (x, y) source position\n", + "# AND we fix the PSF model position so that it does not vary in the fit\n", + "# (only flux will be fit)\n", + "init_params = Table()\n", + "init_params['x'] = [x]\n", + "init_params['y'] = [y]\n", + "\n", + "# This requires photutils 1.11.0\n", + "psf_model_forced = psf_model.copy()\n", + "psf_model_forced.x_0.fixed = True\n", + "psf_model_forced.y_0.fixed = True\n", + "psf_model_forced.fixed" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "a13930cb-fd2a-447e-a68f-357eb7014ad3", + "metadata": {}, + "outputs": [], + "source": [ + "fit_shape = 5\n", + "localbkg_estimator = LocalBackground(5, 10, bkg_estimator=MMMBackground())\n", + "psfphot = PSFPhotometry(psf_model_forced, fit_shape, finder=None, aperture_radius=fit_shape, \n", + " localbkg_estimator=localbkg_estimator, progress_bar=True)\n", + "\n", + "phot = psfphot(data, error=error, init_params=init_params)\n", + "phot" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "3206847e-2e2e-4bc3-a51e-3e5877f87b28", + "metadata": {}, + "outputs": [], + "source": [ + "# To calculate upper limit, multiply the flux_err by your desired sigma\n", + "sigma = 3.0\n", + "limit = sigma * phot['flux_err']\n", + "limit.to(u.ABmag)" + ] + }, + { + "cell_type": "markdown", + "id": "4e4996d2-6274-473d-b13c-f3848c27ad78", + "metadata": {}, + "source": [ + "## Note: you can go significantly deeper with the Level 3 combined data product" + ] + }, + { + "cell_type": "markdown", + "id": "9a969717-bbef-40b9-ac9b-f83dec99dc09", + "metadata": {}, + "source": [ + "# 5. Stellar Field (LMC) " + ] + }, + { + "cell_type": "markdown", + "id": "da877310-fd47-41d6-afea-7fa725a546af", + "metadata": {}, + "source": [ + "In this case, we are going to do the same steps as in Section 3, but for multiple stars. The purpose is to illustrate the workflow and runtime for using Photutils on a large number of stars." + ] + }, + { + "cell_type": "markdown", + "id": "32bdafe6-db19-4080-9587-b9785c2f7fa7", + "metadata": {}, + "source": [ + "## 5.1 Multiple Stars, Single Level 2 File " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "72aabedd-8d80-4c6e-8d3d-9287577665bf", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "# Find stars in Level 3 File\n", + "path = './mast/01171/'\n", + "level3 = os.path.join(path, 'jw01171-o004_t001_miri_f560w_i2d.fits')\n", + "level3" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "4eb6155b-5f1a-4f92-9c34-5cd0876e439d", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "# Get rough estimate of background (there are better ways to do background subtraction)\n", + "bkgrms = MADStdBackgroundRMS()\n", + "mmm_bkg = MMMBackground()\n", + "\n", + "with ImageModel(level3) as model:\n", + " wcs_l3 = model.meta.wcs\n", + " std = bkgrms(model.data)\n", + " bkg = mmm_bkg(model.data)\n", + " data_bkgsub = model.data.copy()\n", + " data_bkgsub -= bkg \n", + "\n", + "# Find stars\n", + "# F560W FWHM = 1.882 pix\n", + "fwhm_psf = 1.882\n", + "threshold = 5.0\n", + "daofind = DAOStarFinder(threshold=threshold * std, fwhm=fwhm_psf, exclude_border=True, min_separation=10)\n", + "found_stars = daofind(data_bkgsub)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "59337919-ee59-4f09-a716-1f978ab5c775", + "metadata": {}, + "outputs": [], + "source": [ + "found_stars.pprint_all(max_lines=10)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "73aba802", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "# plot the found stars\n", + "norm = simple_norm(data_bkgsub, 'sqrt', percent=99)\n", + "fig, ax = plt.subplots(figsize=(10, 10))\n", + "ax.imshow(data_bkgsub, origin='lower', norm=norm)\n", + "\n", + "xypos = zip(found_stars['xcentroid'], found_stars['ycentroid'])\n", + "aper = CircularAperture(xypos, r=10)\n", + "aper.plot(ax, color='red');" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "84411432-3270-47de-a546-a81935c19c5c", + "metadata": {}, + "outputs": [], + "source": [ + "fig, ax = plt.subplots(nrows=2, figsize=(10, 8))\n", + "\n", + "ax[0].scatter(found_stars['mag'], found_stars['sharpness'], s=10, color='k')\n", + "ax[0].set_xlabel('mag')\n", + "ax[0].set_ylabel('sharpness')\n", + "\n", + "ax[1].scatter(found_stars['mag'], found_stars['roundness2'], s=10, color='k')\n", + "ax[1].set_xlabel('mag')\n", + "ax[1].set_ylabel('roundness')\n", + "\n", + "mag0 = -3.0\n", + "mag1 = -5.0\n", + "for ax_ in ax:\n", + " ax_.axvline(mag0, color='red', linestyle='dashed')\n", + " ax_.axvline(mag1, color='red', linestyle='dashed')\n", + "\n", + "sh0 = 0.40\n", + "sh1 = 0.82\n", + "ax[0].axhline(sh0, color='red', linestyle='dashed')\n", + "ax[0].axhline(sh1, color='red', linestyle='dashed')\n", + "\n", + "rnd0 = -0.40\n", + "rnd1 = 0.40\n", + "ax[1].axhline(rnd0, color='red', linestyle='dashed')\n", + "ax[1].axhline(rnd1, color='red', linestyle='dashed');" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "04f2dbae-1b1c-4f67-902f-eb91f73b2d6d", + "metadata": {}, + "outputs": [], + "source": [ + "mask = ((found_stars['mag'] < mag0) & (found_stars['mag'] > mag1) & (found_stars['roundness2'] > rnd0)\n", + " & (found_stars['roundness2'] < rnd1) & (found_stars['sharpness'] > sh0) \n", + " & (found_stars['sharpness'] < sh1) & (found_stars['xcentroid'] > 100) & (found_stars['xcentroid'] < 700)\n", + " & (found_stars['ycentroid'] > 100) & (found_stars['ycentroid'] < 700))\n", + "\n", + "found_stars_sel = found_stars[mask]\n", + "\n", + "print('Number of stars found originally:', len(found_stars))\n", + "print('Number of stars in final selection:', len(found_stars_sel))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "948f0a96-8bdc-4c94-95db-8edf513e4094", + "metadata": {}, + "outputs": [], + "source": [ + "# plot the selected stars\n", + "norm = simple_norm(data_bkgsub, 'sqrt', percent=99)\n", + "fig, ax = plt.subplots(figsize=(10, 10))\n", + "ax.imshow(data_bkgsub, origin='lower', norm=norm)\n", + "\n", + "xypos = zip(found_stars_sel['xcentroid'], found_stars_sel['ycentroid'])\n", + "aper = CircularAperture(xypos, r=10)\n", + "aper.plot(ax, color='red');" + ] + }, + { + "cell_type": "markdown", + "id": "a221246a-9acb-4c79-9831-2f6cb72bb90a", + "metadata": { + "execution": { + "iopub.execute_input": "2024-02-10T02:52:51.761591Z", + "iopub.status.busy": "2024-02-10T02:52:51.761230Z", + "iopub.status.idle": "2024-02-10T02:52:51.765256Z", + "shell.execute_reply": "2024-02-10T02:52:51.764441Z", + "shell.execute_reply.started": "2024-02-10T02:52:51.761562Z" + } + }, + "source": [ + "## 5.2 Generate empirical PSF grid for MIRI F560W using WebbPSF " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "b5fdcf15-7f72-40b8-b1aa-2b5fb0375ed8", + "metadata": {}, + "outputs": [], + "source": [ + "psfgrid_filename = 'miri_mirim_f560w_fovp101_samp4_npsf4.fits'\n", + "\n", + "if not os.path.exists(psfgrid_filename):\n", + " miri = webbpsf.MIRI()\n", + " miri.filter = 'F560W'\n", + " psf_model = miri.psf_grid(num_psfs=4, all_detectors=True, verbose=True, save=True)\n", + "else:\n", + " psf_model = GriddedPSFModel.read(psfgrid_filename)\n", + "\n", + "psf_model" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "29fa22f3-0331-4966-9bb4-c75f325011c8", + "metadata": {}, + "outputs": [], + "source": [ + "psf_model.plot_grid();" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "b2350526-2bba-4322-9092-88b6acdc117b", + "metadata": {}, + "outputs": [], + "source": [ + "# get the level 2 image\n", + "# here, we'll use the PID 1171 files\n", + "path = \"./mast/01171/\"\n", + "level2 = sorted(glob.glob(os.path.join(path, 'jw01171004*cal.fits')))\n", + "filename = level2[0]\n", + "print(filename)\n", + "\n", + "# load data and convert units from MJy/sr to uJy\n", + "with ImageModel(filename) as model:\n", + " unit = u.Unit(model.meta.bunit_data)\n", + " model.data[model.dq >= 10] = np.nan\n", + " data = model.data << unit\n", + " error = model.err << unit\n", + " \n", + " pixel_area = pixel_area = model.meta.photometry.pixelarea_steradians * u.sr\n", + " data *= pixel_area\n", + " error *= pixel_area\n", + " \n", + " data = data.to(u.uJy)\n", + " error = error.to(u.uJy)\n", + "\n", + " wcs = model.meta.wcs\n", + "\n", + "data.unit, error.unit" + ] + }, + { + "cell_type": "markdown", + "id": "9aa77cc2-26e2-4ac3-8384-90b9f8567f22", + "metadata": {}, + "source": [ + "## 5.3 PSF Photometry " + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "3b416d45-1081-4a16-b189-cc0e5e9d4e4e", + "metadata": {}, + "outputs": [], + "source": [ + "# translate (x, y) positions from the level 3 image to the level 2 image\n", + "xc = found_stars_sel['xcentroid']\n", + "yc = found_stars_sel['ycentroid']\n", + "sc = wcs_l3.pixel_to_world(xc, yc)\n", + "\n", + "x, y = wcs.world_to_pixel(sc)\n", + "init_params = Table()\n", + "init_params['x'] = x\n", + "init_params['y'] = y\n", + "\n", + "# we need to remove stars in the masked region of\n", + "# the level-2 data\n", + "mask = x > 400\n", + "init_params = init_params[mask]" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "f17e393a-9f0b-4b83-83b9-4bb8497dd835", + "metadata": {}, + "outputs": [], + "source": [ + "mask" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "511183e0-0adc-4966-b7f4-fcd912b054e9", + "metadata": {}, + "outputs": [], + "source": [ + "# plot the selected stars\n", + "norm = simple_norm(data.value, 'sqrt', percent=99)\n", + "fig, ax = plt.subplots(figsize=(10, 10))\n", + "ax.imshow(data.value, origin='lower', norm=norm)\n", + "\n", + "xypos = zip(init_params['x'], init_params['y'])\n", + "aper = CircularAperture(xypos, r=10)\n", + "aper.plot(ax, color='red');" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "84812421-f447-4b75-9fe4-e377fba54b17", + "metadata": {}, + "outputs": [], + "source": [ + "fit_shape = 5\n", + "localbkg_estimator = LocalBackground(5, 10, bkg_estimator=MMMBackground())\n", + "psfphot = PSFPhotometry(psf_model, fit_shape, finder=None, aperture_radius=fit_shape, \n", + " localbkg_estimator=localbkg_estimator, progress_bar=True)\n", + "phot = psfphot(data, error=error, init_params=init_params)\n", + "phot" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "8d5899cf-8956-4b3c-941a-67b0fedc0c51", + "metadata": {}, + "outputs": [], + "source": [ + "# convert fit flux from uJy to ABmag\n", + "flux = phot['flux_fit']\n", + "flux_err = phot['flux_err']\n", + "mag = phot['flux_fit'].to(u.ABmag)\n", + "magerr = 2.5 * np.log10(1.0 + (flux_err / flux))\n", + "magerr = magerr.value * u.ABmag\n", + "mag, magerr" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "a5ca496b-9826-450e-982c-f671c009c23f", + "metadata": {}, + "outputs": [], + "source": [ + "# Write to File\n", + "\n", + "df = DataFrame({\"RA\": sc.ra.deg[mask], \"DEC\": sc.dec.deg[mask], \"Mag\": mag.value, \"Mag Err\": magerr.value})\n", + "df.to_csv('miri_photometry_photutils.txt', index=False) " + ] + }, + { + "cell_type": "markdown", + "id": "5630029f-31d1-42cd-8454-225e86cabc48", + "metadata": {}, + "source": [ + "
" + ] + }, + { + "cell_type": "markdown", + "id": "843b5201-6f57-46f0-9da0-b738714178d3", + "metadata": {}, + "source": [ + "\"Space" + ] + } + ], + "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" + }, + "toc-showcode": false + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/notebooks/psf_photometry/MIRI/miri_spacephot.ipynb b/notebooks/psf_photometry/MIRI/miri_spacephot.ipynb new file mode 100644 index 000000000..de881cd53 --- /dev/null +++ b/notebooks/psf_photometry/MIRI/miri_spacephot.ipynb @@ -0,0 +1,2256 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "b58d6954", + "metadata": {}, + "source": [ + "# MIRI PSF Photometry With Space_Phot\n", + "\n", + "**Author**: Ori Fox
\n", + "\n", + "**Submitted**: August, 2023
\n", + "**Updated**: November, 2023
\n", + "\n", + "**Use case**: PSF Photometry on Level3 data using dedicated package Space_Phot (https://github.com/jpierel14/space_phot). Space_phot is built on Astropy's Photutils package. Unlike photutils, space_phot can be used on Level3 data. This is because has built in functionality that can generate a resampled Level3 PSF at a given detector position. Such use cases are particularly useful for faint, point source targets or deep upper limits where observers need to take advantage of the combined image stack. For a large number of bright sources, users may find space_phot to be too slow and should consider other packages, such as DOLPHOT and/or Photutils. **NOTE:** A companion notebook exists that illustrates how to use Photutils for the same Level2 data set.
\n", + "**Important Note**: When not to use. Due to the sensitivity of the space_phot parameters, this tool is not meant to be used for a large sample of stars (i.e., Section 5 below). If a user would like to use space_phot on more than one source, they should carefully construct a table of parameters that are carefully refined for each source.\n", + "**Data**: MIRI Data PID 1028 (Calibration Program; Single Star Visit 006 A5V dwarf 2MASSJ17430448+6655015) and MIRI Data PID 1171 (LMC; Multiple Stars).
\n", + "**Tools**: photutils, space_phot drizzlepac, jupyter
\n", + "**Cross-Instrument**: NIRCam, MIRI.
\n", + "**Documentation**: This notebook is part of a STScI's larger post-pipeline Data Analysis Tools Ecosystem and can be downloaded directly from the JDAT Notebook Github directory.
\n", + "**Pipeline Version**: JWST Pipeline
\n" + ] + }, + { + "cell_type": "markdown", + "id": "88c61bcf-1c4d-407a-b80c-aa13a01fd746", + "metadata": { + "tags": [] + }, + "source": [ + "## Table of contents\n", + "1. [Introduction](#intro)
\n", + " 1.1 [Setup](#webbpsf)
\n", + " 1.2 [Python imports](#py_imports)
\n", + "2. [Download Data](#data)
\n", + "3. [Bright, Single Object](#bso)
\n", + " 3.1 [Multiple, Level2 Files](#bso2)
\n", + " 3.2 [Single, Level3 Mosaicked File](#bso3)
\n", + "4. [Faint/Upper Limit, Single Object](#fso)
\n", + " 4.1 [Multiple, Level2 Files](#fso2)
\n", + " 4.2 [Single, Level3 Mosaicked File](#fso3)
\n", + "5. [Stellar Field (LMC)](#lmv)
\n", + " 5.1 [Multiple, Level2 Files](#lmc2)
\n", + " 5.2 [Single, Level3 Mosaicked File](#lmc3)
" + ] + }, + { + "cell_type": "markdown", + "id": "4f572688", + "metadata": {}, + "source": [ + "1.-Introduction \n", + "------------------" + ] + }, + { + "cell_type": "markdown", + "id": "95891849", + "metadata": {}, + "source": [ + "GOALS:
\n", + "\n", + "PSF Photometry can be obtained using:
\n", + "\n", + "* grid of PSF models from WebbPSF
\n", + "* single effective PSF (ePSF) NOT YET AVAILABLE
\n", + "* grid of effective PSF NOT YET AVAILABLE
\n", + "\n", + "The notebook shows:
\n", + "\n", + "* how to obtain the PSF model from WebbPSF (or build an ePSF)
\n", + "* how to perform PSF photometry on the image
\n", + "\n", + "**Data**:
\n", + "\n", + "MIRI Data PID 1028 (Calibration Program), F770W
\n", + "MIRI Data PID 1171 (LMC), F560W/F770W" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "991de542-26f4-4847-a6c8-1e6e1b39e757", + "metadata": {}, + "outputs": [], + "source": [ + "%matplotlib inline" + ] + }, + { + "cell_type": "markdown", + "id": "5b762602", + "metadata": {}, + "source": [ + "### 1.1-Setup WebbPSF and Synphot Directories ###" + ] + }, + { + "cell_type": "raw", + "id": "d7df5c79-da6b-4d50-bca3-23256e9afc80", + "metadata": {}, + "source": [ + "import space_phot\n", + "from importlib.metadata import version\n", + "version('space_phot')" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "ced803b2-d0f5-454f-bce2-6628aa905a4d", + "metadata": {}, + "outputs": [], + "source": [ + "%matplotlib inline" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "8c50eace", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "import sys,os,glob,shutil\n", + "import tarfile, urllib.request\n", + "\n", + "# Set environmental variables\n", + "os.environ[\"WEBBPSF_PATH\"] = \"./webbpsf-data/webbpsf-data\"\n", + "os.environ[\"PYSYN_CDBS\"] = \"./grp/redcat/trds/\"\n", + "\n", + "# WEBBPSF Data\n", + "boxlink = 'https://stsci.box.com/shared/static/qxpiaxsjwo15ml6m4pkhtk36c9jgj70k.gz' \n", + "boxfile = './webbpsf-data/webbpsf-data-1.0.0.tar.gz'\n", + "synphot_url = 'http://ssb.stsci.edu/trds/tarfiles/synphot5.tar.gz'\n", + "synphot_file = './synphot5.tar.gz'\n", + "\n", + "webbpsf_folder = './webbpsf-data'\n", + "synphot_folder = './grp'\n", + "\n", + "# Gather webbpsf files\n", + "psfExist = os.path.exists(webbpsf_folder)\n", + "if not psfExist:\n", + " os.makedirs(webbpsf_folder)\n", + " urllib.request.urlretrieve(boxlink, boxfile)\n", + " gzf = tarfile.open(boxfile)\n", + " gzf.extractall(webbpsf_folder)\n", + "\n", + "# Gather synphot files\n", + "synExist = os.path.exists(synphot_folder)\n", + "if not synExist:\n", + " os.makedirs(synphot_folder)\n", + " urllib.request.urlretrieve(synphot_url, synphot_file)\n", + " gzf = tarfile.open(synphot_file)\n", + " gzf.extractall('./')" + ] + }, + { + "cell_type": "markdown", + "id": "5e534877-5c31-4020-9263-4f234f19e1cd", + "metadata": {}, + "source": [ + "### 1.2-Python Imports ###" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "93e91ff7-ea71-4507-b38d-c2067cae3a8f", + "metadata": {}, + "outputs": [], + "source": [ + "from astropy.io import fits\n", + "from astropy.table import Table\n", + "from astropy.nddata import extract_array\n", + "from astropy.coordinates import SkyCoord\n", + "from astropy import wcs\n", + "from astropy.wcs.utils import skycoord_to_pixel\n", + "from astropy import units as u\n", + "import numpy as np\n", + "import matplotlib.pyplot as plt\n", + "from astroquery.mast import Observations\n", + "from astropy.visualization import (simple_norm,LinearStretch)\n", + "from mpl_toolkits.axes_grid1 import make_axes_locatable\n", + "import time\n", + "import math\n", + "from importlib.metadata import version\n", + "\n", + "import space_phot\n", + "\n", + "# JWST models\n", + "#\n", + "from jwst import datamodels, associations\n", + "from jwst.datamodels import ImageModel, dqflags\n", + "\n", + "# Background and PSF Functions\n", + "#\n", + "from photutils.background import MMMBackground, MADStdBackgroundRMS, Background2D, LocalBackground, MADStdBackgroundRMS\n", + "\n", + "from photutils.detection import DAOStarFinder\n", + "from photutils import EPSFBuilder, GriddedPSFModel\n", + "from photutils.psf import DAOGroup, extract_stars, IterativelySubtractedPSFPhotometry\n", + "\n", + "# Photutils library and tools\n", + "#\n", + "import photutils\n", + "from photutils.aperture import CircularAperture, CircularAnnulus, aperture_photometry\n", + "from photutils import Background2D, MedianBackground, ModeEstimatorBackground, MMMBackground" + ] + }, + { + "cell_type": "raw", + "id": "a65873e6-a9c7-45ba-99df-b8f6d9fc1faf", + "metadata": {}, + "source": [ + "print(np.log(10)*0.1/2.5)\n", + "print(np.log10(10**37.1))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "a5ff3213-2d53-4b6c-ab14-bd5f56643b71", + "metadata": {}, + "outputs": [], + "source": [ + "print(np.sqrt(5.8e41/(4.*np.pi*5.67e-5*650**4.)))\n", + "print(np.sqrt(3.7e40/(4.*np.pi*5.67e-5*300**4.)))\n", + "#5.8e41, 650\n", + "#3.7e40, 220" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "e6082768-3e91-4c83-8e1a-e98235e50bc4", + "metadata": {}, + "outputs": [], + "source": [ + "!pip install space_phot --upgrade" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "7bc12d37-4dc6-41e8-a1a4-756e143ea2eb", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "# Space_Phot Version Number for Reference \n", + "from importlib.metadata import version\n", + "version('space_phot')" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "52505408-5e55-4c82-857b-d6c6baf7c7c9", + "metadata": {}, + "outputs": [], + "source": [ + "version('poppy')" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "d7a88f4b-c40a-4dd0-88ef-4fe6638159e8", + "metadata": {}, + "outputs": [], + "source": [ + "version('jwst')" + ] + }, + { + "cell_type": "markdown", + "id": "68f0b2d7-45a1-4511-858e-51425a50de00", + "metadata": {}, + "source": [ + "2.-Download Data\n", + "------------------" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "e6629878-a5d4-4e29-a56e-0f12271016a5", + "metadata": {}, + "outputs": [], + "source": [ + "# Query the MAST (Mikulski Archive for Space Telescopes) database for observations\n", + "# with proposal ID 1028 and a specific filter 'F770W'\n", + "\n", + "obs = Observations.query_criteria(proposal_id=1028, filters=['F770W'])\n", + "\n", + "# Locate a specific observation by its unique observation ID 'jw01537-o024_t001_nircam_clear-f444w-sub160'\n", + "#row = obs[obs['obs_id']=='jw01028-o006_t001_miri_f770w']\n", + "\n", + "# Get a list of products associated with the located observation\n", + "plist = Observations.get_product_list(obs)\n", + "\n", + "# Filter the product list to include only specific product subgroups: 'RATE', 'CAL', 'I2D', and 'ASN'\n", + "fplist = Observations.filter_products(plist, productSubGroupDescription=['CAL', 'I2D', 'ASN'])\n", + "\n", + "# Download the selected products from the MAST database (UNCOMMENT TO DOWNLOAD)\n", + "Observations.download_products(fplist)\n", + "\n", + "# Define source and destination directories\n", + "source_dir = 'mastDownload/JWST/'\n", + "destination_dir = 'mast/01028/'\n", + "\n", + "# Create the destination directory if it doesn't exist\n", + "if not os.path.exists(destination_dir):\n", + " os.makedirs(destination_dir)\n", + "\n", + "# Use glob to find all files matching the pattern 'mastDownload/JWST/j*/jw01537*cal.fits'\n", + "files_to_copy = glob.glob(os.path.join(source_dir, 'j*/jw01028*'))\n", + "\n", + "# Copy the matching files to the destination directory\n", + "for file_path in files_to_copy:\n", + " shutil.copy(file_path, destination_dir)\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "b649b6f4-aa18-4760-a7d0-979b4e3caec2", + "metadata": {}, + "outputs": [], + "source": [ + "# Query the MAST (Mikulski Archive for Space Telescopes) database for observations\n", + "# with proposal ID 1171 and a specific filters 'F550W' and 'F770W'\n", + "\n", + "obs = Observations.query_criteria(proposal_id=1171, filters=['F560W','F770W'])\n", + "\n", + "# Locate a specific observation by its unique observation ID 'jw01537-o024_t001_nircam_clear-f444w-sub160'\n", + "#row = obs[obs['obs_id']=='jw01028-o006_t001_miri_f770w']\n", + "\n", + "# Get a list of products associated with the located observation\n", + "plist = Observations.get_product_list(obs)\n", + "\n", + "# Filter the product list to include only specific product subgroups: 'RATE', 'CAL', 'I2D', and 'ASN'\n", + "fplist = Observations.filter_products(plist, productSubGroupDescription=['CAL', 'I2D', 'ASN'])\n", + "fplist\n", + "\n", + "# Download the selected products from the MAST database (UNCOMMENT TO DOWNLOAD)\n", + "Observations.download_products(fplist)\n", + "\n", + "# Define source and destination directories\n", + "source_dir = 'mastDownload/JWST/'\n", + "destination_dir = 'mast/01171/'\n", + "\n", + "# Create the destination directory if it doesn't exist\n", + "if not os.path.exists(destination_dir):\n", + " os.makedirs(destination_dir)\n", + "\n", + "# Use glob to find all files matching the pattern 'mastDownload/JWST/j*/jw01537*cal.fits'\n", + "files_to_copy = glob.glob(os.path.join(source_dir, 'j*/jw01171*'))\n", + "\n", + "# Copy the matching files to the destination directory\n", + "for file_path in files_to_copy:\n", + " shutil.copy(file_path, destination_dir)" + ] + }, + { + "cell_type": "markdown", + "id": "5611799d", + "metadata": {}, + "source": [ + "3.-Bright, Single Object\n", + "------------------" + ] + }, + { + "cell_type": "markdown", + "id": "6d052f0e-dcb4-4c2a-bc11-4b467dad07c2", + "metadata": {}, + "source": [ + "The purpose of this section is to illustrate how to perform PSF photometry on a single, bright object. While aperture photometry is feasible in isolated cases, the user may find PSF photometry preferable in crowded fields or complicated backgrounds." + ] + }, + { + "cell_type": "markdown", + "id": "55c52f95", + "metadata": {}, + "source": [ + "### 3.1-Multiple, Level2 Files ###" + ] + }, + { + "cell_type": "markdown", + "id": "058a14ad-be89-4d7e-934e-1e0a909319c8", + "metadata": {}, + "source": [ + "Generally, PSF photometry for data from a space telescope is most accurately performed on pre-mosaiced data. In the case of HST, that corresponds to FLT files rather than DRZ. And in the case of JWST, this corresponds to Level2 files rather than Level3. The reason is that a mosaiced PSF changes the inherent PSF as a function of position on the detector so that there is no adequate model (theoretical or empirical) to use.
\n", + "\n", + "In this example, we aim to fit a source simultaneously across multiple Level 2 images. A more basic approach would be to fit each Level 2 file individually and then average together the measured fluxes. However, this approach more easily corrects for bad pixels or cosmic rays that are only in one image and allows for a more accurate photometric solution by reducing the number of free parameters per source.
\n", + "\n", + "Useful references:
\n", + "HST Documentation on PSF Photometry: https://www.stsci.edu/hst/instrumentation/wfc3/data-analysis/psf
\n", + "WFPC2 Stellar Photometry with HSTPHOT: https://ui.adsabs.harvard.edu/abs/2000PASP..112.1383D/abstract
\n", + "Space-Phot Documentation on Level2 Fitting: https://space-phot.readthedocs.io/en/latest/examples/plot_a_psf.html#jwst-images
" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "8f44a6cf", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "### Define Level 3 File\n", + "lvl3 = ['./mast/01028/jw01028-o006_t001_miri_f770w_i2d.fits']\n", + "#lvl3 = ['./mast/01028/obsnum06/jw01028-o006_t001_miri_f770w_i2d.fits']\n", + "lvl3" + ] + }, + { + "cell_type": "raw", + "id": "c33e5985", + "metadata": {}, + "source": [ + "asnfile = 'jw02079-o004_20230622t175524_spec2_00001_asn.json'\n", + "asn_data = json.load(open(asnfile))\n", + "print(asn_data['products'][0]['members'][0]['expname'])" + ] + }, + { + "cell_type": "raw", + "id": "ec9baac8", + "metadata": {}, + "source": [ + "from jwst.associations import load_asn\n", + "asn_data = load_asn(open(\"my_asn.json\"))\n", + "for member in asn_data['products'][0]['members']:\n", + " print(member['expname'])" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "d35e67aa", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "### Create Level 2 Data List from ASN files\n", + "\n", + "prefix = \"./mast/01028/\"\n", + "asn = glob.glob(prefix+'jw01028-o006_*_image3_00004_asn.json')\n", + "with open(asn[0],\"r\") as fi:\n", + " lvl2 = []\n", + " for ln in fi:\n", + " #print(ln)\n", + " if ln.startswith(' \"expname\":'):\n", + " x = ln[2:].split(':')\n", + " y = x[1].split('\"')\n", + " lvl2.append(prefix+y[1])\n", + "print(lvl2)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "ab0ac799-a832-40af-9cd4-b5b3934cfee4", + "metadata": {}, + "outputs": [], + "source": [ + "# Examine the First Image (Before DQ Flags Set)\n", + "ref_image = lvl2[0]\n", + "print(ref_image)\n", + "ref_fits = ImageModel(ref_image)\n", + "ref_data = ref_fits.data\n", + "\n", + "# The scale should highlight the background noise so it is possible to see all faint sources.\n", + "norm1 = simple_norm(ref_data,stretch='log',min_cut=4.5,max_cut=5)\n", + "\n", + "plt.figure(figsize=(20,12))\n", + "plt.imshow(ref_data,origin='lower',norm=norm1,cmap='gray')\n", + "clb = plt.colorbar()\n", + "clb.set_label('MJy/Str', labelpad=-40, y=1.05, rotation=0)\n", + "plt.gca().tick_params(axis='both',color='none')\n", + "plt.xlabel('Pixels')\n", + "plt.ylabel('Pixels')\n", + "plt.show()" + ] + }, + { + "cell_type": "raw", + "id": "b89fed70-2f65-4fa6-9242-87f98372e887", + "metadata": {}, + "source": [ + "# Examine the First Image (Before DQ Flags Set)\n", + "ref_image = lvl2[0]\n", + "print(ref_image)\n", + "ref_fits = fits.open(ref_image)\n", + "ref_data = fits.open(ref_image)['SCI',1].data\n", + "norm1 = simple_norm(ref_data,stretch='linear',min_cut=-1,max_cut=10)\n", + "\n", + "plt.imshow(ref_data, origin='lower',norm=norm1,cmap='gray')\n", + "plt.gca().tick_params(labelcolor='none',axis='both',color='none')\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "4ab947da-be0a-4cde-88eb-2d947adf2b81", + "metadata": {}, + "outputs": [], + "source": [ + "# Change all DQ flagged pixels to NANs\n", + "\n", + "# Reference for JWST DQ Flag Definitions: https://jwst-pipeline.readthedocs.io/en/latest/jwst/references_general/references_general.html\n", + "# In this case, we choose all DQ > 10, but users are encouraged to choose their own values accordingly.\n", + "for file in lvl2:\n", + " ref_fits = ImageModel(file)\n", + " data = ref_fits.data\n", + " dq = ref_fits.dq\n", + " data[dq >= 10]=np.nan\n", + " ref_fits.data=data\n", + " ref_fits.save(file)" + ] + }, + { + "cell_type": "raw", + "id": "21d3d16a-84a8-44af-84b2-c7eb0dbf230e", + "metadata": {}, + "source": [ + "# Change all DQ flagged pixels to NANs\n", + "\n", + "# Reference for JWST DQ Flag Definitions: https://jwst-pipeline.readthedocs.io/en/latest/jwst/references_general/references_general.html\n", + "# In this case, we choose all DQ > 10, but users are encouraged to choose their own values accordingly.\n", + "for file in lvl2:\n", + " hdul = fits.open(file, mode='update')\n", + " data = fits.open(file)['SCI',1].data\n", + " dq = fits.open(file)['DQ',1].data\n", + " data[dq >= 10]=np.nan\n", + " hdul['SCI',1].data=data\n", + " hdul.flush()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "87673dd8-1bc6-4663-bebc-8e2434974ceb", + "metadata": {}, + "outputs": [], + "source": [ + "# Examine the First Image (After DQ Flags Set)\n", + "ref_image = lvl2[0]\n", + "print(ref_image)\n", + "ref_fits = ImageModel(ref_image)\n", + "ref_data = ref_fits.data\n", + "\n", + "# The scale should highlight the background noise so it is possible to see all faint sources.\n", + "norm1 = simple_norm(ref_data,stretch='log',min_cut=4.5,max_cut=5)\n", + "\n", + "plt.figure(figsize=(20,12))\n", + "plt.imshow(ref_data,origin='lower',norm=norm1,cmap='gray')\n", + "clb = plt.colorbar()\n", + "clb.set_label('MJy/Str', labelpad=-40, y=1.05, rotation=0)\n", + "plt.gca().tick_params(axis='both',color='none')\n", + "plt.xlabel('Pixels')\n", + "plt.ylabel('Pixels')\n", + "plt.show()" + ] + }, + { + "cell_type": "raw", + "id": "30f99a12-b3c9-4cf7-a8cc-a98f848dc6d1", + "metadata": { + "tags": [] + }, + "source": [ + "# Examine the First Image (After DQ Flags Set)\n", + "ref_image = lvl2[0]\n", + "print(ref_image)\n", + "ref_fits = fits.open(ref_image)\n", + "ref_data = fits.open(ref_image)['SCI',1].data\n", + "norm1 = simple_norm(ref_data,stretch='linear',min_cut=-1,max_cut=10)\n", + "\n", + "plt.imshow(ref_data, origin='lower',norm=norm1,cmap='gray')\n", + "plt.gca().tick_params(labelcolor='none',axis='both',color='none')\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "4858f4d9-dbc9-40f1-a00c-80339cc69fae", + "metadata": {}, + "outputs": [], + "source": [ + "# Zoom in to see the source. In this case, our source is from MIRI Program ID #1028, a Calibration Program.\n", + "# We are using Visit 006, which targets the A5V dwarf 2MASSJ17430448+6655015\n", + "# Reference Link: http://simbad.cds.unistra.fr/simbad/sim-basic?Ident=2MASSJ17430448%2B6655015&submit=SIMBAD+search\n", + "\n", + "source_location = SkyCoord('17:43:04.4879','+66:55:01.837',unit=(u.hourangle,u.deg))\n", + "ref_wcs = ref_fits.get_fits_wcs()\n", + "ref_y,ref_x = skycoord_to_pixel(source_location,ref_wcs)\n", + "ref_cutout = extract_array(ref_data,(21,21),(ref_x,ref_y))\n", + "\n", + "# The scale should highlight the background noise so it is possible to see all faint sources.\n", + "norm1 = simple_norm(ref_cutout,stretch='log',min_cut=4.3,max_cut=15)\n", + "plt.imshow(ref_cutout,origin='lower',norm=norm1,cmap='gray')\n", + "clb = plt.colorbar()\n", + "clb.set_label('MJy/Str', labelpad=-40, y=1.05, rotation=0)\n", + "plt.title('PID1028,Obs006')\n", + "plt.xlabel('Pixels')\n", + "plt.ylabel('Pixels')\n", + "plt.gca().tick_params(axis='both',color='none')\n", + "plt.show()" + ] + }, + { + "cell_type": "raw", + "id": "31bcb979-870c-4a98-8d89-eb825d05e45b", + "metadata": { + "tags": [] + }, + "source": [ + "# Zoom in to see the source\n", + "source_location = SkyCoord('17:43:04.4879','+66:55:01.837',unit=(u.hourangle,u.deg))\n", + "ref_y,ref_x = skycoord_to_pixel(source_location,wcs.WCS(ref_fits['SCI',1],ref_fits))\n", + "ref_cutout = extract_array(ref_data,(11,11),(ref_x,ref_y))\n", + "norm1 = simple_norm(ref_cutout,stretch='linear',min_cut=-1,max_cut=10)\n", + "plt.imshow(ref_cutout, origin='lower',\n", + " norm=norm1,cmap='gray')\n", + "plt.title('PID1028,Obs006')\n", + "plt.gca().tick_params(labelcolor='none',axis='both',color='none')\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "id": "51fd4846-8cd6-4c60-8438-c6352488a3cc", + "metadata": {}, + "source": [ + "#### DEV NOTES:
\n", + "1. Add documentation to space_phot describing the get_jwst_psf wrapper.
\n", + "2. Add documentation on using webbpsf.calc_psf instead of space_phot.get_jwst_psf to have more control over your PSF.
" + ] + }, + { + "cell_type": "raw", + "id": "bd083573-528d-4f2a-a6b0-e6c762dcbb77", + "metadata": {}, + "source": [ + "# Get PSF from WebbPSF\n", + "\n", + "# The function get_jwst_psf is a space_phot wrapper for the WebbPSF calc_psf and CreatePSFLibrary functions. It is intended to provide easy, but limited, access to WebbPSF.\n", + "# The source code for psfs = space_phot.get_jwst_psf(jwst_obs,source_location) is as follows:\n", + "c = webbpsf.gridded_library.CreatePSFLibrary(inst,inst.filter, num_psfs = 1, psf_location = (x,y), fov_pixels = psf_width,\n", + " detectors=st_obs.detector,save=False)\n", + "grid = c.create_grid()\n", + "\n", + "# In most cases, we assume the user wishes to simply grab the default WebbPSF. Advanced users should consider using webbpsf.\n", + "# In that case, you would use the code below.
\n", + "\n", + "psf_model = inst.calc_psf(oversample=4,normalize='last') # creates a psf model in webbpsf\n", + "psfs = photutils.psf.FittableImageModel(psf[0].data*16,normalize=False,oversampling=oversampling) # turns the psf into a fittableimage model to improve processing time\n", + "\n", + "# There are more advanced methods for generating your WebbPSF, but those are beyond the scope of this notebook.\n", + "# Useful reference: https://webbpsf.readthedocs.io/en/latest/api/webbpsf.JWInstrument.html#webbpsf.JWInstrument.calc_psf\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "d67d57b9", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "# Get the PSF from WebbPSF using defaults.\n", + "\n", + "jwst_obs = space_phot.observation2(lvl2)\n", + "psfs = space_phot.get_jwst_psf(jwst_obs,source_location)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "cee4ba3b-2a71-4685-ae1f-5beb60c7d4ad", + "metadata": {}, + "outputs": [], + "source": [ + "# The scale should highlight the background noise so it is possible to see all faint sources.\n", + "\n", + "ref_cutout = extract_array(psfs[0].data,(41,41),(122,122))\n", + "#ref_cutout = psfs[0].data\n", + "norm1 = simple_norm(ref_cutout,stretch='log',min_cut=0.0,max_cut=0.2)\n", + "plt.imshow(ref_cutout,origin='lower',norm=norm1,cmap='gray')\n", + "clb = plt.colorbar()\n", + "clb.set_label('MJy/Str', labelpad=-40, y=1.05, rotation=0)\n", + "plt.title('WebbPSF Model')\n", + "plt.xlabel('Pixels')\n", + "plt.ylabel('Pixels')\n", + "plt.gca().tick_params(axis='both',color='none')\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "id": "2f35f730-ce0f-450a-8606-dcc0eb8643ad", + "metadata": {}, + "source": [ + "#### DEV NOTES:
\n", + "\n", + "1. Add improved documentation on how the fitting is performed, including packages, background, parameters, and how a calibrated flux is calculated, etc.\n", + "1. Add improved documentation on reading the plots and metrics from psf_photometry. And details about the Table Columns.\n", + "2. Add improved documentation on the various parameters and how to optimize them based on the metrics." + ] + }, + { + "cell_type": "markdown", + "id": "bf53473e-9671-4f80-9d59-29c4ed8bab09", + "metadata": {}, + "source": [ + "#### Notes on the PSF Fitting in Space_Phot:
\n", + "\n", + "https://st-phot.readthedocs.io/en/latest/examples/plot_a_psf.html#jwst-images\n", + "As noted above, improved documentation will be coming. For now, here are some important points to consider.\n", + "\n", + "All fitting is performed with Astropy's Photutils. As with any photometry program, the printed statistical errors are good indicators of your success.\n", + "\n", + "There are different fitting techniques, but when the fit_flux parameter is set to 'single', the source is fit simultaneously in all Level2 images. There is good reason for this outlined in a paper for PSF fitting in Hubble: https://iopscience.iop.org/article/10.1086/316630/pdf\n", + "\n", + "As a result, the flux and corresponding error take into account a single overall fit. As part of this, the fitting therefore assumes a constant zero point across all images. While this is not exactly true, it is typically true to within 1\\% and good enough for our purposes. Users can alternatively fit with the fit_flux parameter set to 'multi', which treats each image independently. The final flux must therefore be averaged.\n", + "\n", + "When you run space_phot, you will see some additional diagnositics displayed in the cell. At the top, the % value printed is the fraction of flux remaining in the residual and can be considered a good indicator of a successful model and subtraction. Next are three columns displaying the original, the model, and the residual, respectively, for each Level2 image. Finally, there are corner plots suggesting the success of the fits (more documentation and explanation of these plots is coming).\n", + "\n", + "In this case, you will notice a systematic trend in the residuals. The PSF is oversubtracted in the centermost pixel and undersubtracted in the wings. The cause is unknown. It may be due to a poor PSF model. The user should consider generating more complex PSF models, but that is beyond the scope of this notebook. Nonetheless, the residual value is pretty good so the overall statistical error is relatively small." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "c9a1b447", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "# Do PSF Photometry using space_phot (details of fitting are in documentation)\n", + "\n", + "jwst_obs.psf_photometry(psfs,source_location,bounds={'flux':[-10000,10000],\n", + " 'centroid':[-2,2],\n", + " 'bkg':[0,50]},\n", + " fit_width=7,\n", + " fit_centroid='pixel',\n", + " fit_bkg=True,\n", + " fit_flux='single')\n", + "jwst_obs.plot_psf_fit()\n", + "plt.show()\n", + "\n", + "jwst_obs.plot_psf_posterior(minweight=.0005)\n", + "plt.show()\n", + "\n", + "print(jwst_obs.psf_result.phot_cal_table)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "46b39817-d45e-4504-9e0a-281ef405f695", + "metadata": {}, + "outputs": [], + "source": [ + "jwst_obs.psf_result.phot_cal_table" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "31bd70f0", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "# Print Magnitude from Table\n", + "# As noted above, As a result, the flux and corresponding error take into account a single overall fit. \n", + "# Therefore, there is no need to average the resulting magnitudes or errors. They should all be the same to within their individual zero-point differences (typically <1%).\n", + "\n", + "mag_lvl2_arr = jwst_obs.psf_result.phot_cal_table['mag']\n", + "magerr_lvl2_arr = jwst_obs.psf_result.phot_cal_table['magerr']\n", + "\n", + "print(mag_lvl2_arr,magerr_lvl2_arr)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "af070d39-4c07-4801-b98f-dd9d3be0d7fe", + "metadata": {}, + "outputs": [], + "source": [ + "jwst_obs_fast = space_phot.observation2(lvl2[0])" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "8c4785b0-4caf-4fef-848f-4842e91a6eac", + "metadata": {}, + "outputs": [], + "source": [ + "ref_x" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "ed225d52-9d2c-4641-a37f-922e8fa696fc", + "metadata": {}, + "outputs": [], + "source": [ + "centers = [ref_x, ref_y]\n", + "jwst_obs_fast.fast_psf(psfs[0],centers)" + ] + }, + { + "cell_type": "markdown", + "id": "a57f9275", + "metadata": {}, + "source": [ + "### 3.2-Single, Level3 Mosaicked File ###" + ] + }, + { + "cell_type": "markdown", + "id": "62a15eb8-2a9d-4f5a-895a-b18dc33b24e4", + "metadata": {}, + "source": [ + "Despite the above discussion on performing PSF photometry on the pre-mosaiced data products, space_phot has the functionality to create a mosaiced Level3 PSF at a given single position on the detector based on the Level2 images. The advantage to this is the ability to perform PSF photometry on the deep, stacked data in cases where faint sources are expected to have prohibitively low signal-to-noise in Level2 data. The disadvantage is the amount of time required to make mosaiced Level3 PSF, so that this method is most useful when dealing with a small number of low signal-to-noise sources.
\n", + "\n", + "Useful references:
\n", + "Space-Phot Documentation on Level3 Fitting: https://space-phot.readthedocs.io/en/latest/examples/plot_a_psf.html#level-3-psf
" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "7bd91d15", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "# Level3 data file the same as above.\n", + "\n", + "lvl3" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "31eb83ad-b752-4a51-ac9a-6e44f59495f9", + "metadata": {}, + "outputs": [], + "source": [ + "# Now do the same photometry on the Level 3 Data\n", + "ref_image = lvl3[0]\n", + "ref_fits = ImageModel(ref_image)\n", + "ref_data = ref_fits.data\n", + "\n", + "# The scale should highlight the background noise so it is possible to see all faint sources.\n", + "norm1 = simple_norm(ref_data,stretch='log',min_cut=4.5,max_cut=5)\n", + "\n", + "plt.figure(figsize=(20,12))\n", + "plt.imshow(ref_data,origin='lower',norm=norm1,cmap='gray')\n", + "clb = plt.colorbar()\n", + "clb.set_label('MJy/Str', labelpad=-40, y=1.05, rotation=0)\n", + "plt.gca().tick_params(axis='both',color='none')\n", + "plt.xlabel('Pixels')\n", + "plt.ylabel('Pixels')\n", + "plt.show()" + ] + }, + { + "cell_type": "raw", + "id": "375143a8-e3a5-4a44-b8bd-ba91875d1aa4", + "metadata": { + "tags": [] + }, + "source": [ + "# Now do the same photometry on the Level 3 Data\n", + "ref_image = lvl3[0]\n", + "ref_fits = fits.open(ref_image)\n", + "ref_data = fits.open(ref_image)['SCI',1].data\n", + "norm1 = simple_norm(ref_data,stretch='linear',min_cut=-1,max_cut=10)\n", + "\n", + "plt.imshow(ref_data, origin='lower',\n", + " norm=norm1,cmap='gray')\n", + "plt.gca().tick_params(labelcolor='none',axis='both',color='none')\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "0cd06c0a", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "source_location = SkyCoord('17:43:04.4879','+66:55:01.837',unit=(u.hourangle,u.deg))\n", + "\n", + "ref_wcs = ref_fits.get_fits_wcs()\n", + "ref_y,ref_x = skycoord_to_pixel(source_location,ref_wcs)\n", + "ref_cutout = extract_array(ref_data,(21,21),(ref_x,ref_y))\n", + "\n", + "# The scale should highlight the background noise so it is possible to see all faint sources.\n", + "norm1 = simple_norm(ref_cutout,stretch='log',min_cut=4.5,max_cut=30)\n", + "plt.imshow(ref_cutout,origin='lower',norm=norm1,cmap='gray')\n", + "clb = plt.colorbar()\n", + "clb.set_label('MJy/Str', labelpad=-40, y=1.05, rotation=0)\n", + "plt.title('PID1028,Obs006')\n", + "plt.xlabel('Pixels')\n", + "plt.ylabel('Pixels')\n", + "plt.gca().tick_params(axis='both',color='none')\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "7bea076d", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "# Get PSF from WebbPSF\n", + "\n", + "# The function get_jwst_psf is a space_phot wrapper for the WebbPSF calc_psf function and uses a lot of the same keywords.\n", + "# There are more advanced methods for generating your WebbPSF, but those are beyond the scope of this notebook.\n", + "# The defaults used by get_jwst_psf in this notebook are:\n", + " # oversample=4\n", + " # normalize='last'\n", + " # Non-distorted PSF\n", + "# Useful reference: https://webbpsf.readthedocs.io/en/latest/api/webbpsf.JWInstrument.html#webbpsf.JWInstrument.calc_psf\n", + "\n", + "jwst3_obs = space_phot.observation3(lvl3[0])\n", + "psf3 = space_phot.get_jwst3_psf(jwst_obs,jwst3_obs,source_location)#,num_psfs=4)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "d126d410-7a21-45e9-b0ad-8b9a8308a174", + "metadata": {}, + "outputs": [], + "source": [ + "# The scale should highlight the background noise so it is possible to see all faint sources.\n", + "ref_cutout = extract_array(psf3.data,(161,161),(200,200))\n", + "norm1 = simple_norm(ref_cutout,stretch='log',min_cut=0.0,max_cut=0.01)\n", + "plt.imshow(ref_cutout,origin='lower',norm=norm1,cmap='gray')\n", + "clb = plt.colorbar()\n", + "clb.set_label('MJy/Str', labelpad=-40, y=1.05, rotation=0)\n", + "plt.title('WebbPSF Model (Mosaiced)')\n", + "plt.xlabel('Pixels')\n", + "plt.ylabel('Pixels')\n", + "plt.gca().tick_params(axis='both',color='none')\n", + "plt.show()" + ] + }, + { + "cell_type": "markdown", + "id": "990505bf-184e-456f-9a7e-14df8316414c", + "metadata": {}, + "source": [ + "#### Notes on the PSF Fitting in Space_Phot:
\n", + "\n", + "https://st-phot.readthedocs.io/en/latest/examples/plot_a_psf.html#jwst-images\n", + "As noted above, improved documentation will be coming. For now, here are some important points to consider.\n", + "See detailed notes in Section 3.1 above about the fitting process and diagnostics\n", + "\n", + "In addition, consider here that jwst3_obs is generating a Level3 PSF by using the JWST pipeline to resample and combine multiple Level2 PSFs. The Level2 PSFs are generated at the precise location of the source in each Level2 file to account for detector level effects. The resampling uses default resampling paramters. However, users should be aware that if they performed customized resampling for their Level2 data products, they should use similar resampling steps for their PSF below. \n", + "\n", + "### DEV NOTE: Create a resampling step with more resampling paramaters." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "27525a0a", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "# Do PSF Photometry using space_phot (details of fitting are in documentation)\n", + "# See detailed notes in Section 3.1 above about the fitting process and diagnostics\n", + "\n", + "jwst3_obs.psf_photometry(psf3,source_location,bounds={'flux':[-10000,10000],\n", + " 'centroid':[-2,2],\n", + " 'bkg':[0,50]},\n", + " fit_width=9,\n", + " fit_bkg=True,\n", + " fit_flux=True)\n", + "\n", + "jwst_obs.plot_psf_fit()\n", + "plt.show()\n", + "\n", + "jwst_obs.plot_psf_posterior(minweight=.0005)\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "dc1f930a", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "mag_lvl3psf = jwst3_obs.psf_result.phot_cal_table['mag'][0]\n", + "magerr_lvl3psf = jwst3_obs.psf_result.phot_cal_table['magerr'][0]\n", + "print(round(mag_lvl2_arr[0],4),round(magerr_lvl2_arr[0],4))\n", + "print(round(mag_lvl3psf,5),round(magerr_lvl3psf,5))" + ] + }, + { + "cell_type": "markdown", + "id": "4920c8b1-da54-434c-a8b9-4427c3157a62", + "metadata": {}, + "source": [ + "## Good agreement between Level2 and level3 results!" + ] + }, + { + "cell_type": "markdown", + "id": "5b5f0ad5-b59e-4eff-8687-f6a2199d8bd9", + "metadata": {}, + "source": [ + "4.-Faint/Upper Limit, Single Object\n", + "------------------" + ] + }, + { + "cell_type": "markdown", + "id": "1dc60da1-0f4d-4e5e-a109-f7352dfd0fdc", + "metadata": {}, + "source": [ + "The purpose of this section is to illustrate how to calculate an upper limit using PSF photometry a blank part of the sky. " + ] + }, + { + "cell_type": "markdown", + "id": "8af6f83a-5925-44d5-be0f-facd2316d1ca", + "metadata": {}, + "source": [ + "### 4.1-Multiple, Level2 Files ###" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "29ac6f64", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "### Level 3 Files\n", + "lvl3 = ['mast/01028/jw01028-o006_t001_miri_f770w_i2d.fits']\n", + "lvl3" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "3ebedb35", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "### Create Level 2 Data List from ASN files\n", + "prefix = \"./mast/01028/\"\n", + "asn = glob.glob(prefix+'jw01028-o006_*_image3_00004_asn.json')\n", + "with open(asn[0],\"r\") as fi:\n", + " lvl2 = []\n", + " for ln in fi:\n", + " #print(ln)\n", + " if ln.startswith(' \"expname\":'):\n", + " x = ln[2:].split(':')\n", + " y = x[1].split('\"')\n", + " lvl2.append(prefix+y[1])\n", + "print(lvl2)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "322b5a27", + "metadata": {}, + "outputs": [], + "source": [ + "# Change all DQ flagged pixels to NANs\n", + "\n", + "# Reference for JWST DQ Flag Definitions: https://jwst-pipeline.readthedocs.io/en/latest/jwst/references_general/references_general.html\n", + "# In this case, we choose all DQ > 10, but users are encouraged to choose their own values accordingly.\n", + "for file in lvl2:\n", + " ref_fits = ImageModel(file)\n", + " data = ref_fits.data\n", + " dq = ref_fits.dq\n", + " data[dq >= 10]=np.nan\n", + " ref_fits.data=data\n", + " ref_fits.save(file)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "2b5c7fea-3df3-4434-92cc-c7d8afa531dc", + "metadata": {}, + "outputs": [], + "source": [ + "# Examine the First Image (After DQ Flags Set)\n", + "ref_image = lvl2[0]\n", + "print(ref_image)\n", + "ref_fits = ImageModel(ref_image)\n", + "ref_data = ref_fits.data\n", + "\n", + "# The scale should highlight the background noise so it is possible to see all faint sources.\n", + "norm1 = simple_norm(ref_data,stretch='log',min_cut=4.5,max_cut=5)\n", + "\n", + "plt.figure(figsize=(20,12))\n", + "plt.imshow(ref_data,origin='lower',norm=norm1,cmap='gray')\n", + "clb = plt.colorbar()\n", + "clb.set_label('MJy/Str', labelpad=-40, y=1.05, rotation=0)\n", + "plt.gca().tick_params(axis='both',color='none')\n", + "plt.xlabel('Pixels')\n", + "plt.ylabel('Pixels')\n", + "plt.show()" + ] + }, + { + "cell_type": "raw", + "id": "eac6b4b7-0d02-45f8-9f96-8ea82c589126", + "metadata": { + "tags": [] + }, + "source": [ + "# Examine the First Image\n", + "ref_image = lvl2[0]\n", + "print(ref_image)\n", + "ref_fits = fits.open(ref_image)\n", + "ref_data = fits.open(ref_image)['SCI',1].data\n", + "norm1 = simple_norm(ref_data,stretch='linear',min_cut=-1,max_cut=10)\n", + "\n", + "plt.imshow(ref_data, origin='lower',norm=norm1,cmap='gray')\n", + "plt.gca().tick_params(labelcolor='none',axis='both',color='none')\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "6f4ed7dc-bb03-4dca-9908-2130d96f7c63", + "metadata": {}, + "outputs": [], + "source": [ + "# Pick a blank part of the sky to calculate the upper limit\n", + "\n", + "source_location = SkyCoord('17:43:00.0332','+66:54:42.677',unit=(u.hourangle,u.deg))\n", + "ref_wcs = ref_fits.get_fits_wcs()\n", + "ref_y,ref_x = skycoord_to_pixel(source_location,ref_wcs)\n", + "ref_cutout = extract_array(ref_data,(21,21),(ref_x,ref_y))\n", + "\n", + "# The scale should highlight the background noise so it is possible to see all faint sources.\n", + "norm1 = simple_norm(ref_cutout,stretch='log',min_cut=4.5,max_cut=5)\n", + "plt.imshow(ref_cutout,origin='lower',norm=norm1,cmap='gray')\n", + "clb = plt.colorbar()\n", + "clb.set_label('MJy/Str', labelpad=-40, y=1.05, rotation=0)\n", + "plt.title('PID1028,Obs006')\n", + "plt.xlabel('Pixels')\n", + "plt.ylabel('Pixels')\n", + "plt.gca().tick_params(axis='both',color='none')\n", + "plt.show()" + ] + }, + { + "cell_type": "raw", + "id": "7b2c277e-43d8-4429-9612-2b4440cbdbb9", + "metadata": { + "tags": [] + }, + "source": [ + "# Pick a blank part of the sky to calculate the upper limit\n", + "\n", + "source_location = SkyCoord('17:43:00.0332','+66:54:42.677',unit=(u.hourangle,u.deg))\n", + "ref_y,ref_x = skycoord_to_pixel(source_location,wcs.WCS(ref_fits['SCI',1],ref_fits))\n", + "ref_cutout = extract_array(ref_data,(11,11),(ref_x,ref_y))\n", + "norm1 = simple_norm(ref_cutout,stretch='linear',min_cut=-1,max_cut=10)\n", + "plt.imshow(ref_cutout, origin='lower',\n", + " norm=norm1,cmap='gray')\n", + "plt.title('PID1028,Obs006')\n", + "plt.gca().tick_params(labelcolor='none',axis='both',color='none')\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "72b8c907", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "# Get PSF from WebbPSF\n", + "# The function get_jwst_psf is a space_phot wrapper for the WebbPSF calc_psf function and uses a lot of the same keywords.\n", + "# There are more advanced methods for generating your WebbPSF, but those are beyond the scope of this notebook.\n", + "# The defaults used by get_jwst_psf in this notebook are:\n", + " # oversample=4\n", + " # normalize='last'\n", + " # Non-distorted PSF\n", + "# Useful reference: https://webbpsf.readthedocs.io/en/latest/api/webbpsf.JWInstrument.html#webbpsf.JWInstrument.calc_psf\n", + "\n", + "jwst_obs = space_phot.observation2(lvl2)\n", + "psfs = space_phot.get_jwst_psf(jwst_obs,source_location)#,num_psfs=4)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "04793a1d-5406-4516-9db8-5845a9f97194", + "metadata": {}, + "outputs": [], + "source": [ + "# The scale should highlight the background noise so it is possible to see all faint sources.\n", + "ref_cutout = extract_array(psf3.data,(161,161),(200,200))\n", + "norm1 = simple_norm(ref_cutout,stretch='log',min_cut=0.0,max_cut=0.01)\n", + "plt.imshow(ref_cutout,origin='lower',norm=norm1,cmap='gray')\n", + "clb = plt.colorbar()\n", + "clb.set_label('MJy/Str', labelpad=-40, y=1.05, rotation=0)\n", + "plt.title('WebbPSF Model')\n", + "plt.xlabel('Pixels')\n", + "plt.ylabel('Pixels')\n", + "plt.gca().tick_params(axis='both',color='none')\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "a2615569", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "# Do PSF Photometry using space_phot (details of fitting are in documentation)\n", + "# https://st-phot.readthedocs.io/en/latest/examples/plot_a_psf.html#jwst-images\n", + "\n", + "jwst_obs.psf_photometry(psfs,source_location,bounds={'flux':[-10,1000],\n", + " #'centroid':[-2,2],\n", + " 'bkg':[0,50]},\n", + " fit_width=5,\n", + " fit_bkg=True,\n", + " fit_centroid='fixed',\n", + " fit_flux='single')\n", + "jwst_obs.plot_psf_fit()\n", + "plt.show()\n", + "\n", + "jwst_obs.plot_psf_posterior(minweight=.0005)\n", + "plt.show()\n", + "\n", + "print(jwst_obs.psf_result.phot_cal_table)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "234c8a45", + "metadata": {}, + "outputs": [], + "source": [ + "# Print Upper Limits\n", + "# As noted above, As a result, the flux and corresponding error take into account a single overall fit. \n", + "# Therefore, there is no need to average the resulting magnitudes or errors. They should all be the same to within their individual zero-point differences (typically <1%).\n", + "\n", + "\n", + "magupper_lvl2psf = jwst_obs.upper_limit(nsigma=5)\n", + "magupper_lvl2psf" + ] + }, + { + "cell_type": "markdown", + "id": "e325db03-6e9b-4f04-be7b-5e80063dd9b8", + "metadata": { + "tags": [] + }, + "source": [ + "### 4.2-Single, Level3 Mosaicked File ###" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "2e7ce46d", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "# Level3 data file the same as above.\n", + "\n", + "lvl3" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "ffe34f40-6e67-4357-af85-908f92deb889", + "metadata": {}, + "outputs": [], + "source": [ + "# Now do the same photometry on the Level 3 Data\n", + "ref_image = lvl3[0]\n", + "ref_fits = ImageModel(ref_image)\n", + "ref_data = ref_fits.data\n", + "\n", + "# The scale should highlight the background noise so it is possible to see all faint sources.\n", + "norm1 = simple_norm(ref_data,stretch='log',min_cut=4.5,max_cut=5)\n", + "\n", + "plt.figure(figsize=(20,12))\n", + "plt.imshow(ref_data,origin='lower',norm=norm1,cmap='gray')\n", + "clb = plt.colorbar()\n", + "clb.set_label('MJy/Str', labelpad=-40, y=1.05, rotation=0)\n", + "plt.gca().tick_params(axis='both',color='none')\n", + "plt.xlabel('Pixels')\n", + "plt.ylabel('Pixels')\n", + "plt.show()" + ] + }, + { + "cell_type": "raw", + "id": "e06609c2-a895-4666-b7c5-d3d3697e1923", + "metadata": { + "tags": [] + }, + "source": [ + "# Now do the same photometry on the Level 3 Data\n", + "ref_image = lvl3[0]\n", + "ref_fits = fits.open(ref_image)\n", + "ref_data = fits.open(ref_image)['SCI',1].data\n", + "norm1 = simple_norm(ref_data,stretch='linear',min_cut=-1,max_cut=10)\n", + "\n", + "plt.imshow(ref_data, origin='lower',\n", + " norm=norm1,cmap='gray')\n", + "plt.gca().tick_params(labelcolor='none',axis='both',color='none')\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "26c66637-1ee7-44b2-92a0-c042ac91e10d", + "metadata": {}, + "outputs": [], + "source": [ + "# Pick a blank part of the sky to calculate the upper limit\n", + "\n", + "source_location = SkyCoord('17:43:00.0332','+66:54:42.677',unit=(u.hourangle,u.deg))\n", + "ref_wcs = ref_fits.get_fits_wcs()\n", + "ref_y,ref_x = skycoord_to_pixel(source_location,ref_wcs)\n", + "ref_cutout = extract_array(ref_data,(21,21),(ref_x,ref_y))\n", + "\n", + "# The scale should highlight the background noise so it is possible to see all faint sources.\n", + "norm1 = simple_norm(ref_cutout,stretch='log',min_cut=4.5,max_cut=5)\n", + "plt.imshow(ref_cutout,origin='lower',norm=norm1,cmap='gray')\n", + "clb = plt.colorbar()\n", + "clb.set_label('MJy/Str', labelpad=-40, y=1.05, rotation=0)\n", + "plt.title('PID1028,Obs006')\n", + "plt.xlabel('Pixels')\n", + "plt.ylabel('Pixels')\n", + "plt.gca().tick_params(axis='both',color='none')\n", + "plt.show()" + ] + }, + { + "cell_type": "raw", + "id": "b652b4cf-0d19-49b7-8b09-ecc9c24c6b23", + "metadata": { + "tags": [] + }, + "source": [ + "# Pick a blank part of the sky to calculate the upper limit\n", + "\n", + "source_location = SkyCoord('17:43:00.0332','+66:54:42.677',unit=(u.hourangle,u.deg))\n", + "\n", + "ref_y,ref_x = skycoord_to_pixel(source_location,wcs.WCS(ref_fits['SCI',1],ref_fits))\n", + "ref_cutout = extract_array(ref_data,(11,11),(ref_x,ref_y))\n", + "norm1 = simple_norm(ref_cutout,stretch='linear',min_cut=-1,max_cut=10)\n", + "plt.imshow(ref_cutout, origin='lower',\n", + " norm=norm1,cmap='gray')\n", + "plt.title('PID1028,Obs006 (level 3)')\n", + "plt.gca().tick_params(labelcolor='none',axis='both',color='none')\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "50fbb856", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "# Get PSF from WebbPSF\n", + "# The function get_jwst_psf is a space_phot wrapper for the WebbPSF calc_psf function and uses a lot of the same keywords.\n", + "# There are more advanced methods for generating your WebbPSF, but those are beyond the scope of this notebook.\n", + "# The defaults used by get_jwst_psf in this notebook are:\n", + " # oversample=4\n", + " # normalize='last'\n", + " # Non-distorted PSF\n", + "# Useful reference: https://webbpsf.readthedocs.io/en/latest/api/webbpsf.JWInstrument.html#webbpsf.JWInstrument.calc_psf\n", + "\n", + "\n", + "jwst3_obs = space_phot.observation3(lvl3[0])\n", + "psf3 = space_phot.get_jwst3_psf(jwst_obs,jwst3_obs,source_location)#,num_psfs=4)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "37eaaa85-da28-4dd6-a1c5-727b718f1831", + "metadata": {}, + "outputs": [], + "source": [ + "# The scale should highlight the background noise so it is possible to see all faint sources.\n", + "ref_cutout = extract_array(psf3.data,(161,161),(200,200))\n", + "norm1 = simple_norm(ref_cutout,stretch='log',min_cut=0.0,max_cut=0.01)\n", + "plt.imshow(ref_cutout,origin='lower',norm=norm1,cmap='gray')\n", + "clb = plt.colorbar()\n", + "clb.set_label('MJy/Str', labelpad=-40, y=1.05, rotation=0)\n", + "plt.title('WebbPSF Model')\n", + "plt.xlabel('Pixels')\n", + "plt.ylabel('Pixels')\n", + "plt.gca().tick_params(axis='both',color='none')\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "fe699262", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "jwst3_obs.psf_photometry(psf3,source_location,bounds={'flux':[-1000,1000],\n", + " #'centroid':[-2,2],\n", + " 'bkg':[0,50]},\n", + " fit_width=9,\n", + " fit_bkg=True, \n", + " fit_centroid=False,\n", + " fit_flux=True)\n", + "\n", + "jwst3_obs.plot_psf_fit()\n", + "plt.show()\n", + "\n", + "jwst3_obs.plot_psf_posterior(minweight=.0005)\n", + "plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "ecaec8db", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "magupper_lvl3psf = jwst3_obs.upper_limit(nsigma=5)\n", + "print(round(magupper_lvl2psf[0],4))\n", + "print(round(magupper_lvl3psf[0],5))" + ] + }, + { + "cell_type": "markdown", + "id": "4e4996d2-6274-473d-b13c-f3848c27ad78", + "metadata": {}, + "source": [ + "## Note you can go significantly deeper with the Level3 combined data product" + ] + }, + { + "cell_type": "markdown", + "id": "9a969717-bbef-40b9-ac9b-f83dec99dc09", + "metadata": {}, + "source": [ + "5.-Stellar Field (LMC)\n", + "------------------" + ] + }, + { + "cell_type": "markdown", + "id": "da877310-fd47-41d6-afea-7fa725a546af", + "metadata": {}, + "source": [ + "#### In this case, we are going to do the same steps as in Section 3, but for multiple stars. The purpose is to illustrate the workflow and runtime for using space_phot on a large number of stars. We suggest that space_phot may be less optimal for large numbers of bright stars. Other programs, such as DOLPHOT or Photutils, may be better suited for this use case. The primary advantage to space_phot is on faint, single sources. But it can be extended to a larger number if desired." + ] + }, + { + "cell_type": "markdown", + "id": "32bdafe6-db19-4080-9587-b9785c2f7fa7", + "metadata": {}, + "source": [ + "### 5.1-Multiple, Level2 Files ###" + ] + }, + { + "cell_type": "markdown", + "id": "b618756f", + "metadata": {}, + "source": [ + "##### Now do the same thing for a larger group of stars and test for speed" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "838bd76d", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "### Level 3 Files\n", + "#lvl3 = ['mastDownload/JWST/jw01171-o004_t001_miri_f560w/jw01171-o004_t001_miri_f560w_i2d.fits']\n", + "lvl3 = [\"./mast/01171/jw01171-o004_t001_miri_f560w_i2d.fits\"]\n", + "lvl3" + ] + }, + { + "cell_type": "raw", + "id": "76e791f7", + "metadata": {}, + "source": [ + "asnfile = 'jw02079-o004_20230622t175524_spec2_00001_asn.json'\n", + "asn_data = json.load(open(asnfile))\n", + "print(asn_data['products'][0]['members'][0]['expname'])" + ] + }, + { + "cell_type": "raw", + "id": "f39e7500-0909-4823-acaf-50a04f5fda30", + "metadata": {}, + "source": [ + "lvl2 = glob.glob('mastDownload/JWST/jw01171004*/*cal.fits')\n", + "lvl2" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "73aba802", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "### Level 2 Files\n", + "#lvl2 = glob.glob('mastDownload/JWST/jw01171004*/*cal.fits')\n", + "\n", + "lvl2 = glob.glob('./mast/01171/jw01171004*cal.fits')\n", + "lvl2" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "57f9d790", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "# Find Stars in Level 3 File\n", + "\n", + "# Get rough estimate of background (There are Better Ways to Do Background Subtraction)\n", + "bkgrms = MADStdBackgroundRMS()\n", + "mmm_bkg = MMMBackground()\n", + "\n", + "#im = fits.open(lvl3[0]) \n", + "#w = wcs.WCS(im['SCI',1])\n", + "\n", + "ref_fits = ImageModel(lvl3[0])\n", + "w = ref_fits.get_fits_wcs()\n", + "\n", + "#std = bkgrms(im[1].data)\n", + "#bkg = mmm_bkg(im[1].data)\n", + "#data_bkgsub = im[1].data.copy()\n", + "\n", + "std = bkgrms(ref_fits.data)\n", + "bkg = mmm_bkg(ref_fits.data)\n", + "data_bkgsub = ref_fits.data.copy()\n", + "data_bkgsub -= bkg \n", + "#sigma_psf = 1.636 #pixls for F770W\n", + "fwhm_psf = 1.882 #pixels for F560W\n", + "threshold = 5.\n", + "\n", + "#daofind = DAOStarFinder(threshold=threshold * std, fwhm=fwhm_psf, exclude_border=True)\n", + "daofind = DAOStarFinder(threshold=threshold * std, fwhm=fwhm_psf, exclude_border=True, min_separation=10)\n", + "\n", + "found_stars = daofind(data_bkgsub)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "e4cee97c", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "found_stars.pprint_all(max_lines=10)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "cb42faab-f55f-41fb-9145-65adc8cb28b5", + "metadata": {}, + "outputs": [], + "source": [ + "# plot the found stars\n", + "norm = simple_norm(data_bkgsub, 'sqrt', percent=99)\n", + "fig, ax = plt.subplots(figsize=(10, 10))\n", + "ax.imshow(data_bkgsub, origin='lower', norm=norm)\n", + "\n", + "xypos = zip(found_stars['xcentroid'], found_stars['ycentroid'])\n", + "aper = CircularAperture(xypos, r=10)\n", + "aper.plot(ax, color='red');" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "7c7d793b", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "# Filter out only stars you want\n", + "\n", + "plt.figure(figsize=(12, 8))\n", + "plt.clf()\n", + "\n", + "ax1 = plt.subplot(2, 1, 1)\n", + "\n", + "ax1.set_xlabel('mag')\n", + "ax1.set_ylabel('sharpness')\n", + "\n", + "xlim0 = np.min(found_stars['mag']) - 0.25\n", + "xlim1 = np.max(found_stars['mag']) + 0.25\n", + "ylim0 = np.min(found_stars['sharpness']) - 0.15\n", + "ylim1 = np.max(found_stars['sharpness']) + 0.15\n", + "\n", + "ax1.set_xlim(xlim0, xlim1)\n", + "ax1.set_ylim(ylim0, ylim1)\n", + "\n", + "#ax1.xaxis.set_major_locator(ticker.AutoLocator())\n", + "#ax1.xaxis.set_minor_locator(ticker.AutoMinorLocator())\n", + "#ax1.yaxis.set_major_locator(ticker.AutoLocator())\n", + "#ax1.yaxis.set_minor_locator(ticker.AutoMinorLocator())\n", + "\n", + "ax1.scatter(found_stars['mag'], found_stars['sharpness'], s=10, color='k')\n", + "\n", + "sh_inf = 0.40\n", + "sh_sup = 0.82\n", + "#mag_lim = -5.0\n", + "lmag_lim = -3.0\n", + "umag_lim = -5.0\n", + "\n", + "ax1.plot([xlim0, xlim1], [sh_sup, sh_sup], color='r', lw=3, ls='--')\n", + "ax1.plot([xlim0, xlim1], [sh_inf, sh_inf], color='r', lw=3, ls='--')\n", + "ax1.plot([lmag_lim, lmag_lim], [ylim0, ylim1], color='r', lw=3, ls='--')\n", + "ax1.plot([umag_lim, umag_lim], [ylim0, ylim1], color='r', lw=3, ls='--')\n", + "\n", + "ax2 = plt.subplot(2, 1, 2)\n", + "\n", + "ax2.set_xlabel('mag')\n", + "ax2.set_ylabel('roundness')\n", + "\n", + "ylim0 = np.min(found_stars['roundness2']) - 0.25\n", + "ylim1 = np.max(found_stars['roundness2']) - 0.25\n", + "\n", + "ax2.set_xlim(xlim0, xlim1)\n", + "ax2.set_ylim(ylim0, ylim1)\n", + "\n", + "#ax2.xaxis.set_major_locator(ticker.AutoLocator())\n", + "#ax2.xaxis.set_minor_locator(ticker.AutoMinorLocator())\n", + "#ax2.yaxis.set_major_locator(ticker.AutoLocator())\n", + "#ax2.yaxis.set_minor_locator(ticker.AutoMinorLocator())\n", + "\n", + "round_inf = -0.40\n", + "round_sup = 0.40\n", + "\n", + "ax2.scatter(found_stars['mag'], found_stars['roundness2'], s=10, color='k')\n", + "\n", + "ax2.plot([xlim0, xlim1], [round_sup, round_sup], color='r', lw=3, ls='--')\n", + "ax2.plot([xlim0, xlim1], [round_inf, round_inf], color='r', lw=3, ls='--')\n", + "ax2.plot([lmag_lim, lmag_lim], [ylim0, ylim1], color='r', lw=3, ls='--')\n", + "ax2.plot([umag_lim, umag_lim], [ylim0, ylim1], color='r', lw=3, ls='--')\n", + "\n", + "plt.tight_layout()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "6ac852af", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "mask = ((found_stars['mag'] < lmag_lim) & (found_stars['mag'] > umag_lim) & (found_stars['roundness2'] > round_inf)\n", + " & (found_stars['roundness2'] < round_sup) & (found_stars['sharpness'] > sh_inf) \n", + " & (found_stars['sharpness'] < sh_sup) & (found_stars['xcentroid'] > 100) & (found_stars['xcentroid'] < 700)\n", + " & (found_stars['ycentroid'] > 100) & (found_stars['ycentroid'] < 700))\n", + "\n", + "found_stars_sel = found_stars[mask]\n", + "\n", + "print('Number of stars found originally:', len(found_stars))\n", + "print('Number of stars in final selection:', len(found_stars_sel))\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "567f81f5", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "found_stars_sel" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "f7fd293e-a8f4-4c34-ab99-23d71adee080", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "# Convert pixel to wcs coords\n", + "from astropy.wcs.utils import skycoord_to_pixel\n", + "skycoords = w.pixel_to_world(found_stars_sel['xcentroid'], found_stars_sel['ycentroid'])\n", + "#skycoords = skycoords[9:] #Cutting down the list artificially to make things run faster\n", + "len(skycoords)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "e6c46e19", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "# Change all DQ flagged pixels to NANs\n", + "for file in lvl2:\n", + " #hdul = fits.open(file, mode='update')\n", + " #data = fits.open(file)['SCI',1].data\n", + " #dq = fits.open(file)['DQ',1].data\n", + " #data[dq == 262657]=np.nan\n", + " #data[dq == 262661]=np.nan\n", + " #hdul['SCI',1].data=data\n", + " #hdul.flush()\n", + "\n", + " ref_fits = ImageModel(file)\n", + " data = ref_fits.data\n", + " dq = ref_fits.dq\n", + " data[dq >= 10]=np.nan\n", + " ref_fits.data=data\n", + " ref_fits.save(file)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "5516a64f", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "# Create a grid for fast lookup using WebbPSF. The larger the number of grid points, the better the photometric precision.\n", + "# Developer note. Would be great to have a fast/approximate look up table. \n", + "jwst_obs = space_phot.observation2(lvl2)\n", + "grid = space_phot.util.get_jwst_psf_grid(jwst_obs,num_psfs=16)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "38760103-e702-4b6b-bd5e-7ed12c67a6d8", + "metadata": {}, + "outputs": [], + "source": [ + "from astropy.io import ascii\n", + "from astropy.table import QTable\n", + "\n", + "t = QTable([skycoords], names=[\"skycoord\"])\n", + "t.write('skycoord.ecsv',overwrite=True)" + ] + }, + { + "cell_type": "markdown", + "id": "abe8c20c-087a-4b59-9f0c-26ccdb3dcbec", + "metadata": {}, + "source": [ + "#### DEV NOTE: allow space_phot to take pixel coordinates as input" + ] + }, + { + "cell_type": "markdown", + "id": "0199d630-d0e9-415b-a664-5725b45fa2c0", + "metadata": {}, + "source": [ + "#### DEV NOTE: Open WebbPSF ticket for webbpsf.gridded_library.CreatePSFLibrary, which is called in get_jwst_psf_from_grid below. The reason for this is that the output from that command currently prints: Position 2/4 centroid: (201.3017584451758, 201.31980818321827), but this is the lower corner of the array. We want a grid across the entire array. The notebook can currently be run as is." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "ebc298b0-633b-4915-8f03-b510db3d81cc", + "metadata": {}, + "outputs": [], + "source": [ + "lvl2" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "562f9d50-0884-47c7-923b-9d4e019c247b", + "metadata": {}, + "outputs": [], + "source": [ + "!pip freeze | grep space_phot" + ] + }, + { + "cell_type": "markdown", + "id": "1953b4f4-4ec3-4322-abbd-d53a1954c1f5", + "metadata": {}, + "source": [ + "#### Dev Note: Need Grid Capability for Level3 Data for fater PSF Generation" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "b85e222f", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "# Now Loop Through All Stars and Build Photometry Table\n", + "# Readers should refer to all diagnostics discussed above. \n", + "# It should be noted that empty plots correspond to LVL2 files with dither positions that do not cover that particular coordinate.\n", + "\n", + "import pandas as pd\n", + "import warnings\n", + "warnings.simplefilter('ignore')\n", + "counter = 0.\n", + "badindex = []\n", + "\n", + "jwst_obs = space_phot.observation2(lvl2)\n", + "localbkg_estimator = LocalBackground(5, 10, bkg_estimator=MMMBackground())\n", + "\n", + "for source_location in skycoords:\n", + " #print(source_location)\n", + " tic = time.perf_counter()\n", + " print('Starting',counter+1., ' of', len(skycoords), ':', source_location)\n", + " #psfs = space_phot.get_jwst_psf(jwst_obs,source_location)#,num_psfs=4)\n", + " psfs = space_phot.util.get_jwst_psf_from_grid(jwst_obs,source_location,grid)\n", + " xys = [jwst_obs.wcs_list[i].world_to_pixel(source_location) for i in range(jwst_obs.n_exposures)]\n", + " bkg = [localbkg_estimator(jwst_obs.data_arr_pam[i],xys[i][0],xys[i][0]) for i in range(jwst_obs.n_exposures)]\n", + " print(bkg)\n", + " jwst_obs.psf_photometry(psfs,source_location,bounds={'flux':[-100000,100000],\n", + " 'centroid':[-2.,2.]},\n", + " #'bkg':[0,50]},\n", + " fit_width=5,\n", + " fit_bkg=False,\n", + " fit_centroid='wcs',\n", + " background=bkg,\n", + " fit_flux='single',\n", + " maxiter=None)\n", + " \n", + " jwst_obs.plot_psf_fit()\n", + " plt.show()\n", + " ra = jwst_obs.psf_result.phot_cal_table['ra'][0]\n", + " dec = jwst_obs.psf_result.phot_cal_table['dec'][0]\n", + "\n", + " fit_location = SkyCoord(ra,dec,unit=u.deg)\n", + "\n", + " \n", + " \n", + " jwst_obs.psf_photometry(psfs,fit_location,bounds={'flux':[-100000,100000],\n", + " 'centroid':[-2.,2.]},\n", + " #'bkg':[0,50]},\n", + " fit_width=5,\n", + " fit_bkg=False,\n", + " background=bkg,\n", + " fit_centroid='fixed',\n", + " fit_flux='single',\n", + " maxiter=None)\n", + "#\n", + "# jwst_obs.plot_psf_posterior(minweight=.0005)\n", + "# plt.show()\n", + " jwst_obs.plot_psf_fit()\n", + " plt.show()\n", + " ra = jwst_obs.psf_result.phot_cal_table['ra'][0]\n", + " dec = jwst_obs.psf_result.phot_cal_table['dec'][0]\n", + " mag_arr = jwst_obs.psf_result.phot_cal_table['mag']\n", + " magerr_arr = jwst_obs.psf_result.phot_cal_table['magerr']\n", + " mag_lvl2psf = np.mean(mag_arr)\n", + " magerr_lvl2psf = math.sqrt(sum(p**2 for p in magerr_arr))\n", + "\n", + " #print('Finished',counter, ' of', len(skycoords), ':', source_location)\n", + " if counter == 0:\n", + " df = pd.DataFrame(np.array([[ra, dec, mag_lvl2psf, magerr_lvl2psf]]), columns=['ra','dec','mag','magerr'])\n", + " else:\n", + " df = pd.concat([df, pd.DataFrame(np.array([[ra, dec, mag_lvl2psf, magerr_lvl2psf]]), columns=['ra','dec','mag','magerr'])], ignore_index=True)\n", + " counter = counter + 1.\n", + " \n", + " toc = time.perf_counter()\n", + " print(\"Elapsed Time for Photometry:\", toc - tic)\n", + " #sys.exit()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "13518464-e91b-4fa0-9e64-9cc84b0a56ef", + "metadata": {}, + "outputs": [], + "source": [ + "xys" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "d4fbc1e5-96c5-48bb-b463-5af8970f1a2f", + "metadata": {}, + "outputs": [], + "source": [ + "df" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "dcd3bc1b-0263-4fb6-8217-108f4a88be4e", + "metadata": {}, + "outputs": [], + "source": [ + "# Write to File\n", + "\n", + "df.to_csv('miri_photometry_space_phot_lvl2.txt', index=False) " + ] + }, + { + "cell_type": "markdown", + "id": "3604e260-2da4-4f43-b306-fb7cd65e738b", + "metadata": {}, + "source": [ + "### 5.2-Single, Level3 Mosaicked File ###" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "a57e893d-92cb-4de6-8c64-69911b691246", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "80.705186\t-69.437532" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "24dbbba6-6d1a-40b2-9028-de916cdc76e4", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "# Now do the same photometry on the Level 3 Data\n", + "ref_image = lvl3[0]\n", + "\n", + "#ref_fits = fits.open(ref_image)\n", + "#ref_data = fits.open(ref_image)['SCI',1].data\n", + "\n", + "ref_fits = ImageModel(ref_image)\n", + "ref_data = ref_fits.data\n", + "norm1 = simple_norm(ref_data,stretch='linear',min_cut=0.5,max_cut=5)\n", + "\n", + "plt.figure(figsize=(20,12))\n", + "plt.imshow(ref_data,origin='lower',norm=norm1,cmap='gray')\n", + "clb = plt.colorbar()\n", + "clb.set_label('MJy/Str', labelpad=-40, y=1.05, rotation=0)\n", + "plt.gca().tick_params(axis='both',color='none')\n", + "plt.xlabel('Pixels')\n", + "plt.ylabel('Pixels')\n", + "plt.show()\n", + "\n", + "#plt.imshow(ref_data, origin='lower',norm=norm1,cmap='gray')\n", + "#plt.gca().tick_params(labelcolor='none',axis='both',color='none')\n", + "#plt.show()" + ] + }, + { + "cell_type": "markdown", + "id": "b6489fb5-e99c-4204-aff6-49402842f5c0", + "metadata": {}, + "source": [ + "#### Dev Note: Need Grid Capability for Level3 Data for fater PSF Generation" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "5be212d8-c43a-478e-98ed-b1877e44a347", + "metadata": { + "tags": [] + }, + "outputs": [], + "source": [ + "# Get PSF from WebbPSF and drizzle it to the source location\n", + "jwst3_obs = space_phot.observation3(lvl3[0])\n", + "#psf3 = space_phot.get_jwst3_psf(jwst_obs,source_location,num_psfs=4)\n", + "#grid = space_phot.util.get_jwst_psf_grid(jwst_obs,num_psfs=4)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "e4c85327-4b59-4228-8610-4a85909fb397", + "metadata": {}, + "outputs": [], + "source": [ + "lvl3[0]" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "717c2fc5-1e72-48c7-98aa-4b8216aac567", + "metadata": {}, + "outputs": [], + "source": [ + "skycoords" + ] + }, + { + "cell_type": "markdown", + "id": "1881aeea-d5c7-4e1c-9669-19351f9c1157", + "metadata": {}, + "source": [ + "#### Readers should refer to all diagnostics discussed above. In general, this loop shows the difficulty in doing PSF photometry on a wide variety of stars (brightness, distribution on the detector, etc) without visual inspection. Especially when dealing with low SNR sources.This is true for all photometry packages. Users should inspect the associated metrics and consider optimizing the parameters for specific stars of interest. Nonetheless, the success of the fits can always be quantified in with the associated error bars.\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "6f414c09-c62a-4645-9a29-b1042288661d", + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "922accc4-2179-4e03-ad60-2beeb594faea", + "metadata": { + "scrolled": true, + "tags": [] + }, + "outputs": [], + "source": [ + "# Now Loop Through All Stars and Build Photometry Table\n", + "\n", + "import pandas as pd\n", + "counter = 0.\n", + "badindex = []\n", + "\n", + "for source_location in skycoords[1:3]:\n", + " #print(source_location)\n", + " tic = time.perf_counter()\n", + " print('Starting',counter+1., ' of', len(skycoords), ':', source_location)\n", + " #psf3 = space_phot.util.get_jwst_psf_from_grid(jwst3_obs,source_location,grid)\n", + " psf3 = space_phot.get_jwst3_psf(jwst_obs,jwst3_obs,source_location,num_psfs=4)\n", + " xys = [jwst3_obs.wcs.world_to_pixel(source_location)]\n", + " bkg = [localbkg_estimator(jwst3_obs.data,xys[0][0],xys[0][1])]\n", + " print(bkg)\n", + " jwst3_obs.psf_photometry(psf3,source_location,bounds={'flux':[-100000,100000],\n", + " 'centroid':[-2.,2.]},\n", + " #'bkg':[0,50]},\n", + " fit_width=5,\n", + " fit_bkg=False,\n", + " background=bkg,\n", + " fit_flux=True)\n", + "\n", + " ra = jwst3_obs.psf_result.phot_cal_table['ra'][0]\n", + " dec = jwst3_obs.psf_result.phot_cal_table['dec'][0]\n", + " fit_location = SkyCoord(ra,dec,unit=u.deg)\n", + " jwst3_obs.aperture_photometry(fit_location,encircled_energy=70)\n", + " print(jwst3_obs.aperture_result.phot_cal_table['mag'])\n", + " \n", + " jwst3_obs.psf_photometry(psf3,fit_location,bounds={'flux':[-100000,100000],\n", + " 'centroid':[-2.,2.]},\n", + " #'bkg':[0,50]},\n", + " fit_width=5,\n", + " fit_bkg=False,\n", + " fit_centroid=False,\n", + " background=bkg,\n", + " fit_flux=True)\n", + " # fit_flux='single',\n", + " # maxiter=10000)\n", + " #jwst3_obs.psf_photometry(psf3,source_location,bounds={'flux':[-10000,10000],\n", + " # 'centroid':[-2,2],\n", + " # 'bkg':[0,50]},\n", + " # fit_width=9,\n", + " # fit_bkg=True,\n", + " # fit_flux=True)\n", + "\n", + " jwst3_obs.plot_psf_fit()\n", + " plt.show()\n", + "\n", + " ra = jwst3_obs.psf_result.phot_cal_table['ra'][0]\n", + " dec = jwst3_obs.psf_result.phot_cal_table['dec'][0]\n", + " mag_lvl3psf = jwst3_obs.psf_result.phot_cal_table['mag'][0]\n", + " magerr_lvl3psf = jwst3_obs.psf_result.phot_cal_table['magerr'][0]\n", + "\n", + " #print('Finished',counter, ' of', len(skycoords), ':', source_location)\n", + " if counter == 0:\n", + " df = pd.DataFrame(np.array([[ra, dec, mag_lvl3psf, magerr_lvl3psf]]), columns=['ra','dec','mag','magerr'])\n", + " else:\n", + " df = pd.concat([df, pd.DataFrame(np.array([[ra, dec, mag_lvl3psf, magerr_lvl3psf]]), columns=['ra','dec','mag','magerr'])], ignore_index=True)\n", + " counter = counter + 1.\n", + " toc = time.perf_counter()\n", + " print(\"Elapsed Time for Photometry:\", toc - tic)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "db250a73-d9a6-4f38-a42b-335f66d73bde", + "metadata": {}, + "outputs": [], + "source": [ + "lvl2" + ] + }, + { + "cell_type": "markdown", + "id": "4cb0557a-dee2-4be3-9513-83bb9671d71e", + "metadata": {}, + "source": [ + "**Important Note**: When not to use. Due to the sensitivity of the space_phot parameters, this tool is not meant to be used for a large sample of stars (i.e., Section 5 below). If a user would like to use space_phot on more than one source, they should carefully construct a table of parameters that are carefully refined for each source." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "5af6df8c-85ea-4aac-b6b7-53ac630fa3e0", + "metadata": {}, + "outputs": [], + "source": [ + "df" + ] + }, + { + "cell_type": "markdown", + "id": "5630029f-31d1-42cd-8454-225e86cabc48", + "metadata": {}, + "source": [ + "
" + ] + }, + { + "cell_type": "markdown", + "id": "843b5201-6f57-46f0-9da0-b738714178d3", + "metadata": {}, + "source": [ + "\"Space" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "798fd5bf-3269-4f5f-8cc0-9d3cb82cbcaf", + "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" + }, + "toc-showcode": false + }, + "nbformat": 4, + "nbformat_minor": 5 +}