EPSF Builder making non-continuous PSF

[moira-andrews]

We’re having an issue with EPSFBuilder where our resultant epsf appears to not be resampling properly when we have an oversampling greater than 1. We are inputing 32 stars, and have an oversampling of 5. When looking at the 3d and top down view, it is clear that pixel to pixel is not continuous. We’ve tried increasing the number of stars to as many as 123, but the ending epsf has the same issue.

Here is the code making the epsf and the 3d and 2d plots:

from astropy.io import fits
import matplotlib.pyplot as plt
import numpy as np
from astropy.nddata import NDData
from astropy.wcs import WCS
from astropy import units as u
from astropy.coordinates import SkyCoord
from photutils.psf import EPSFBuilder
from photutils.psf import extract_stars
from astropy.table import Table

obj_filepath ='cpt1m013-fa14-20210412-0125-e91.fits'
hdu = fits.open(obj_filepath)
data = hdu[0].data
cat_data = hdu[1].data
cat_header = hdu[1].header
header = hdu[0].header
hdu.close()

wcs = WCS(header)

cat_ra = cat_data['ra']
cat_dec = cat_data['dec']
cat_fwhm = cat_data['fwhm']
cat_peak = cat_data['peak']

stars_wcscoords = SkyCoord(ra=cat_ra*u.degree, dec=cat_dec*u.degree, frame='icrs')
stars_pixcoords_x = wcs.world_to_pixel(stars_wcscoords)[0]
stars_pixcoords_y = wcs.world_to_pixel(stars_wcscoords)[1]

size = 25
hsize = (size - 1) / 2
mask = ((stars_pixcoords_x > hsize) & (stars_pixcoords_x < (data.shape[1] -1 - hsize)) &
        (stars_pixcoords_y > hsize) & (stars_pixcoords_y < (data.shape[0] -1 - hsize)) & 
        (cat_fwhm > 4) & (cat_fwhm < 5) & (cat_peak > 250) & (cat_peak < 100000))

stars_tbl = Table()
stars_tbl['x'] = stars_pixcoords_x[mask]
stars_tbl['y'] = stars_pixcoords_y[mask]

nddata = NDData(data=data)  

stars = extract_stars(nddata, stars_tbl, size=25)  

epsf_builder = EPSFBuilder(oversampling=5, maxiters=1, progress_bar=False)
epsf, fitted_stars = epsf_builder.build_epsf(stars) 

x = np.arange(epsf.shape[1])
y = np.arange(epsf.shape[0])
x, y = np.meshgrid(x, y)

fig = plt.figure(1)
ax = fig.add_subplot(111, projection='3d')
ax.plot_surface(x, y, epsf.data, cmap='viridis')
ax.set_xlabel('X')
ax.set_ylabel('Y')
ax.set_title('32 stars')
plt.show()

fig = plt.figure(2)
plt.imshow(epsf.data, origin='lower')
plt.title('32 stars')
plt.colorbar()
plt.show()

cpt1m013-fa14-20210412-0125-e91.fits.fz.zip

[larrybradley]

@moira-andrews With an oversampling factor of 5, you need a lot more stars to get a good EPSF. For an oversampling=5, you would need to ensure that at least 1 star's true center is located in each of the 5x5 subpixel locations. That means for a "perfect" sampling a star centers, you need at least 25 stars to fill in the oversampled grid. Otherwise, you will get weird artifacts like holes in the resulting EPSF. But even with 25 perfectly subsampled stars, each star would only contribute once to each pixel, so you will end up with a very noisy EPSF. Ideally, you should have many hundreds of stars that can fill the oversampled grid. I would recommend lowing the oversampling factor. Your stars are very well resolved, so you can probably even use oversampling=1.