negative pixel grid in ePSF

[jryon]

I am attempting to make an ePSF from a set of 18 stars in the field of an ACS/SBC F125LP drizzled image. I'm following the example here. The resulting ePSF has a strange grid pattern of negative pixels, which appears to me as some kind of aliasing. I've explored a range of values for oversampling, maxiters, and norm_radius, but the grid pattern remains (though the appearance changes). I looked through related issues here, but I couldn't find a solution for this problem, nor am I confident it is not a bug. There may be poor subpixel sampling of the peak locations of the PSF stars, since there are relatively few bright, isolated stars in the image.

My code and the resulting ePSF are below, and I can share the drizzled image if needed. I'm using version 0.7.2 of photutils.

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from astropy.io import fits
from astropy.visualization import simple_norm
from astropy.table import Table
from astropy.nddata import NDData
from photutils.psf import extract_stars
from photutils import EPSFBuilder

if __name__ == '__main__':

    # Make table out of PSF star positions
    tab = Table.read('psf_stars.coo', format='ascii')

    tab['col1'].name = 'x'
    tab['col2'].name = 'y'

    # Read in SBC image
    hdu = fits.open('../images/native/F125LP_NGC5253_ACS_SBC_drz.fits')
    data = hdu[1].data

    nddata = NDData(data=data)

    # cut out selected stars
    stars = extract_stars(nddata, tab, size=25)

    # Create EPSFBuilder instance and run
    epsf_builder = EPSFBuilder(oversampling=4, norm_radius=10.5)

    epsf, fitted_stars = epsf_builder(stars)

    # Plotting
    norm = simple_norm(epsf.data, 'log', percent=99.)

    plt.ion()
    plt.figure()

    plt.imshow(epsf.data, norm=norm, origin='lower', cmap='viridis')
    plt.colorbar()

[eteq]

Hmm, curious. Is there any way you can try this with one of the non-drizzled images that fed into the drizzled image? I realize it's the drizzled one you want, but this might help distinguish between a problem that the drizzling is creating somehow vs the epsf fitter itself.

[larrybradley]

@jryon Typically you need many more than 18 stars to create a good ePSF. At an absolute minimum, for 4x oversampling you need 16 stars, but those stars would need to be centered such that they uniformly sample the 4x4 subpixel grid (not likely with random stars) -- and even in that case the ePSF would be very noisy because you'd only get one sample at each ePSF pixel. The result you are seeing is very likely due to holes in the 4x oversampled ePSF grid due to the lack of stars. Is it possible for you to add more input stars?

[jryon]

@larrybradley I increased the number of stars to 31, but that is nearly the max available in the field that aren't contaminated by other stars, and some of them are fairly dim. I also tried 2x oversampling (then I'd only need 8 stars, if uniformly sampled?), but I'm still seeing a grid pattern, though it is less extreme with certain choices of parameter values. If I use quadratic smoothing, the grid pattern is less noticeable around the middle of the PSF, but then a row or column of negative pixels appears at the edge of the PSF frame.

I'm working on testing with an FLT image instead of the DRZ I was using. Thanks for your help!

[jryon]

Using the FLT image jb7h06eaq_flt.fits and the list of stars here (x and y coordinates), I get the following PSFs. I think I still see a grid pattern in the PSF with quartic smoothing, but I don't see one in the quadratic smoothing case. Both have strong edge effects, however.

[jryon]

I obtained a very deep drizzled SBC image of a star cluster and attempted to make a PSF from ~50 high S/N, isolated stars in the image. I used the same parameters as for the FLT PSFs above: 2x oversampling, norm_radius=10, and maxiters=10. The PSF with quadratic smoothing again looks pretty good, though a bit of a grid pattern appears in the top and bottom of the frame. The PSF with quartic smoothing again has a very strong grid pattern. Both have edge effects.

At this point, it seems like using drizzled frames is not the issue, and neither is the number of stars provided to EPSFBuilder, so I'd hazard a guess that something is wrong with the code.

[larrybradley]

The edge effects you are seeing is because of the fixed cutout size of the input stars. The input star cutouts (using extract_stars) are centered on the initial guess for the center. During the ePSF building process, the star centers are continually updated. When the star centers are improved, the cutouts effectively need to be shifted in x/y to stack them in the oversampled grid. That can cause a lack of data on the edges, which gets worse if there aren't enough stars (or if the initial guesses aren't very accurate). You can simply crop the edges. Most of the PSF flux is in the central few pixels, so that's the most critical region for PSF fitting.

