add function to create GriddedPSFModel from EPSFs

CHANGES.rst

@@ -29,6 +29,8 @@ New Features
 
   - Added a ``GriddedPSFModel`` ``fill_value`` attribute, [#1583]
 
+  - Added function to make GriddedPSFModel from EPSFs. [#1584]
+
 Bug Fixes
 ^^^^^^^^^
 

photutils/psf/tests/test_utils.py

@@ -7,15 +7,19 @@
 import pytest
 from astropy.convolution.utils import discretize_model
 from astropy.modeling.models import Gaussian2D
+from astropy.nddata import NDData
 from astropy.table import Table
 from astropy.utils.exceptions import AstropyDeprecationWarning
 from numpy.testing import assert_allclose
 
+from photutils.detection import find_peaks
+from photutils import datasets
+from photutils.psf import EPSFBuilder, extract_stars
 from photutils.psf.groupstars import DAOGroup
 from photutils.psf.models import IntegratedGaussianPRF
 from photutils.psf.photometry_depr import BasicPSFPhotometry
-from photutils.psf.utils import (get_grouped_psf_model, prepare_psf_model,
-                                 subtract_psf)
+from photutils.psf.utils import (get_grouped_psf_model, grid_from_epsfs,
+                                 prepare_psf_model, subtract_psf)
 from photutils.utils._optional_deps import HAS_SCIPY
 
 PSF_SIZE = 11
@@ -260,3 +264,104 @@ def test_subtract_psf():
             posflux.rename_column(n, n.split('_')[0] + '_fit')
         residuals = subtract_psf(image, psf, posflux)
         assert np.max(np.abs(residuals)) < 0.0052
+
+
+@pytest.mark.remote_data
+class TestGridFromEPSFs:
+    """Tests for `photutils.psf.utils.grid_from_epsfs`."""
+
+    def setup_class(self, cutout_size=25):
+        # make a set of 4 EPSF models
+
+        self.cutout_size = cutout_size
+
+        # make simulated image
+        hdu = datasets.load_simulated_hst_star_image()
+        data = hdu.data
+
+        # break up the image into four quadrants
+        q1 = data[0:500, 0:500]
+        q2 = data[0:500, 500:1000]
+        q3 = data[500:1000, 0:500]
+        q4 = data[500:1000, 500:1000]
+
+        # select some starts from each quadrant to use to build the epsf
+        quad_stars = {'q1': {'data': q1, 'fiducial': (0., 0.), 'epsf': None},
+                      'q2': {'data': q2, 'fiducial': (1000., 1000.), 'epsf': None},
+                      'q3': {'data': q3, 'fiducial': (1000., 0.), 'epsf': None},
+                      'q4': {'data': q4, 'fiducial': (0., 1000.), 'epsf': None}}
+
+        for q in ['q1', 'q2', 'q3', 'q4']:
+            quad_data = quad_stars[q]['data']
+            peaks_tbl = find_peaks(quad_data, threshold=500.)
+
+            # filter out sources near edge
+            size = cutout_size
+            hsize = (size - 1) / 2
+            x = peaks_tbl['x_peak']
+            y = peaks_tbl['y_peak']
+            mask = ((x > hsize) & (x < (quad_data.shape[1] - 1 - hsize))
+                    & (y > hsize) & (y < (quad_data.shape[0] - 1 - hsize)))
+
+            stars_tbl = Table()
+            stars_tbl['x'] = peaks_tbl['x_peak'][mask]
+            stars_tbl['y'] = peaks_tbl['y_peak'][mask]
+
+            stars = extract_stars(NDData(quad_data), stars_tbl,
+                                  size=cutout_size)
+
+            epsf_builder = EPSFBuilder(oversampling=4, maxiters=3,
+                                       progress_bar=False)
+            epsf, fitted_stars = epsf_builder(stars)
+
+            # set x_0, y_0 to fiducial point
+            epsf.y_0 = quad_stars[q]['fiducial'][0]
+            epsf.x_0 = quad_stars[q]['fiducial'][1]
+
+            quad_stars[q]['epsf'] = epsf
+
+        self.epsfs = [quad_stars[x]['epsf'] for x in quad_stars]
+        self.grid_xypos = [quad_stars[x]['fiducial'] for x in quad_stars]
+
+    def test_basic_test_grid_from_epsfs(self):
+
+        psf_grid = grid_from_epsfs(self.epsfs)
+
+        assert np.all(psf_grid.oversampling == self.epsfs[0].oversampling)
+        assert psf_grid.data.shape == (4, psf_grid.oversampling * 25 + 1,
+                                       psf_grid.oversampling * 25 + 1)
+
+    def test_grid_xypos(self):
+        """Test both options for setting PSF locations"""
+
+        # default option x_0 and y_0s on input EPSFs
+        psf_grid = grid_from_epsfs(self.epsfs)
+
+        assert psf_grid.meta['grid_xypos'] == [(0.0, 0.0), (1000.0, 1000.0),
+                                               (0.0, 1000.0), (1000.0, 0.0)]
+
+        # or pass in a list
+        grid_xypos = [(250.0, 250.0), (750.0, 750.0),
+                      (250.0, 750.0), (750.0, 250.0)]
+
+        psf_grid = grid_from_epsfs(self.epsfs, grid_xypos=grid_xypos)
+        assert psf_grid.meta['grid_xypos'] == grid_xypos
+
+    def test_meta(self):
+        """Test the option for setting 'meta'"""
+
+        keys = ['grid_xypos', 'oversampling', 'fill_value']
+
+        # when 'meta' isn't provided, there should be just three keys
+        psf_grid = grid_from_epsfs(self.epsfs)
+        assert list(psf_grid.meta.keys()) == keys
+
+        # when meta is provided, those new keys should exist and anything
+        # in the list above should be overwritten
+        meta = {'grid_xypos': 0.0, 'oversampling': 0.0,
+                'fill_value': -999, 'extra_key': 'extra'}
+        psf_grid = grid_from_epsfs(self.epsfs, meta=meta)
+        assert list(psf_grid.meta.keys()) == keys + ['extra_key']
+        assert psf_grid.meta['grid_xypos'].sort() == self.grid_xypos.sort()
+        assert psf_grid.meta['oversampling'] == 4
+        assert psf_grid.meta['fill_value'] == 0.0

photutils/psf/utils.py

@@ -5,11 +5,14 @@
 
 import numpy as np
 from astropy.modeling.models import Const2D, Identity, Shift
+from astropy.nddata import NDData
 from astropy.nddata.utils import add_array, extract_array
 from astropy.table import QTable
 from astropy.utils.decorators import deprecated
+from photutils.psf.models import EPSFModel, GriddedPSFModel
 
-__all__ = ['prepare_psf_model', 'get_grouped_psf_model', 'subtract_psf']
+__all__ = ['prepare_psf_model', 'get_grouped_psf_model', 'subtract_psf',
+           'grid_from_epsfs']
 
 
 class _InverseShift(Shift):
@@ -335,3 +338,156 @@ def subtract_psf(data, psf, posflux, *, subshape=None):
             subbeddata = add_array(subbeddata, -psf(x, y), (y_0, x_0))
 
     return subbeddata
+
+
+def grid_from_epsfs(epsfs, grid_xypos=None, meta=None):
+
+    """
+    Create a GriddedPSFModel from a list of EPSFModels.
+
+    Given a list of EPSFModels, this function will return a GriddedPSFModel.
+    The fiducial points for each input EPSFModel can either be set on each
+    individual model by setting the 'x_0' and 'y_0' attributes, or provided as
+    a list of tuples (``grid_xypos``). If a ``grid_xypos`` list is provided, it
+    must match the length of input EPSFs. In either case, the fiducial points
+    must be on a grid.
+
+    Optionally, a ``meta`` dictionary may be provided for the
+    output GriddedPSFModel. If this dictionary contains the keys 'grid_xypos',
+    'oversampling', or 'fill_value', they will be overridden.
+
+    Note: If set on the input EPSFModel (x_0, y_0), then ``origin`` must be the
+    same for each input EPSF. Additionally data units and dimensions must be
+    for each input EPSF, and values for ``flux`` and ``oversampling``, and
+    ``fill_value`` must match as well.
+
+    Parameters
+    ----------
+    epsfs : list of `photutils.psf.models.EPSFModel`
+        A list of EPSFModels representing the individual PSFs.
+    grid_xypos : list, optional
+        A list of fiducial points (x_0, y_0) for each PSF.
+        If not provided, the x_0 and y_0 of each input EPSF will be considered
+        the fiducial point for that PSF. Default is None.
+    meta : dict, optional
+        Additional metadata for the GriddedPSFModel. Note that, if they exist
+        in the supplied ``meta``, any values under the keys ``grid_xypos`` ,
+        ``oversampling``, or ``fill_value`` will be overridden. Default is None.
+
+    Returns
+    -------
+    GriddedPSFModel: `photutils.psf.models.GriddedPSFModel`
+        The gridded PSF model created from the input EPSFs.
+    """
+
+    x_0s = []  # optional, to store fiducial from input if `grid_xypos` is None
+    y_0s = []
+    data_arrs = []
+    oversampling = None
+    fill_value = None
+    dat_unit = None
+    origin = None
+    flux = None
+
+    # make sure, if provided, that ``grid_xypos`` is the same length as ``epsfs``
+    if grid_xypos is not None:
+        if len(grid_xypos) != len(epsfs):
+            raise ValueError('``grid_xypos`` must be the same length as ``epsfs``.')
+
+    # loop over input once
+    for i, epsf in enumerate(epsfs):
+
+        # check input type
+        if not isinstance(epsf, EPSFModel):
+            raise ValueError('All input `epsfs` must be of type '
+                             '`photutils.psf.models.EPSFModel`.')
+
+        # get data array from EPSF
+        data_arrs.append(epsf.data)
+
+        if i == 0:
+            # EPSFModel allows a tuple for oversampling factor in x, y,
+            # but GriddedPSFModel requires it to be a single scalar value.
+            # Keep this condition for now by checking that x and y match
+            if np.isscalar(epsf.oversampling):
+                oversampling = epsf.oversampling
+            else:
+                if epsf.oversampling[0] != epsf.oversampling[1]:
+                    raise ValueError('Oversampling must be the same in x and y.')
+                oversampling = epsf.oversampling[0]
+
+            # same for fill value and flux, grid will have a single value
+            # so it should be the same for all input, and error if not.
+            fill_value = epsf.fill_value
+
+            # check that origins are the same
+            if grid_xypos is None:
+                origin = epsf.origin
+
+            flux = epsf.flux
+
+            # if theres a unit, those should also all be the same
+            try:
+                dat_unit = epsf.data.unit
+            except AttributeError:
+                pass  # just keep as None
+
+        else:
+            if np.isscalar(epsf.oversampling):
+                if epsf.oversampling != oversampling:
+                    raise ValueError('All input EPSFModels must have the same '
+                                     'value for ``oversampling``.')
+                else:
+                    if epsf.oversampling[0] != epsf.oversampling[1] != oversampling:
+                        raise ValueError('All input EPSFModels must have the '
+                                         'same value for ``oversampling``.')
+
+            if epsf.fill_value != fill_value:
+                raise ValueError('All input EPSFModels must have the same value '
+                                 'for ``fill_value``.')
+
+            if epsf.data.ndim != data_arrs[0].ndim:
+                raise ValueError('All input EPSFModels must have data with the '
+                                 'same dimensions.')
+
+            try:
+                unitt = epsf.data_unit
+                if unitt != dat_unit:
+                    raise ValueError('All input data must have the same unit.')
+            except AttributeError:
+                if dat_unit is not None:
+                    raise ValueError('All input data must have the same unit.')
+
+            if epsf.flux != flux:
+                raise ValueError('All input EPSFModels must have the same value '
+                                 'for ``flux``.')
+
+        if grid_xypos is None:  # get gridxy_pos from x_0, y_0 if not provided
+            x_0s.append(epsf.x_0.value)
+            y_0s.append(epsf.y_0.value)
+
+            # also check that origin is the same, if using x_0s and y_0s from input
+            if epsf.origin != origin:
+                raise ValueError('If using ``x_0``, ``y_0`` as fiducial point,'
+                                 '``origin`` must match for each input EPSF.')
+
+    # if not supplied, use from x_0, y_0 of input EPSFs as fiducuals
+    # these are checked when GriddedPSFModel is created to make sure they
+    # are actually on a grid.
+    if grid_xypos is None:
+        grid_xypos = list(zip(x_0s, y_0s))
+
+    data_cube = np.stack(data_arrs, axis=0)
+
+    if meta is None:
+        meta = {}
+    # add required keywords to meta
+    meta['grid_xypos'] = grid_xypos
+    meta['oversampling'] = oversampling
+    meta['fill_value'] = fill_value
+
+    data = NDData(data_cube, meta=meta)
+
+    grid = GriddedPSFModel(data, fill_value=fill_value)
+
+    return grid