<jupyter_start><jupyter_text>Creating a simulated starfield image and performing photometry using `webbpsf` and `photutils.GriddedPSFModel`  Imports* *matplotlib inline* for creating non-interactive plots* *numpy* for numerical computations* *astropy.io.fits* for accessing FITS files* *astropy.stats.gaussian_sigma_to_fwhm* for converting $\sigma$ to FWHM* *astropy.table.Table* for constructing Astropy table* *astropy.utils.misc.NumpyRNGContext* for reproducible random result* *astropy.visualization* for normalized image display* *photutils* for photometry related classes and functions* *webbpsf* for JWST PSF related functions<jupyter_code>%matplotlib inline

import numpy as np
import matplotlib.pyplot as plt
from astropy.io import fits
from astropy.stats import gaussian_sigma_to_fwhm
from astropy.table import Table
from astropy.utils.misc import NumpyRNGContext
from astropy.visualization import imshow_norm, PercentileInterval, SqrtStretch
from photutils import BasicPSFPhotometry
from photutils import DBSCANGroup, MMMBackground
from photutils.datasets import make_noise_image
from webbpsf.gridded_library import display_psf_grid
from webbpsf.utils import to_griddedpsfmodel<jupyter_output><empty_output><jupyter_text>IntroductionThe simulated image will be made using a spatially-dependent JWST PSF generated by `webbpsf`.  `photutils` will be then used to perform PSF-fitting photometry on the simulated JWST NIRCam image.The intended audience are astronomers already familiar with PSF photometry, the FITS file format, and basic Python usage. In-depth scientific and technical backgrounds will not be discussed nor taught in this notebook.**Background on JWST PSF simulation:** https://jwst.stsci.edu/science-planning/proposal-planning-toolbox/psf-simulation-tool-webbpsf**Background on PSF photometry:** https://photutils.readthedocs.io/en/stable/psf.html**Some notes about the file:** It is pre-generated grid of JWST PSFs for NIRCam A1 detector in F090W filter. The file is hosted on STScI Box file sharing system.**Defining some terms:*** **F090W:** A NIRCam filter* **FITS:** Flexible Image Transport System (https://fits.gsfc.nasa.gov/fits_primer.html)* **HDU:** Header/Data Unit in a FITS file* **JWST:** James Webb Space Telescope* **NIRCam:** Near-Infrared Camera on JWST (https://jwst.stsci.edu/instrumentation/nircam)* **NRCA1:** NIRCam A1 detector* **PSF:** Point-spread function  Loading dataWe start by loading into FITS HDU a pre-generated grid of JWST PSFs for NIRCam A1 detector in F090W filter.<jupyter_code>url = 'https://stsci.box.com/shared/static/6h8wsr2ts0t24s79cxmnyk8bldt0vv3i.fits'
hdu = fits.open(url)<jupyter_output><empty_output><jupyter_text>File information<jupyter_code>hdu.info()<jupyter_output><empty_output><jupyter_text>Astropy has saved a cached copy of this file under your local cache directory for Astropy; default is `~/.astropy/cache/download`.This file only has a single extension that holds a data cube. The cube consists of 16 PSF images, each with dimension of 404 by 404 pixels. Each PSF is oversampled by a factor of 4, so the actual detector pixel dimension covered by the PSF is 101 by 101.Now, let's look at the first PSF in the cube.<jupyter_code>data = hdu[0].data[0]
plt.figure(figsize=(9, 9))
imshow_norm(data, interval=PercentileInterval(99.), stretch=SqrtStretch())<jupyter_output><empty_output><jupyter_text>Creating a PSF modelWe now create a `photutils.GriddedPSFModel` from the FITS HDU.<jupyter_code>psfmodel = to_griddedpsfmodel(hdu, ext=0)
psfmodel<jupyter_output><empty_output><jupyter_text>`psfmodel` is an instance of a `photutils.GriddedPSFModel`.  Let's view some of it's attributes. The oversampled PSF grid is stored internally as a 3D `numpy.ndarray`.  The first axis represents the number of PSFs, while the second and third axes represent the `(ny, nx)` shape of the 2D PSFs.Here there are 16 PSFs (from a 4x4 reference grid) and the shape of each is 404x404 pixels.<jupyter_code>psfmodel.data.shape<jupyter_output><empty_output><jupyter_text>The `oversampling` attribute is the PSF oversampling factor.<jupyter_code>psfmodel.oversampling<jupyter_output><empty_output><jupyter_text>The `grid_xypos` attribute contains a list of the `(x, y)` positions of reference PSFs.  The PSF at an arbitrary `(x, y)` position is interpolated from the grid of reference PSFs.<jupyter_code>psfmodel.grid_xypos<jupyter_output><empty_output><jupyter_text>This model contains a 4x4 grid of reference PSFs.<jupyter_code>len(psfmodel.grid_xypos)<jupyter_output><empty_output><jupyter_text>The `meta` attribute is dictionary holding detailed information how the PSF model grid was generated by `webbpsf`. For example, `psfmodel.meta['fovpixel']` would indicate that "field of view in pixels (full array)" is 101 pixels.<jupyter_code>psfmodel.meta
fovpixel = psfmodel.meta['fovpixel']
print(fovpixel[1], 'is', fovpixel[0])<jupyter_output><empty_output><jupyter_text>We can use the `webbpsf.gridded_library.display_psf_grid` function to visualize the PSF grid.<jupyter_code>display_psf_grid(psfmodel)<jupyter_output><empty_output><jupyter_text>Creating simulated starfieldNow let's use the `psfmodel` grid to create a simulated starfield image. First, we'll define 25 stars with random positions and fluxes. We are using a small image with limited number of stars for faster computation; feel free to use a larger image and/or more stars when you run this notebook locally, if desired.<jupyter_code>shape = (512, 512)
data = np.zeros(shape)
nstars = 25
flux_max = 200000.