In your last post, the "quadratic" smoothing looks pretty good to me, considering you only have 50 stars. I think Jay Anderson would tell you that you need at least several hundred well-dithered stars to build a good ePSF. I'm a bit surprised the "quartic" looks significantly worse, but IIRC that particular smoothing kernel was generated by Jay for 4x oversampling. I've had conversations with Jay about the smoothing kernels and unfortunately, the best kernel depends on the shape of the PSF. It's a bit of an art.

Have you talked to Jay about this? He may have generated a library of ePSFs for the ACS/SBC.

[jryon]

Ok, the edge effects make sense now, thanks for the explanation. The PSF I'm attempting to make is for GALFIT fitting, not PSF photometry. In this case, I think the wings are more important.

In many cases, it is very difficult or impossible to find hundreds of well-dithered, isolated, bright stars in an image. If an ePSF is not considered good unless that criterion is satisfied, that may limit EPSFBuilder's usefulness to the average user.

Is there anything like IRAF or DAOPHOT's psf in photutils? Or plans to implement it? In Anderson & King 2000, there's some discussion that DAOPHOT PSFs are poor for astrometry, but fine for PSF photometry. And from my experience with IRAF's psf, a reasonable PSF model can be built from an image with 25-30 stars. While I appreciate that ePSFs are more accurate in general, a reliable, less accurate PSF may work fine for some users, and wouldn't have the artifacts like the ones I've seen here.

I'm on the ACS instrument team, and we have ePSF models for WFC and HRC only. The lack of an SBC ePSF library is possibly because most fields have a small number of UV-bright stars.

[keflavich]

It would be great to include some advisory documentation in the docs about the EPSF builder. I just encountered the same issue with JWST data with 2800 stars:

There's still some gridding in the quadratic version, but it's further out and not so bad.

[keflavich]

You can simply crop the edges.

What's the appropriate way to do this? The naive approach

epsf.data = epsf.data[2:-2, 2:-2]

doesn't work because .data is an unwriteable attribute. One could do

epsf._data = epsf.data[2:-2, 2:-2]

but that's using a private variable and not encouraged. Would it be best to do

epsf = EPSFModel(data=epsf.data[2:-2, 2:-2])

or does that lose important information about, e.g., normalization?

[andywinhold]

I'm having a similar issue of a grid in the final ePSF. This thread has been helpful so I wanted to share where I'm at. If I've created an act against all that is sacred let me know.

tl;dr: making the kernel not have values below zero seems to help.

I went into my local copy of epsf.py and added some code to accept/create a 'custom' kernel. After trying different oversampling values and array sizes (quadratic and quartic are 5x5) I still had the grid issue. It seems like it tracked with the corners of the quartic or quadratic kernels having negative values in the corners so I just adjusted my custom kernel to keep everything zero or greater.

kernel[kernel < 0] = 0

Here is a before and after example of the same ePSF. Oversampling: 4x, polynomial2D degree: 2, kernel array size: 7x7. One image, 72 stars found (some may be hot pixels but removing them doesn't change the grid artefact).

Before:

After:

        elif self.smoothing_kernel == 'custom':
            from astropy.modeling import models, fitting
            edge = 7
            array = np.zeros((edge,edge))
            mid = int(np.floor(edge/2))
            array[mid,mid] = 1
            gy, gx = np.mgrid[:edge,:edge]
            pinit = models.Polynomial2D(degree=2)
            fitp = fitting.LevMarLSQFitter()
            
            p = fitp(pinit,gx,gy,array)
            kernel = np.asarray(p(gx,gy))
            kernel[kernel < 0] = 0
        else:
            raise TypeError('Unsupported kernel.')
            
        return kernel, convolve(epsf_data, kernel)

[ghost]

I am trying to perform Quasar PSF subtraction after getting my epsf from this code. But I am not getting a good subtracted image which I presume is due to some issues in my psf.
This is what I am getting as a result of generating the psf for my image:

Though the psf looks good the only issue is that my quasar subtracted image isn't what it is supposed to look like. It looks something like:

The black region clearly shows that it hasn't been properly subtracted. Does anyone have any idea as to what could be the issue?
Any help would be appreciated

[Astro-Sean]

HI @AuroraBorealis30 - I am running into a similar issue with the PSF subtractions - did you find a solution to this?