PSF photometry just uses circular PSFs

[StuartLittlefair]

As pointed out by @coclar in #80, the PSF photometry option simply fits circularly symmetric PSFs to the stars in the images.

The parameters for these PSFs are simply copied across from those found in moveApers without an extra fitting step. This causes bad results in the case of elliptical PSFs.

In #105 , I introduced new code that uses PSF functions that are - in theory - capable of having different widths in the x and y directions. To implement elliptical PSFs the following steps need doing:

• [ ] the PSF models need a theta argument to allow for rotation
• [ ] we need an extra fitting step before PSF photometry to set the ellipticity

Note the second part cannot be done by simply setting the parameters on the PSF model to variable instead of fixed. This is because if you do that, the photutils PSF fitting code will fit each star independently and obviously the shape parameters cannot be constrained by faint stars.

Instead we will need a separate step using the reference stars to set the PSF model parameters before the PSF photometry step.

[StuartLittlefair]

Some progress on this, based on some experiments I have done with PSF photometry on HiPERCAM data which was marginally successful (see attached script).

The issue with the PSF models as described above is still current and not solved, and if we want to use asymmetric Moffat or Gaussian profiles we still need an additional fitting step. However, there is an alternative that might work, making an effective PSF (EPSF) from the data itself.

Based on this I one way to make PSF photometry useful within the pipeline is this:

1) Make a shiftadd script, which is a version of combine that aligns images based on star positions (I have a rough draft already). shiftadd would be used to make a "master" frame (combine would do for this in the short term).

2) Select bright, isolated reference stars in each image to make a an Effective PSF. This is probably best done selecting those stars with setaper and then running a new command make_epsf that will allow users to tweak the parameters of EPSF creation and view the result. The epsf star list and successful settings are saved to a file.

3) Next another new script make_source_list is run, using the EPSF settings and the master frame. This is based on the attached script and uses Iterative PSF Photometry to create a master source list, including detecting blended stars.

4) Finally, using the EPSF info and the master source list of stars from (3), we can run the reduce where for each frame a new EPSF is created using the EPSF star list, and forced photometry is carried out using the master source list (with a position tweak calculated from the EPSF stars).

[StuartLittlefair]

import sys
from astropy.io import fits
from astropy.table import QTable
from photutils.background import LocalBackground, MMMBackground
from astropy.nddata import NDData
from photutils.psf import IterativePSFPhotometry, PSFPhotometry
from photutils.detection import DAOStarFinder
from astropy.visualization import simple_norm
from astropy.stats import sigma_clipped_stats
from photutils.psf import EPSFBuilder
from photutils.psf import extract_stars
from photutils.utils import CutoutImage
import matplotlib.pyplot as plt
from photutils.psf import SourceGrouper
from photutils.centroids import centroid_1dg, centroid_com

stamp_shape = (21, 21)
fit_shape = (5, 5)
fwhm_estimate = 12

data = fits.getdata(sys.argv[1])
# subtract BKG
mean_val, median_val, std_val = sigma_clipped_stats(data, sigma=2.0)
bkgsub_data = data - median_val

# find stars for EPSF building
finder = DAOStarFinder(50.0, fwhm_estimate)
peaks_tbl = finder(bkgsub_data)
stars_tbl = QTable()
stars_tbl["x"] = peaks_tbl["xcentroid"]
stars_tbl["y"] = peaks_tbl["ycentroid"]

# CHOOSE GOOD STARS
use_idx = []
for i in range(len(stars_tbl)):
    x, y = stars_tbl[i]["x"], stars_tbl[i]["y"]
    try:
        co = CutoutImage(bkgsub_data, (y, x), shape=stamp_shape, mode="strict")
    except:
        continue
    norm = simple_norm(co.data, "sqrt", percent=99.0)
    plt.imshow(co.data, norm=norm, origin="lower", cmap="viridis")
    use = input("Use? ")
    if int(use):
        use_idx.append(i)
    plt.close()
psfstar_tbl = QTable(rows=[stars_tbl[i] for i in use_idx], names=("x", "y"))
stars = extract_stars(NDData(bkgsub_data), psfstar_tbl, size=stamp_shape)

# BUILD EPSF
epsf_builder = EPSFBuilder(
    oversampling=1,
    maxiters=10,
    smoothing_kernel=None,
    norm_radius=fwhm_estimate * 1.5,
    recentering_boxsize=(fwhm_estimate + 3, fwhm_estimate + 3),
    progress_bar=False,
    recentering_func=centroid_com,
)
epsf, fitted_stars = epsf_builder(stars)
norm = simple_norm(epsf.data, "linear", percent=98.0)
plt.imshow(epsf.data, norm=norm, origin="lower", cmap="viridis")
plt.colorbar()
plt.show()

continue_ = input("Continue? ")
if continue_.lower() != "y":
    sys.exit()

########## PSF Photometry

# crop around the transient
data = data[255:350, 525:655]
# subtract BKG
mean_val, median_val, std_val = sigma_clipped_stats(data, sigma=2.0)
bkgsub_data = data - median_val
# create NDData object
nddata = NDData(data=bkgsub_data)

# First pass for initial params
psfphot = PSFPhotometry(
    epsf, fit_shape, finder=finder, aperture_radius=fwhm_estimate * 1.2
)
phot = psfphot(nddata)
init_params = QTable()
init_params["x"] = phot["x_fit"]
init_params["y"] = phot["y_fit"]

# Iterative PSF to find all stars
grouper = SourceGrouper(min_separation=25)
bkgstat = MMMBackground()
localbkg_estimator = LocalBackground(25, 35, bkgstat)
psfphot2 = IterativePSFPhotometry(
    epsf,
    fit_shape,
    finder=finder,
    maxiters=5,
    localbkg_estimator=localbkg_estimator,
    aperture_radius=fwhm_estimate * 1.5,
    # grouper=grouper,
)
phot2 = psfphot2(nddata, init_params=init_params)

print(phot2[["x_fit", "y_fit", "flux_fit", "flux_err", "iter_detected", "group_id"]])

fig, ax = plt.subplots(nrows=1, ncols=2, figsize=(15, 5), sharex=True, sharey=True)
resid = psfphot2.make_residual_image(nddata, stamp_shape)
norm = simple_norm(nddata.data, "linear", percent=98)
ax[0].imshow(nddata.data, origin="lower", norm=norm)
ax[0].plot(phot2["x_fit"], phot2["y_fit"], "r.")
im = ax[1].imshow(resid.data, origin="lower")
ax[0].set_title("Data")
ax[1].set_title("Residual Image")
plt.tight_layout()
plt.show()