HSTPSF

Bayesian PSF posterior sampling for HST data. Currently only works for HST FLT images for the F814W band. I am open to expanding this if you want to get involved!

How to use:

Step 1: Install

Create a fresh virtual environment and then install the requirements for the package found in requirements.txt. There's a good chance you can just run this:

python3 -m venv newvenv
source newvenv/bin/activate
pip install -r requirements.txt

Step 2: Collect data

You need to collect the relevant data from your HST FLT F814W image. This includes a 32x32 cutout of: data ("SCI" extension), mask ("DQ" extension), and err ("ERR" extension). You will also need an estimate of the sky level, which you can just use the median of all pixel not flagged by the DQ image. You will also need the exposure time, which you can get from the header. For most FLT images it would look something like this:

from astropy.io import fits
from astropy.wcs import WCS
from astropy.coordinates import SkyCoord
from astropy.nddata import Cutout2D
import numpy as np

ver = 2 # layer with your data
_ra = 0.0 # location of your point source
_dec = 0.0

with fits.open("j9op01l8q_flt.fits") as hdu:
    exposure_time = hdu[0].header.get("EXPTIME")
    data = hdu["SCI", ver].data
    mask = np.array(hdu["DQ", ver].data != 0, dtype=bool)
    err = hdu["ERR", ver].data
    sky = np.nanmedian(data[~mask])
    wcs = WCS(hdu["SCI", ver].header, fobj=hdu)

sky_coord = SkyCoord(ra=_ra, dec=_dec, unit="deg", frame="icrs")
data = Cutout2D(data, sky_coord, size=32, wcs=wcs, mode="partial", fill_value=np.nan).data
mask = Cutout2D(mask, sky_coord, size=32, wcs=wcs, mode="partial", fill_value=np.nan).data
err = Cutout2D(err, sky_coord, size=32, wcs=wcs, mode="partial", fill_value=np.nan).data

np.savez("my_data.npz", data=data, mask=mask, err=err, sky=sky, exposure_time=exposure_time)

The above code is available in the load_from_fits.py file included here. We also include a demo file with some example data in it at demo_data.npz.

Step 3: Run the sampler

You will need to run the sample_cutout_posterior function from the sample_psf_posterior file. A demo will run if you call the file with python from the command line:

python sample_psf_posterior.py

Ths will demonstrate what the results look like. To run with your own data, lets presume you have a npz file from the previous step, here is how you would set the analysis going:

import numpy as np
from sample_psf_posterior import sample_psf_posterior

D = np.load("my_data.npz")

sample_psf_posterior(
    data = D["data"],
    mask = D["mask"],
    err = D["err"],
    sky = D["sky"],
    exposure_time = D["exposure_time"],
)

Example outputs

First, caskade makes a graph to represent the simulator. This is relatively simple since it is just a pixelated PSF that will be scaled and interpolated into the data space.

Next, an initial fit for the nuisance parameters will be done with a gaussian, and then the mean PSF prior image. This is what the fit looks like for the mean prior image:

Then, the full posterior sampling will begin with an iterative process to anneal and settle the nuisance parameters. Here is one of the intermediate steps:

Finally, a run with more SDE steps and corrector steps will be run to allow for the best posterior sampling possible:

This is just one example (and not even an especially good example), the output .h5 file contains all the fitted parameters, including 32 posterior PSF samples. Here's what the collection looks like for the demo:

Which you can extract from the .h5 file like this:

import h5py
import numpy as np
import matplotlib.pyplot as plt

with h5py.File("psf_posterior_samples_final.h5", "r") as f:
    D = f["cutout/psf/image/value"][()]

fig, axarr = plt.subplots(4, 8, figsize=(32,16))

for i, ax in enumerate(axarr.flatten()):
    ax.imshow(np.log10(D[i]), cmap="Greys_r", vmin=-5, vmax = 1)
    ax.axis("off")
plt.savefig("psf_posterior_samples.png", dpi = 300, bbox_inches = "tight")
plt.close()