GriddedPSFModel evaluated PSFs with the same input flux have very different output fluxes

[york-stsci]

This has been tested with a GriddedPSFModel created out of a set of webbpsf Roman WFI apertures, using the sample code below. Note that instead of using webbpsf directly to create the PSF grid, instead the grid is created image-by-image and put into an NDData array after creation. This is because webbpsf currently does not preserve parity for even oversamples, and the simulator that I'm using these models for requires odd parity to build its EPSF models.

[york-stsci]

import numpy as np

from astropy.io import fits
from astropy.nddata import NDData
from photutils.psf import GriddedPSFModel
from webbpsf import wfirst as roman

def create_gridded_psf(fov_pixels, oversample, psf_grid):
    meta = {'oversampling': oversample, 'grid_xypos': []}
    pixel_scale = 0.11
    num_psf = psf_grid*psf_grid
    fov = fov_pixels * oversample
    if fov%2 == 0:
        fov += 1
    data = np.zeros((num_psf, fov, fov), dtype=np.float64)
    if num_psf == 1:
        psf_xgrid = [2048.]
        psf_ygrid = [2048.]
    else:
        pix_low, pix_high = 4., 4092.
        psf_xgrid = np.linspace(4., 4092., psf_grid)
        psf_ygrid = np.linspace(4., 4092., psf_grid)
    i = 0
    for yp in psf_ygrid:
        for xp in psf_xgrid:
            meta['grid_xypos'].append((xp, yp))
            ins = roman.WFI()
            ins.options['parity'] = 'odd'
            ins.filter = 'F087'
            ins.detector = 'SCA01'
            ins.detector_position = (xp, yp)
            ins.pixelscale = pixel_scale/oversample
            psf_fits = ins.calc_psf(fov_pixels=fov, oversample=1)
            psf = psf_fits["DET_SAMP"].data
            print("PSF at ({},{}) has total flux {}".format(xp, yp, np.sum(psf)))
            print("PSF Shape is {}".format(psf.shape))
            data[i][:] = psf * oversample**2
            i += 1
    return GriddedPSFModel(NDData(data, meta=meta))

grid = create_gridded_psf(65, 10, 5)
half_grid = 32
x_loc = [1234.56, 2198.75, 3114.03, 155.603, 1610.203, 2048.0, 2048.0]
y_loc = [1234.56, 729.14, 2541.89, 1338.08, 1338.04, 1338.04, 1338.08]
fluxes = [1000.0, 1000.0, 1000.0, 39.482, 39.482, 100.0, 100.0]

for (x,y,flux) in zip(x_loc, y_loc, fluxes):
    ix, iy = int(round(x)), int(round(y))
    lx = ix - half_grid
    hx = ix + half_grid + 1
    ly = iy - half_grid
    hy = iy + half_grid + 1
    xg, yg = np.mgrid[lx:hx, ly:hy]
    psf_img = grid.evaluate(xg, yg, flux=flux, x_0=x, y_0=y)
    print("Image at ({},{}) with flux {} has sum {}".format(x, y, flux, np.sum(psf_img)))

[larrybradley]

@york-stsci Many thanks for the code example. I'll try to take a look soon.

[larrybradley]

The issue here is that GriddedPSFImage (and FittableImageModel, EPSFMode) sample from the interpolation function instead of integrating over bins of size oversampling * oversampling. For large oversampling factors, this can lead to significant differences, especially when there are large gradients in the model (e.g., near the peak) within blocks of size oversampling*oversampling. The flux differences you are seeing are almost all due to the PSF peak. In cases where your desired x_0/y_0 PSF position is near the center of a pixel (e.g., integer pixel values), you will get more output flux than expected. This is because the brightest pixel will be sampled, but the mean of the 10x10 surrounding region will be smaller. On the other hand, when x_0/y_0 are closer to half-pixel values, you will get less flux. I have ideas how to address this issue, but it will require a completely different and more complicated approach to implementing an image-based PRF model class. In the meantime, if you lower your oversampling factor (to something like 4, typically used for HST), you should get much better results.