with NumpyRNGContext(12345):  # Seed for repeatability
    xx = np.random.uniform(low=0, high=shape[1], size=nstars)
    yy = np.random.uniform(low=0, high=shape[0], size=nstars)
    zz = np.random.uniform(low=0, high=flux_max, size=nstars)<jupyter_output><empty_output><jupyter_text>Now we'll evaluate the model at these positions and fluxes.<jupyter_code>eval_xshape = np.int(np.ceil(psfmodel.data.shape[2] / psfmodel.oversampling))
eval_yshape = np.int(np.ceil(psfmodel.data.shape[1] / psfmodel.oversampling))

for xxi, yyi, zzi in zip(xx, yy, zz):
    x0 = np.int(np.floor(xxi - (eval_xshape - 1) / 2.))
    y0 = np.int(np.floor(yyi - (eval_yshape - 1) / 2.))
    x1 = x0 + eval_xshape
    y1 = y0 + eval_yshape
    
    if x0 < 0:
        x0 = 0
    if y0 < 0:
        y0 = 0
    if x1 > shape[1]:
        x1 = shape[1]
    if y1 > shape[0]:
        y1 = shape[0]
        
    y, x = np.mgrid[y0:y1, x0:x1]
    data[y, x] += psfmodel.evaluate(x=x, y=y, flux=zzi, x_0=xxi, y_0=yyi)<jupyter_output><empty_output><jupyter_text>Let's add some noise to the image.<jupyter_code>noise = make_noise_image(data.shape, 'gaussian', mean=0, stddev=2, random_state=123)
data += noise<jupyter_output><empty_output><jupyter_text>Let's display the simulated JWST NIRCam A1 detector image in F090W filter.<jupyter_code>plt.figure(figsize=(9, 9))
imshow_norm(data, interval=PercentileInterval(99.), stretch=SqrtStretch())<jupyter_output><empty_output><jupyter_text>Performing PSF-fitting photometry with `GriddPSFModel`We can also use `GriddedPSFModel` to perform PSF-fitting photometry. For this simple example, we'll use `photutils.BasicPSFPhotometry` to perform the photometry.  Here we're inputting a table of the initial positions (instead of inputting a star finder object). First, we'll create an Astropy table defining the initial guesses of star positions.<jupyter_code>init_tbl = Table()
init_tbl['x_0'] = xx.astype(int)
init_tbl['y_0'] = yy.astype(int)
init_tbl['flux_0'] = zz.astype(int)<jupyter_output><empty_output><jupyter_text>Then, we define the parameters to initialize a `BasicPSFPhotometry` instance.<jupyter_code>sigma_psf = 3.
daogroup = DBSCANGroup(2.0 * sigma_psf * gaussian_sigma_to_fwhm)
mmm_bkg = MMMBackground()
fit_shape = (eval_yshape, eval_xshape)
phot = BasicPSFPhotometry(daogroup, mmm_bkg, psfmodel, fit_shape, finder=None, 
                          aperture_radius=3.)<jupyter_output><empty_output><jupyter_text>Now, we call the instance on the data. This step might take a few seconds. See https://github.com/astropy/photutils/issues/631 for more information.<jupyter_code>tbl = phot(data, init_guesses=init_tbl)<jupyter_output><empty_output><jupyter_text>The result is an Astropy table containing the initial and fitted values of the star positions and fluxes. We can format the displayed values by setting `table[column_name].format`.<jupyter_code>tbl['x_fit'].format = '%.1f'
tbl['y_fit'].format = '%.1f'
tbl['flux_fit'].format = '%.4e'
tbl['flux_unc'].format = '%.4e'
tbl['x_0_unc'].format = '%.4e'
tbl['y_0_unc'].format = '%.4e'
tbl<jupyter_output><empty_output><jupyter_text>We can also view the residual image of the best-fit model PSF image subtracted from the data.<jupyter_code>diff = phot.get_residual_image()
plt.figure(figsize=(9, 9))
imshow_norm(diff, interval=PercentileInterval(99.), stretch=SqrtStretch())
plt.colorbar()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.10087326466636812, 'fraction_numerical': 0.017913121361397222, 'mean_word_length': 6.028324154209284, 'pattern_counts': {'":': 0, '<': 83, '<?xml version=': 0, '>': 83, 'https://': 6, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '20', 'dir': 'jupyter-structured-clean-dedup', 'id': '532479', 'license': 'BSD-3-Clause', 'n_tokens_mistral': 3203, 'n_tokens_neox': 2983, 'n_words': 866, 'path': 'notebooks/JWST/NIRCam/GriddedPSFModel_simulation_and_photometry/GriddedPSFModel_simulation_and_photometry.ipynb', 'repo_name': 'cmbennett01/notebooks', 'size': '8932'}