%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
ValueErrorTraceback (most recent call last) <ipython-input-1-c9861cb722ec> in <module> 8 from astropy.utils.misc import NumpyRNGContext 9 from astropy.visualization import imshow_norm, PercentileInterval, SqrtStretch ---> 10 from photutils import BasicPSFPhotometry 11 from photutils import DBSCANGroup, MMMBackground 12 from photutils.datasets import make_noise_image /opt/conda/envs/notebooks_env/lib/python3.6/site-packages/photutils/__init__.py in <module> 17 18 if not _ASTROPY_SETUP_: # noqa ---> 19 from .aperture import * # noqa 20 from .background import * # noqa 21 from .centroids import * # noqa /opt/conda/envs/notebooks_env/lib/python3.6/site-packages/photutils/aperture/__init__.py in <module> 5 6 from .bounding_box import * # noqa ----> 7 from .circle import * # noqa 8 from .core import * # noqa 9 from .ellipse import * # noqa /opt/conda/envs/notebooks_env/lib/python3.6/site-packages/photutils/aperture/circle.py in <module> 13 from .core import PixelAperture, SkyAperture 14 from .mask import ApertureMask ---> 15 from ..geometry import circular_overlap_grid 16 17 __all__ = ['CircularMaskMixin', 'CircularAperture', 'CircularAnnulus', /opt/conda/envs/notebooks_env/lib/python3.6/site-packages/photutils/geometry/__init__.py in <module> 4 """ 5 ----> 6 from .circular_overlap import * # noqa 7 from .elliptical_overlap import * # noqa 8 from .rectangular_overlap import * # noqa __init__.pxd in init photutils.geometry.circular_overlap() ValueError: numpy.ufunc size changed, may indicate binary incompatibility. Expected 216 from C header, got 192 from PyObject
The simulated image will be made using a spatially-dependent JWST PSF generated by
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:
We start by loading into FITS HDU a pre-generated grid of JWST PSFs for NIRCam A1 detector in F090W filter.
url = 'https://stsci.box.com/shared/static/6h8wsr2ts0t24s79cxmnyk8bldt0vv3i.fits' hdu = fits.open(url)
Astropy has saved a cached copy of this file under your local cache directory for Astropy; default is
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.
data = hdu.data plt.figure(figsize=(9, 9)) imshow_norm(data, interval=PercentileInterval(99.), stretch=SqrtStretch())
psfmodel = to_griddedpsfmodel(hdu, ext=0)
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.
oversampling attribute is the PSF oversampling factor.
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.
This model contains a 4x4 grid of reference PSFs.
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.
fovpixel = psfmodel.meta['fovpixel'] print(fovpixel, 'is', fovpixel)
We can use the
webbpsf.gridded_library.display_psf_grid function to visualize the PSF grid.
Now 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.
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, size=nstars) yy = np.random.uniform(low=0, high=shape, size=nstars) zz = np.random.uniform(low=0, high=flux_max, size=nstars)
Now we'll evaluate the model at these positions and fluxes.
eval_xshape = np.int(np.ceil(psfmodel.data.shape / psfmodel.oversampling)) eval_yshape = np.int(np.ceil(psfmodel.data.shape / 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: x1 = shape if y1 > shape: y1 = shape 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)
Let's add some noise to the image.
noise = make_noise_image(data.shape, 'gaussian', mean=0, stddev=2, random_state=123) data += noise
Let's display the simulated JWST NIRCam A1 detector image in F090W filter.
plt.figure(figsize=(9, 9)) imshow_norm(data, interval=PercentileInterval(99.), stretch=SqrtStretch())
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.
init_tbl = Table() init_tbl['x_0'] = xx.astype(int) init_tbl['y_0'] = yy.astype(int) init_tbl['flux_0'] = zz.astype(int)
Then, we define the parameters to initialize a
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.)
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.
tbl = phot(data, init_guesses=init_tbl)
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
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'
We can also view the residual image of the best-fit model PSF image subtracted from the data.
diff = phot.get_residual_image()
plt.figure(figsize=(9, 9)) imshow_norm(diff, interval=PercentileInterval(99.), stretch=SqrtStretch()) plt.colorbar()
The residual image is essentially noise, indicating good PSF model fits to the data.
Author: Data Analysis Tools Branch, STScI
Updated on: Consult GitHub repository for edit history