.. _psf-photometry: PSF Photometry (`photutils.psf`) ================================ The `photutils.psf` subpackage contains tools for model-fitting photometry, often called "PSF photometry". .. _psf-terminology: Terminology ----------- Different astronomy subfields use the terms "PSF", "PRF", or related terms somewhat differently, especially when colloquial usage is taken into account. This package aims to be at the very least internally consistent, following the definitions described here. We take the Point Spread Function (PSF), or instrumental Point Spread Function (iPSF) to be the infinite-resolution and infinite-signal-to-noise flux distribution from a point source on the detector, after passing through optics, dust, atmosphere, etc. By contrast, the function describing the responsivity variations across individual *pixels* is the Pixel Response Function (sometimes called "PRF", but that acronym is not used here for reasons that will soon be apparent). The convolution of the PSF and pixel response function, when discretized onto the detector (i.e., a rectilinear CCD grid), is the effective PSF (ePSF) or Point Response Function (PRF) (this latter terminology is the definition used by `Spitzer `_). In many cases the PSF/ePSF/PRF distinction is unimportant, and the ePSF/PRF are simply called the "PSF", but the distinction can be critical when dealing carefully with undersampled data or detectors with significant intra-pixel sensitivity variations. For a more detailed description of this formalism, see `Anderson & King 2000 `_. All this said, in colloquial usage "PSF photometry" sometimes refers to the more general task of model-fitting photometry (with the effects of the PSF either implicitly or explicitly included in the models), regardless of exactly what kind of model is actually being fit. For brevity (e.g., ``photutils.psf``), we use "PSF photometry" in this way, as a shorthand for the general approach. PSF Photometry -------------- Photutils provides a modular set of tools to perform PSF photometry for different science cases. The tools are implemented as classes that perform various subtasks of PSF photometry. High-level classes are also provided to connect these pieces together. The two main PSF-photometry classes are `~photutils.psf.PSFPhotometry` and `~photutils.psf.IterativePSFPhotometry`. `~photutils.psf.PSFPhotometry` provides the framework for a flexible PSF photometry workflow that can find sources in an image, optionally group overlapping sources, fit the PSF model to the sources, and subtract the fit PSF models from the image. `~photutils.psf.IterativePSFPhotometry` is an iterative version of `~photutils.psf.PSFPhotometry` where after the fit sources are subtracted, the process repeats until no additional sources are detected or a maximum number of iterations has been reached. When used with the `~photutils.detection.DAOStarFinder`, `~photutils.psf.IterativePSFPhotometry` is essentially an implementation of the DAOPHOT algorithm described by Stetson in his `seminal paper `_ for crowded-field stellar photometry. The star-finding step is controlled by the ``finder`` keyword, where one inputs a callable function or class instance. Typically, this would be one of the star-detection classes implemented in the `photutils.detection` subpackage, e.g., `~photutils.detection.DAOStarFinder`, `~photutils.detection.IRAFStarFinder`, or `~photutils.detection.StarFinder`. After finding sources, one can optionally apply a clustering algorithm to group overlapping sources using the ``grouper`` keyword. Usually, groups are formed by a distance criterion, which is the case of the grouping algorithm proposed by Stetson. Stars that grouped are fit simultaneously. The reason behind the construction of groups and not fitting all stars simultaneously is illustrated as follows: imagine that one would like to fit 300 stars and the model for each star has three parameters to be fitted. If one constructs a single model to fit the 300 stars simultaneously, then the optimization algorithm will have to search for the solution in a 900-dimensional space, which is computationally expensive and error-prone. Having smaller groups of stars effectively reduces the dimension of the parameter space, which facilitates the optimization process. For more details see :ref:`psf-grouping`. The local background around each source can optionally be subtracted using the ``localbkg_estimator`` keyword. This keyword accepts a `~photutils.background.LocalBackground` instance that estimates the local statistics in a circular annulus aperture centered on each source. The size of the annulus and the statistic function can be configured in `~photutils.background.LocalBackground`. The next step is to fit the sources and/or groups. This task is performed using an astropy fitter, for example `~astropy.modeling.fitting.LevMarLSQFitter`, input via the ``fitter`` keyword. The shape of the region to be fitted can be configured using the ``fit_shape`` parameter. In general, ``fit_shape`` should be set to a small size (e.g., (5, 5)) that covers the central star region with the highest flux signal-to-noise. The initial positions are derived from the ``finder`` algorithm. The initial flux values for the fit are derived from measuring the flux in a circular aperture with radius ``aperture_radius``. The initial positions and fluxes can be alternatively input in a table via the ``init_params`` keyword when calling the class. After sources are fitted, a model image of the fit sources or a residual image can be generated using the :meth:`~photutils.psf.PSFPhotometry.make_model_image` and :meth:`~photutils.psf.PSFPhotometry.make_residual_image` methods, respectively. For `~photutils.psf.IterativePSFPhotometry`, the above steps can be repeated until no additional sources are detected (or until a maximum number of iterations). The `~photutils.psf.PSFPhotometry` and `~photutils.psf.IterativePSFPhotometry` classes provide the structure in which the PSF-fitting steps described above are performed, but all the stages can be turned on or off or replaced with different implementations as the user desires. This makes the tools very flexible. One can also bypass several of the steps by directly inputting to ``init_params`` an astropy table containing the initial parameters for the source centers, fluxes, group identifiers, and local backgrounds. This is also useful if one is interested in fitting only one or a few sources in an image. Example Usage ------------- Let's start with a simple example using simulated stars whose PSF is assumed to be Gaussian. We'll create a synthetic image using tools provided by the :ref:`photutils.datasets ` module: .. doctest-requires:: scipy >>> import numpy as np >>> from photutils.datasets import make_test_psf_data, make_noise_image >>> from photutils.psf import IntegratedGaussianPRF >>> psf_model = IntegratedGaussianPRF(flux=1, sigma=2.7 / 2.35) >>> psf_shape = (9, 9) >>> nsources = 10 >>> shape = (101, 101) >>> data, true_params = make_test_psf_data(shape, psf_model, psf_shape, ... nsources, flux_range=(500, 700), ... min_separation=10, seed=0) >>> noise = make_noise_image(data.shape, mean=0, stddev=1, seed=0) >>> data += noise >>> error = np.abs(noise) Let's plot the image: .. plot:: import matplotlib.pyplot as plt from photutils.datasets import make_noise_image, make_test_psf_data from photutils.psf import IntegratedGaussianPRF psf_model = IntegratedGaussianPRF(flux=1, sigma=2.7 / 2.35) psf_shape = (9, 9) nsources = 10 shape = (101, 101) data, true_params = make_test_psf_data(shape, psf_model, psf_shape, nsources, flux_range=(500, 700), min_separation=10, seed=0) noise = make_noise_image(data.shape, mean=0, stddev=1, seed=0) data += noise plt.imshow(data, origin='lower') plt.title('Simulated Data') plt.colorbar() Fitting multiple stars ^^^^^^^^^^^^^^^^^^^^^^ Now let's use `~photutils.psf.PSFPhotometry` to perform PSF photometry on the stars in this image. Note that the input image must be background-subtracted prior to using the photometry classes. See :ref:`background` for tools to subtract a global background from an image. This is not needed for our synthetic image because it does not include background. We'll use the `~photutils.detection.DAOStarFinder` class for source detection. We'll estimate the initial fluxes of each source using a circular aperture with a radius 4 pixels. The central 5x5 pixel region of each star will be fit using an `~photutils.psf.IntegratedGaussianPRF` PSF model. First, let's create an instance of the `~photutils.psf.PSFPhotometry` class: .. doctest-requires:: scipy >>> from photutils.detection import DAOStarFinder >>> from photutils.psf import PSFPhotometry >>> psf_model = IntegratedGaussianPRF(flux=1, sigma=2.7 / 2.35) >>> fit_shape = (5, 5) >>> finder = DAOStarFinder(6.0, 2.0) >>> psfphot = PSFPhotometry(psf_model, fit_shape, finder=finder, ... aperture_radius=4) To perform the PSF fitting, we then call the class instance on the data array, and optionally an error and mask array. A `~astropy.nddata.NDData` object holding the data, error, and mask arrays can also be input into the ``data`` parameter. Note that all non-finite (e.g., NaN or inf) data values are automatically masked. Here we input the data and error arrays: .. doctest-requires:: scipy >>> phot = psfphot(data, error=error) A table of initial PSF model parameter values can also be input when calling the class instance. An example of that is shown later. Equivalently, one can input an `~astropy.nddata.NDData` object with any uncertainty object that can be converted to standard-deviation errors: .. doctest-skip:: >>> from astropy.nddata import NDData, StdDevUncertainty >>> uncertainty = StdDevUncertainty(error) >>> nddata = NDData(data, uncertainty=uncertainty) >>> phot2 = psfphot(nddata) The result is an astropy `~astropy.table.Table` with columns for the source and group identification numbers, the x, y, and flux initial, fit, and error values, local background, number of unmasked pixels fit, the group size, quality-of-fit metrics, and flags. See the `~photutils.psf.PSFPhotometry` documentation for descriptions of the output columns. The full table cannot be shown here as it has many columns, but let's print the source ID along with the fit x, y, and flux values: .. doctest-requires:: scipy >>> phot['x_fit'].info.format = '.4f' # optional format >>> phot['y_fit'].info.format = '.4f' >>> phot['flux_fit'].info.format = '.4f' >>> print(phot[('id', 'x_fit', 'y_fit', 'flux_fit')]) # doctest: +FLOAT_CMP id x_fit y_fit flux_fit --- ------- ------- -------- 1 54.5641 7.7650 645.4552 2 29.0866 25.6111 553.0376 3 79.6285 28.7488 665.8715 4 63.2344 48.6406 626.7983 5 88.8849 54.1203 681.9908 6 79.8762 61.1380 686.3628 7 90.9606 72.0860 610.5911 8 7.8021 78.5730 509.3270 9 5.5349 89.8869 503.9691 10 71.8411 90.5841 621.1757 Let's create the residual image: .. doctest-requires:: scipy >>> resid = psfphot.make_residual_image(data, (9, 9)) and plot it: .. plot:: import matplotlib.pyplot as plt import numpy as np from astropy.visualization import simple_norm from photutils.datasets import make_noise_image, make_test_psf_data from photutils.detection import DAOStarFinder from photutils.psf import IntegratedGaussianPRF, PSFPhotometry psf_model = IntegratedGaussianPRF(flux=1, sigma=2.7 / 2.35) psf_shape = (9, 9) nsources = 10 shape = (101, 101) data, true_params = make_test_psf_data(shape, psf_model, psf_shape, nsources, flux_range=(500, 700), min_separation=10, seed=0) noise = make_noise_image(data.shape, mean=0, stddev=1, seed=0) data += noise error = np.abs(noise) psf_model = IntegratedGaussianPRF(flux=1, sigma=2.7 / 2.35) fit_shape = (5, 5) finder = DAOStarFinder(6.0, 2.0) psfphot = PSFPhotometry(psf_model, fit_shape, finder=finder, aperture_radius=4) phot = psfphot(data, error=error) resid = psfphot.make_residual_image(data, (9, 9)) fig, ax = plt.subplots(nrows=1, ncols=3, figsize=(15, 5)) norm = simple_norm(data, 'sqrt', percent=99) ax[0].imshow(data, origin='lower', norm=norm) ax[1].imshow(data - resid, origin='lower', norm=norm) im = ax[2].imshow(resid, origin='lower') ax[0].set_title('Data') ax[1].set_title('Model') ax[2].set_title('Residual Image') plt.tight_layout() The residual image looks like noise, indicating good fits to the sources. Further details about the PSF fitting can be obtained from attributes on the `~photutils.psf.PSFPhotometry` instance. For example, the results from the ``finder`` instance called during PSF fitting can be accessed using the ``finder_results`` attribute (the ``finder`` returns an astropy table): .. doctest-requires:: scipy >>> psfphot.finder_results['xcentroid'].info.format = '.4f' # optional format >>> psfphot.finder_results['ycentroid'].info.format = '.4f' # optional format >>> psfphot.finder_results['sharpness'].info.format = '.4f' # optional format >>> psfphot.finder_results['peak'].info.format = '.4f' >>> psfphot.finder_results['flux'].info.format = '.4f' >>> psfphot.finder_results['mag'].info.format = '.4f' >>> print(psfphot.finder_results) # doctest: +FLOAT_CMP id xcentroid ycentroid sharpness ... sky peak flux mag --- --------- --------- --------- ... --- ------- ------- ------- 1 54.5300 7.7508 0.5996 ... 0.0 67.0314 8.7012 -2.3490 2 29.0926 25.5993 0.5952 ... 0.0 58.7504 7.7484 -2.2230 3 79.6189 28.7515 0.5956 ... 0.0 70.4621 9.2351 -2.4136 4 63.2472 48.6151 0.5817 ... 0.0 64.9101 8.5701 -2.3325 5 88.8827 54.1299 0.5954 ... 0.0 75.8880 10.3243 -2.5347 6 79.8730 61.1214 0.6204 ... 0.0 78.0913 10.3346 -2.5357 7 90.9621 72.0802 0.6162 ... 0.0 69.1417 9.2325 -2.4133 8 7.7917 78.5454 0.5971 ... 0.0 52.7325 6.6840 -2.0626 9 5.5845 89.8645 0.5721 ... 0.0 50.5544 6.5387 -2.0387 10 71.8284 90.5625 0.6038 ... 0.0 65.7975 8.5373 -2.3283 The ``fit_results`` attribute contains a dictionary with a wealth of detailed information, including the fit models and any information returned from the ``fitter`` for each source: .. doctest-requires:: scipy >>> psfphot.fit_results.keys() dict_keys(['local_bkg', 'init_params', 'fit_infos', 'fit_param_errs', 'fit_error_indices', 'npixfit', 'nmodels', 'psfcenter_indices', 'fit_residuals']) As an example, let's print the covariance matrix of the fit parameters for the first source (note that not all astropy fitters will return a covariance matrix): .. doctest-skip:: >>> psfphot.fit_results['fit_infos'][0]['param_cov'] # doctest: +FLOAT_CMP array([[ 7.27034774e-01, 8.86845334e-04, 3.98593038e-03], [ 8.86845334e-04, 2.92871525e-06, -6.36805464e-07], [ 3.98593038e-03, -6.36805464e-07, 4.29520185e-05]]) Fitting a single source ^^^^^^^^^^^^^^^^^^^^^^^ In some cases, one may want to fit only a single source (or few sources) in an image. We can do that by defining a table of the sources that we want to fit. For this example, let's fit the single star at ``(x, y) = (42, 36)``. We first define a table with this position and then pass that table into the ``init_params`` keyword when calling the PSF photometry class on the data: .. doctest-requires:: scipy >>> from astropy.table import QTable >>> init_params = QTable() >>> init_params['x'] = [63] >>> init_params['y'] = [49] >>> phot = psfphot(data, error=error, init_params=init_params) The PSF photometry class allows for flexible input column names using a heuristic to identify the x, y, and flux columns. See `~photutils.psf.PSFPhotometry` for more details. The output table contains only the fit results for the input source: .. doctest-requires:: scipy >>> phot['x_fit'].info.format = '.4f' # optional format >>> phot['y_fit'].info.format = '.4f' >>> phot['flux_fit'].info.format = '.4f' >>> print(phot[('id', 'x_fit', 'y_fit', 'flux_fit')]) # doctest: +FLOAT_CMP id x_fit y_fit flux_fit --- ------- ------- -------- 1 63.2344 48.6406 626.7983 Finally, let's show the residual image. The red circular aperture shows the location of the star that was fit and subtracted. .. plot:: import matplotlib.pyplot as plt import numpy as np from astropy.table import QTable from photutils.aperture import CircularAperture from photutils.datasets import make_noise_image, make_test_psf_data from photutils.detection import DAOStarFinder from photutils.psf import IntegratedGaussianPRF, PSFPhotometry psf_model = IntegratedGaussianPRF(flux=1, sigma=2.7 / 2.35) psf_shape = (9, 9) nsources = 10 shape = (101, 101) data, true_params = make_test_psf_data(shape, psf_model, psf_shape, nsources, flux_range=(500, 700), min_separation=10, seed=0) noise = make_noise_image(data.shape, mean=0, stddev=1, seed=0) data += noise error = np.abs(noise) psf_model = IntegratedGaussianPRF(flux=1, sigma=2.7 / 2.35) fit_shape = (5, 5) finder = DAOStarFinder(6.0, 2.0) psfphot = PSFPhotometry(psf_model, fit_shape, finder=finder, aperture_radius=4) init_params = QTable() init_params['x'] = [63] init_params['y'] = [49] phot = psfphot(data, error=error, init_params=init_params) resid = psfphot.make_residual_image(data, (9, 9)) plt.imshow(resid, origin='lower') resid = psfphot.make_residual_image(data, (9, 9)) aper = CircularAperture(zip(init_params['x'], init_params['y']), r=4) plt.imshow(resid, origin='lower') aper.plot(color='red') plt.title('Residual Image') plt.colorbar() Forced Photometry ^^^^^^^^^^^^^^^^^ In general, the three parameters fit for each source are the x and y positions and the flux. However, the astropy modeling and fitting framework allows any of these parameters to be fixed during the fitting. Let's say you want to fix the (x, y) position for each source. You can do that by setting the ``fixed`` attribute on the model parameters: .. doctest-requires:: scipy >>> psf_model2 = IntegratedGaussianPRF(flux=1, sigma=2.7 / 2.35) >>> psf_model2.x_0.fixed = True >>> psf_model2.y_0.fixed = True >>> psf_model2.fixed {'flux': False, 'x_0': True, 'y_0': True, 'sigma': True} Now when the model is fit, the flux will be varied but, the (x, y) position will be fixed at its initial position for every source. Let's just fit a single source (defined in ``init_params``): .. doctest-requires:: scipy >>> psfphot = PSFPhotometry(psf_model2, fit_shape, finder=finder, ... aperture_radius=4) >>> phot = psfphot(data, error=error, init_params=init_params) The output table shows that the (x, y) position was unchanged, with the fit values being identical to the initial values. However, the flux was fit: .. doctest-requires:: scipy >>> phot['flux_init'].info.format = '.4f' # optional format >>> phot['flux_fit'].info.format = '.4f' >>> print(phot[('id', 'x_init', 'y_init', 'flux_init', 'x_fit', ... 'y_fit', 'flux_fit')]) # doctest: +FLOAT_CMP id x_init y_init flux_init x_fit y_fit flux_fit --- ------ ------ --------- ----- ----- -------- 1 63 49 619.5049 63.0 49.0 556.8463 Source Grouping ^^^^^^^^^^^^^^^ Source grouping is an optional feature. To turn it on, create a `~photutils.psf.SourceGrouper` instance and input it via the ``grouper`` keyword. Here we'll group stars that are within 20 pixels of each other: .. doctest-requires:: scipy, sklearn >>> from photutils.psf import SourceGrouper >>> grouper = SourceGrouper(min_separation=20) >>> psfphot = PSFPhotometry(psf_model, fit_shape, finder=finder, ... grouper=grouper, aperture_radius=4) >>> phot = psfphot(data, error=error) The ``group_id`` column shows that six groups were identified (each with two stars). The stars in each group were simultaneously fit. .. doctest-requires:: scipy, sklearn >>> print(phot[('id', 'group_id', 'group_size')]) id group_id group_size --- -------- ---------- 1 1 1 2 2 1 3 3 1 4 4 1 5 5 3 6 5 3 7 5 3 8 6 2 9 6 2 10 7 1 Care should be taken in defining the star groups. As noted above, simultaneously fitting very large star groups is computationally expensive and error-prone. A warning will be raised if the number of sources in a group exceeds 25. Local Background Subtraction ^^^^^^^^^^^^^^^^^^^^^^^^^^^^ To subtract a local background from each source, define a `~photutils.background.LocalBackground` instance and input it via the ``localbkg_estimator`` keyword. Here we'll use an annulus with an inner and outer radius of 5 and 10 pixels, respectively, with the `~photutils.background.MMMBackground` statistic (with its default sigma clipping): .. doctest-requires:: scipy, sklearn >>> from photutils.background import LocalBackground, MMMBackground >>> bkgstat = MMMBackground() >>> localbkg_estimator = LocalBackground(5, 10, bkgstat) >>> finder = DAOStarFinder(10.0, 2.0) >>> psfphot = PSFPhotometry(psf_model, fit_shape, finder=finder, ... grouper=grouper, aperture_radius=4, ... localbkg_estimator=localbkg_estimator) >>> phot = psfphot(data, error=error) The local background values are output in the table: .. doctest-requires:: scipy, sklearn >>> phot['local_bkg'].info.format = '.4f' # optional format >>> print(phot[('id', 'local_bkg')]) # doctest: +FLOAT_CMP id local_bkg --- --------- 1 -0.0840 2 0.1784 3 0.2593 4 -0.0455 5 0.3212 6 -0.0836 7 -0.1130 8 -0.2482 9 0.0315 10 0.3926 The local background values can also be input directly using the ``init_params`` keyword. Iterative PSF Photometry ^^^^^^^^^^^^^^^^^^^^^^^^ Now let's use the `~photutils.psf.IterativePSFPhotometry` class to iteratively fit the stars in the image. This class is useful for crowded fields where faint stars are very close to bright stars. The faint stars may not be detected until after the bright stars are subtracted. For this simple example, let's input a table of three stars for the first fit iteration. Subsequent iterations will use the ``finder`` to find additional stars: .. doctest-requires:: scipy >>> from photutils.background import LocalBackground, MMMBackground >>> from photutils.psf import IterativePSFPhotometry >>> fit_shape = (5, 5) >>> finder = DAOStarFinder(10.0, 2.0) >>> bkgstat = MMMBackground() >>> localbkg_estimator = LocalBackground(5, 10, bkgstat) >>> init_params = QTable() >>> init_params['x'] = [54, 29, 80] >>> init_params['y'] = [8, 26, 29] >>> psfphot2 = IterativePSFPhotometry(psf_model, fit_shape, finder=finder, ... localbkg_estimator=localbkg_estimator, ... aperture_radius=4) >>> phot = psfphot2(data, error=error, init_params=init_params) The table output from `~photutils.psf.IterativePSFPhotometry` contains a column called ``iter_detected`` that returns the fit iteration in which the source was detected: .. doctest-requires:: scipy >>> phot['x_fit'].info.format = '.4f' # optional format >>> phot['y_fit'].info.format = '.4f' >>> phot['flux_fit'].info.format = '.4f' >>> print(phot[('id', 'iter_detected', 'x_fit', 'y_fit', 'flux_fit')]) # doctest: +FLOAT_CMP id iter_detected x_fit y_fit flux_fit --- ------------- ------- ------- -------- 1 1 54.5646 7.7648 645.7164 2 1 29.0884 25.6093 550.5295 3 1 79.6278 28.7481 660.1429 4 2 63.2344 48.6411 627.4414 5 2 88.8859 54.1203 676.3059 6 2 79.8764 61.1359 688.1962 7 2 90.9631 72.0879 612.4822 8 2 7.8259 78.5856 516.4237 9 2 5.5348 89.8868 503.4775 10 2 71.8491 90.5828 616.2827 References ---------- `Spitzer PSF vs. PRF `_ `The Kepler Pixel Response Function `_ `Stetson (1987 PASP 99, 191) `_ `Anderson and King (2000 PASP 112, 1360) `_ Reference/API ------------- .. automodapi:: photutils.psf :no-heading: