diff --git a/notebooks/NIRCam/psf_photometry/NIRCam_PSF_Photometry_Example.ipynb b/notebooks/NIRCam/psf_photometry/NIRCam_PSF_Photometry_Example.ipynb
index 5333b3926..5c834bd6e 100644
--- a/notebooks/NIRCam/psf_photometry/NIRCam_PSF_Photometry_Example.ipynb
+++ b/notebooks/NIRCam/psf_photometry/NIRCam_PSF_Photometry_Example.ipynb
@@ -8,7 +8,7 @@
"\n",
"\n",
"**Use case:** PSF photometry, creating a PSF, derive Color-Magnitude Diagram.
\n",
- "**Data:** NIRCam simulated images obtained using [MIRAGE](https://jwst-docs.stsci.edu/jwst-other-tools/mirage-data-simulator) and run through the [JWST pipeline](https://jwst-pipeline.readthedocs.io/en/latest/) of the Large Magellanic Cloud (LMC) Astrometric Calibration Field. Simulations is obtained using a 4-pt subpixel dither for three couples of wide filters: F070W, F115W, and F200W for the SW channel, and F277W, F356W, and F444W for the LW channel. We simulated only 1 NIRCam SW detector (i.e., \"NRCB1\"). For this example, we use Level-2 images (.cal, calibrated but not rectified) for two SW filters (i.e., F115W and F200W) and derive the photometry in each one of them. The images for the other filters are also available and can be used to test the notebook and/or different filters combination.
\n",
+ "**Data:** NIRCam simulated images of the Large Magellanic Cloud (LMC) Astrometric Calibration Field obtained using [MIRAGE](https://jwst-docs.stsci.edu/jwst-other-tools/mirage-jwst-data-simulator) and run through the [JWST pipeline](https://jwst-pipeline.readthedocs.io/en/latest/). The simulations are obtained using a 4-point subpixel dither for three pairs of wide filters: F070W, F115W, and F200W for the SW channel, and F277W, F356W, and F444W for the LW channel. We simulated only 1 NIRCam SW detector (i.e., \"NRCB1\"). For this example, we use Level-2 images (.cal, calibrated but not rectified) for two SW filters (F115W and F200W) and derive the photometry in each of them. The images for the other filters are also available and can be used to test the notebook and/or different filter combinations.
\n",
"**Tools:** photutils.
\n",
"**Cross-intrument:** NIRSpec, NIRISS, 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).
\n",
@@ -45,7 +45,7 @@
"cell_type": "markdown",
"metadata": {},
"source": [
- "**Note on pysynphot**: Data files for pysynphot are distributed separately by Calibration Reference Data System. They are expected to follow a certain directory structure under the root directory, identified by the PYSYN_CDBS environment variable that must be set prior to using this package. In the example below, the root directory is arbitrarily named /my/local/dir/trds/. \\\n",
+ "**Note on pysynphot**: Data files for pysynphot are distributed separately by the Calibration Reference Data System (CRDS). They are expected to follow a certain directory structure under the root directory, identified by the PYSYN_CDBS environment variable that must be set prior to using this package. In the example below, the root directory is arbitrarily named /my/local/dir/trds/. \\\n",
"export PYSYN_CDBS=/my/local/dir/trds/ \\\n",
"See documentation [here](https://pysynphot.readthedocs.io/en/latest/#installation-and-setup) for the configuration and download of the data files."
]
@@ -69,6 +69,7 @@
"import os\n",
"import tarfile\n",
"import time\n",
+ "import warnings\n",
"from urllib import request\n",
"\n",
"import numpy as np\n",
@@ -81,14 +82,15 @@
"from astropy.nddata import NDData\n",
"from astropy.stats import sigma_clipped_stats\n",
"from astropy.table import QTable, Table\n",
+ "from astropy.utils.exceptions import AstropyUserWarning\n",
"from astropy.visualization import simple_norm\n",
"from jwst.datamodels import ImageModel\n",
"from photutils.aperture import (CircularAnnulus, CircularAperture,\n",
" aperture_photometry)\n",
"from photutils.background import MADStdBackgroundRMS, MMMBackground\n",
"from photutils.detection import DAOStarFinder\n",
- "from photutils.psf import (EPSFBuilder, IterativePSFPhotometry, SourceGrouper,\n",
- " extract_stars)"
+ "from photutils.psf import (EPSFBuilder, GriddedPSFModel, IterativePSFPhotometry,\n",
+ " SourceGrouper, extract_stars)"
]
},
{
@@ -118,8 +120,8 @@
"\n",
"plt.rcParams['image.cmap'] = 'viridis'\n",
"plt.rcParams['image.origin'] = 'lower'\n",
- "plt.rcParams['axes.titlesize'] = plt.rcParams['axes.labelsize'] = 30\n",
- "plt.rcParams['xtick.labelsize'] = plt.rcParams['ytick.labelsize'] = 30\n",
+ "plt.rcParams['axes.titlesize'] = plt.rcParams['axes.labelsize'] = 14\n",
+ "plt.rcParams['xtick.labelsize'] = plt.rcParams['ytick.labelsize'] = 14\n",
"\n",
"font1 = {'family': 'helvetica', 'color': 'black', 'weight': 'normal', 'size': '12'}\n",
"font2 = {'family': 'helvetica', 'color': 'black', 'weight': 'normal', 'size': '20'}"
@@ -156,14 +158,14 @@
" os.makedirs(webbpsf_folder)\n",
" request.urlretrieve(boxlink, boxfile)\n",
" gzf = tarfile.open(boxfile)\n",
- " gzf.extractall(webbpsf_folder)\n",
+ " gzf.extractall(webbpsf_folder, filter='data')\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('./')"
+ " gzf.extractall('./', filter='data')"
]
},
{
@@ -213,7 +215,7 @@
" request.urlretrieve(boxlink_images_lev2, boxfile_images_lev2)\n",
"\n",
" tar = tarfile.open(boxfile_images_lev2, 'r')\n",
- " tar.extractall()\n",
+ " tar.extractall(filter='data')\n",
"\n",
" images_dir = './'\n",
" images = sorted(glob.glob(os.path.join(images_dir, \"*cal.fits\")))\n",
@@ -339,21 +341,17 @@
"metadata": {},
"outputs": [],
"source": [
- "plt.figure(figsize=(14, 14))\n",
+ "fig, ax = plt.subplots(ncols=2, figsize=(14, 14))\n",
"\n",
"for det in dets_short:\n",
" for i, filt in enumerate(filts_short):\n",
" image = fits.open(dict_images[det][filt]['images'][0])\n",
" data_sb = image[1].data\n",
- "\n",
- " ax = plt.subplot(1, len(filts_short), i + 1)\n",
- "\n",
- " plt.xlabel(\"X [px]\", fontdict=font2)\n",
- " plt.ylabel(\"Y [px]\", fontdict=font2)\n",
- " plt.title(filt, fontdict=font2)\n",
" norm = simple_norm(data_sb, 'sqrt', percent=99.)\n",
- "\n",
- " ax.imshow(data_sb, norm=norm, cmap='Greys')\n",
+ " ax[i].imshow(data_sb, norm=norm, cmap='Greys')\n",
+ " ax[i].set_xlabel(\"X [px]\", fontdict=font2)\n",
+ " ax[i].set_ylabel(\"Y [px]\", fontdict=font2)\n",
+ " ax[i].set_title(filt, fontdict=font2)\n",
"\n",
"plt.tight_layout()"
]
@@ -371,7 +369,7 @@
"source": [
"## I. Create the PSF model using WebbPSF\n",
"\n",
- "We create a dictionary that contains the PSF created using WebbPSF for the detectors and filters selected above."
+ "We create a dictionary that will contain the PSFs created using WebbPSF for the detectors and filters selected above."
]
},
{
@@ -395,7 +393,7 @@
"cell_type": "markdown",
"metadata": {},
"source": [
- "The function below allows to create a single PSF or a grid of PSFs and allows to save the PSF as a fits file. The model PSF are stored by default in the psf dictionary. For the grid of PSFs, users can select the number of PSFs to be created. The PSF can be created detector sampled or oversampled (the oversample can be changed inside the function).\n",
+ "The function below creates a single PSF or a grid of PSFs and allows one to save the PSF as a fits file. The model PSFs are stored by default in the psf dictionary. For the grid of PSFs, users can select the number of PSFs to be created. The PSF can be created detector sampled or oversampled (the oversampling factor can be changed inside the function).\n",
"\n",
"**Note**: The default source spectrum is, if `pysynphot` is installed, a G2V star spectrum from Castelli & Kurucz (2004). Without `pysynphot`, the default is a simple flat spectrum such that the same number of photons are detected at each wavelength."
]
@@ -406,7 +404,7 @@
"metadata": {},
"outputs": [],
"source": [
- "def create_psf_model(fov=11, create_grid=False, num=9, save_psf=False, detsampled=False):\n",
+ "def create_psf_model(det='NRCB1', fov=11, create_grid=False, num=9, save_psf=False, detsampled=False):\n",
" nrc = webbpsf.NIRCam()\n",
" nrc.detector = det \n",
" nrc.filter = filt\n",
@@ -425,31 +423,33 @@
" print(\"\")\n",
" print(f\"Creating a grid of PSF for filter {filt} and detector {det}\")\n",
" print(\"\")\n",
- " num = num\n",
"\n",
- " if save_psf:\n",
- " outname = f\"./PSF_{filt}_samp4_G5V_fov{fov}_npsfs{num}.fits\"\n",
- " nrc.psf_grid(num_psfs=num, oversample=4, source=src, all_detectors=False, fov_pixels=fov,\n",
- " save=True, outfile=outname, use_detsampled_psf=detsampled)\n",
+ " outname = f'nircam_{det}_{filt}_fovp{fov}_samp4_npsf{num}.fits'\n",
+ " if os.path.exists(outname):\n",
+ " grid_psf = GriddedPSFModel.read(outname)\n",
" else:\n",
" grid_psf = nrc.psf_grid(num_psfs=num, oversample=4, source=src, all_detectors=False,\n",
- " fov_pixels=fov, use_detsampled_psf=detsampled)\n",
- " dict_psfs_webbpsf[det][filt]['psf model grid'] = grid_psf\n",
+ " fov_pixels=fov, use_detsampled_psf=detsampled,\n",
+ " save=save_psf)\n",
+ "\n",
+ " dict_psfs_webbpsf[det][filt]['psf model grid'] = grid_psf\n",
+ " \n",
" else:\n",
" print(\"\")\n",
" print(f\"Creating a single PSF for filter {filt} and detector {det}\")\n",
" print(\"\")\n",
- " num = 1\n",
- " if save_psf:\n",
- " outname = f\"./PSF_{filt}_samp4_G5V_fov{fov}_npsfs{num}.fits\"\n",
- " nrc.psf_grid(num_psfs=num, oversample=4, source=src, all_detectors=False, fov_pixels=fov,\n",
- " save=True, outfile=outname, use_detsampled_psf=detsampled)\n",
+ "\n",
+ " outname = f'nircam_{det}_{filt}_fovp{fov}_samp4_npsf{num}.fits'\n",
+ " if os.path.exists(outname):\n",
+ " single_psf = GriddedPSFModel.read(outname)\n",
" else:\n",
- " single_psf = nrc.psf_grid(num_psfs=num, oversample=4, source=src, all_detectors=False,\n",
- " fov_pixels=fov, use_detsampled_psf=detsampled)\n",
- " dict_psfs_webbpsf[det][filt]['psf model single'] = single_psf\n",
+ " single_psf = nrc.psf_grid(num_psfs=1, oversample=4, source=src, all_detectors=False,\n",
+ " fov_pixels=fov, use_detsampled_psf=detsampled,\n",
+ " save=save_psf)\n",
+ "\n",
+ " dict_psfs_webbpsf[det][filt]['psf model single'] = single_psf\n",
"\n",
- " return "
+ " return dict_psfs_webbpsf "
]
},
{
@@ -467,7 +467,7 @@
"source": [
"for det in dets_short:\n",
" for filt in filts_short:\n",
- " create_psf_model(fov=11, num=25, create_grid=False, save_psf=False, detsampled=False)"
+ " create_psf_model(fov=11, num=1, create_grid=False, save_psf=True, detsampled=False)"
]
},
{
@@ -483,17 +483,17 @@
"metadata": {},
"outputs": [],
"source": [
- "plt.figure(figsize=(14, 14))\n",
+ "fig, ax = plt.subplots(ncols=2, figsize=(14, 14))\n",
"\n",
"for det in dets_short:\n",
" for i, filt in enumerate(filts_short):\n",
- " ax = plt.subplot(1, 2, i + 1)\n",
+ " img = dict_psfs_webbpsf[det][filt]['psf model single'].data[0]\n",
+ " norm_epsf = simple_norm(img, 'log', percent=99.)\n",
+ " ax[i].imshow(img, norm=norm_epsf)\n",
+ " ax[i].set_xlabel('X [px]', fontdict=font2)\n",
+ " ax[i].set_ylabel('Y [px]', fontdict=font2)\n",
+ " ax[i].set_title(filt, fontdict=font2)\n",
"\n",
- " norm_epsf = simple_norm(dict_psfs_webbpsf[det][filt]['psf model single'].data[0], 'log', percent=99.)\n",
- " ax.set_title(filt, fontsize=40)\n",
- " ax.imshow(dict_psfs_webbpsf[det][filt]['psf model single'].data[0], norm=norm_epsf)\n",
- " ax.set_xlabel('X [px]', fontsize=30)\n",
- " ax.set_ylabel('Y [px]', fontsize=30)\n",
"plt.tight_layout()"
]
},
@@ -512,7 +512,7 @@
"source": [
"for det in dets_short:\n",
" for filt in filts_short:\n",
- " create_psf_model(fov=11, num=25, create_grid=True, save_psf=False, detsampled=False)"
+ " create_psf_model(fov=11, num=25, create_grid=True, save_psf=True, detsampled=False)"
]
},
{
@@ -530,8 +530,17 @@
"metadata": {},
"outputs": [],
"source": [
- "webbpsf.gridded_library.display_psf_grid(dict_psfs_webbpsf[dets_short[0]][filts_short[0]]['psf model grid'],\n",
- " zoom_in=False, figsize=(14, 14))"
+ "griddedpsfmodel = dict_psfs_webbpsf[dets_short[0]][filts_short[0]]['psf model grid']\n",
+ "griddedpsfmodel.plot_grid(figsize=(10, 10));"
+ ]
+ },
+ {
+ "cell_type": "code",
+ "execution_count": null,
+ "metadata": {},
+ "outputs": [],
+ "source": [
+ "griddedpsfmodel.plot_grid(figsize=(10, 10), deltas=True, cmap='viridis', vmax_scale=0.3);"
]
},
{
@@ -540,7 +549,7 @@
"source": [
"## II. Create the PSF model building an Effective PSF (ePSF)\n",
"\n",
- "More information on the photutils Effective PSF can be found [here](https://photutils.readthedocs.io/en/stable/epsf.html).\n",
+ "More information on the photutils Effective PSF can be found [here](https://photutils.readthedocs.io/en/stable/user_guide/epsf.html).\n",
"\n",
"* Select the stars from the images we want to use for building the PSF. We use the [DAOStarFinder](https://photutils.readthedocs.io/en/stable/api/photutils.detection.DAOStarFinder.html) function to find bright stars in the images (setting a high detection threshold). [DAOStarFinder](https://photutils.readthedocs.io/en/stable/api/photutils.detection.DAOStarFinder.html#photutils.detection.DAOStarFinder) detects stars in an image using the DAOFIND ([Stetson 1987](https://ui.adsabs.harvard.edu/abs/1987PASP...99..191S/abstract)) algorithm. DAOFIND searches images for local density maxima that have a peak amplitude greater than `threshold` (approximately; threshold is applied to a convolved image) and have a size and shape similar to the defined 2D Gaussian kernel. \\\n",
" **Note**: The threshold and the maximum distance to the closest neighbour depend on the user science case (i.e.; number of stars in the field of view, crowding, number of bright sources, minimum number of stars required to build the ePSF, etc.) and must be modified accordingly. \n",
@@ -572,7 +581,6 @@
" dict_psfs_epsf[det][filt]['epsf grid'] = {}\n",
"\n",
" for i in np.arange(0, len(dict_images[det][filt]['images']), 1):\n",
- "\n",
" dict_psfs_epsf[det][filt]['table psf stars'][i + 1] = None\n",
" dict_psfs_epsf[det][filt]['epsf single'][i + 1] = None\n",
" dict_psfs_epsf[det][filt]['epsf grid'][i + 1] = None"
@@ -582,7 +590,7 @@
"cell_type": "markdown",
"metadata": {},
"source": [
- "Note that the unit of the Level-2 and Level-3 Images from the pipeline is MJy/sr (hence a surface brightness). The actual unit of the image can be checked from the header keyword **BUNIT**. The scalar conversion constant is copied to the header keyword **PHOTMJSR**, which gives the conversion from DN/s to megaJy/steradian. For our analysis we revert back to DN/s."
+ "Note that the unit of the Level-2 and Level-3 Images from the pipeline is MJy/sr (hence a surface brightness). The actual unit of the image can be checked from the header keyword **BUNIT**. The scalar conversion constant is copied to the header keyword **PHOTMJSR**, which gives the conversion from DN/s to MJy/steradian. For our analysis we revert back to DN/s."
]
},
{
@@ -752,11 +760,14 @@
"for det in dets_short:\n",
" for j, filt in enumerate(filts_short):\n",
" for i in np.arange(0, len(dict_images[det][filt]['images']), 1):\n",
- " build_epsf(det=det, filt=filt)\n",
+ " with warnings.catch_warnings():\n",
+ " # ignore warnings about stars close to image edge\n",
+ " warnings.simplefilter(\"ignore\", category=AstropyUserWarning) \n",
+ " build_epsf(det=det, filt=filt)\n",
"\n",
"toc = time.perf_counter()\n",
"\n",
- "print(\"Time to build the Effective PSF:\", toc - tic) "
+ "print(\"Time to build the Effective PSF:\", toc - tic)"
]
},
{
@@ -774,15 +785,14 @@
"metadata": {},
"outputs": [],
"source": [
- "plt.figure(figsize=(14, 14))\n",
+ "fig, ax = plt.subplots(ncols=2, figsize=(14, 14))\n",
"\n",
"for det in dets_short:\n",
" for i, filt in enumerate(filts_short):\n",
- " ax = plt.subplot(1, 2, i + 1)\n",
- "\n",
- " norm_epsf = simple_norm(dict_psfs_epsf[det][filt]['epsf single'][i + 1].data, 'log', percent=99.)\n",
- " plt.title(filt, fontsize=30)\n",
- " ax.imshow(dict_psfs_epsf[det][filt]['epsf single'][i + 1].data, norm=norm_epsf)"
+ " img = dict_psfs_epsf[det][filt]['epsf single'][i + 1].data\n",
+ " norm_epsf = simple_norm(img, 'log', percent=99.)\n",
+ " ax[i].imshow(img, norm=norm_epsf)\n",
+ " ax[i].set_title(filt, fontdict=font2)"
]
},
{
@@ -1003,6 +1013,7 @@
" \n",
" phot = IterativePSFPhotometry(psf_model, psf_shape, daofind,\n",
" grouper=grouper, fitter=fitter,\n",
+ " fitter_maxiters=500,\n",
" maxiters=2, aperture_radius=ap_radius[j])\n",
" result = phot(data_sub)\n",
" \n",
@@ -1012,7 +1023,7 @@
" print(f\"Time needed to perform photometry on image {i + 1}: {dtime:.2f} sec\")\n",
" print(f\"Number of sources detected in image {i + 1} for filter {filt}: {len(result)}\")\n",
" \n",
- " residual_image = phot.make_residual_image(data_sub, psf_shape)\n",
+ " residual_image = phot.make_residual_image(data_sub, psf_shape=psf_shape)\n",
" \n",
" dict_phot[det][filt]['residual images'][i + 1] = residual_image\n",
" dict_phot[det][filt]['output photometry tables'][i + 1] = result\n",
@@ -1094,35 +1105,25 @@
"metadata": {},
"outputs": [],
"source": [
- "plt.figure(figsize=(14, 14))\n",
+ "fig, ax = plt.subplots(ncols=2, nrows=2, figsize=(14, 14))\n",
"\n",
"for det in dets_short:\n",
" for i, filt in enumerate(filts_short):\n",
- "\n",
" image = fits.open(dict_images[det][filt]['images'][0])\n",
" data_sb = image[1].data\n",
- "\n",
- " ax = plt.subplot(2, len(filts_short), i + 1)\n",
- "\n",
- " plt.xlabel(\"X [px]\", fontdict=font2)\n",
- " plt.ylabel(\"Y [px]\", fontdict=font2)\n",
- " plt.title(filt, fontdict=font2)\n",
" norm = simple_norm(data_sb, 'sqrt', percent=99.)\n",
- "\n",
- " ax.imshow(data_sb, norm=norm, cmap='Greys')\n",
+ " ax[0, i].imshow(data_sb, norm=norm, cmap='Greys')\n",
+ " ax[0, i].set_xlabel(\"X [px]\", fontdict=font2)\n",
+ " ax[0, i].set_ylabel(\"Y [px]\", fontdict=font2)\n",
+ " ax[0, i].set_title(filt, fontdict=font2)\n",
"\n",
"for det in dets_short:\n",
" for i, filt in enumerate(filts_short):\n",
- "\n",
" res = dict_phot[det][filt]['residual images'][1]\n",
- "\n",
- " ax = plt.subplot(2, len(filts_short), i + 3)\n",
- "\n",
- " plt.xlabel(\"X [px]\", fontdict=font2)\n",
- " plt.ylabel(\"Y [px]\", fontdict=font2)\n",
- " norm = simple_norm(data_sb, 'sqrt', percent=99.)\n",
- "\n",
- " ax.imshow(res, norm=norm, cmap='Greys')\n",
+ " norm = simple_norm(res, 'sqrt', percent=99.)\n",
+ " ax[1, i].imshow(res, norm=norm, cmap='Greys')\n",
+ " ax[1, i].set_xlabel(\"X [px]\", fontdict=font2)\n",
+ " ax[1, i].set_ylabel(\"Y [px]\", fontdict=font2)\n",
"\n",
"plt.tight_layout()"
]
@@ -1163,14 +1164,14 @@
" request.urlretrieve(boxlink_cat_f115w, boxfile_cat_f115w)\n",
"\n",
" tar = tarfile.open(boxfile_cat_f115w, 'r')\n",
- " tar.extractall()\n",
+ " tar.extractall(filter='data')\n",
"\n",
" boxlink_cat_f200w = 'https://data.science.stsci.edu/redirect/JWST/jwst-data_analysis_tools/stellar_photometry/phot_cat_F200W.tar.gz'\n",
" boxfile_cat_f200w = './phot_cat_F200W.tar.gz'\n",
" request.urlretrieve(boxlink_cat_f200w, boxfile_cat_f200w)\n",
"\n",
" tar = tarfile.open(boxfile_cat_f200w, 'r')\n",
- " tar.extractall()\n",
+ " tar.extractall(filter='data')\n",
"\n",
" cat_dir = './'\n",
" phots_pkl_f115w = sorted(glob.glob(os.path.join(cat_dir, '*F115W*gridPSF*.pkl')))\n",
@@ -1816,7 +1817,7 @@
" request.urlretrieve(boxlink_images_lev3, boxfile_images_lev3)\n",
"\n",
" tar = tarfile.open(boxfile_images_lev3, 'r')\n",
- " tar.extractall()\n",
+ " tar.extractall(filter='data')\n",
"\n",
" images_dir = './'\n",
" files_singles = sorted(glob.glob(os.path.join(images_dir, \"*combined*fits\")))\n",
@@ -2650,7 +2651,7 @@
"## About this Notebook\n",
"\n",
"**Authors**: Matteo Correnti, JWST/NIRCam STScI Scientist II & Larry Bradley, Branch Deputy, Data Analysis Tools Branch\\\n",
- "**Updated on**: 2024-02-26"
+ "**Updated on**: 2024-10-22"
]
},
{
@@ -2679,7 +2680,7 @@
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
- "version": "3.11.8"
+ "version": "3.12.7"
}
},
"nbformat": 4,