Matching PSF sims to in-flight JWST data

import stpsf
import numpy as np
import matplotlib
import matplotlib.pyplot as plt
import astropy, astropy.io.fits as fits

Often one wants to generate PSFs matched to some particular science dataset or file. The convenience function stpsf.setup_sim_to_match_file helps with this, using the file’s FITS header to set up a simulated instrument matched to the appropriate instrument setup and date of observation.

Let’s call that function, providing a filename of a data file already downloaded in the current directory:

filename = 'jw02739010001_02103_00001_nrcalong_cal.fits'

inst = stpsf.setup_sim_to_match_file(filename)

Setting up sim to match jw02739010001_02103_00001_nrcalong_cal.fits
iterating query, tdelta=3.0

MAST OPD query around UTC: 2023-04-05T03:33:47.591
                        MJD: 60039.14846748843

OPD immediately preceding the given datetime:
        URI:     mast:JWST/product/R2023040302-NRCA3_FP1-1.fits
        Date (MJD):      60037.0387
        Delta time:      -2.1098 days

OPD immediately following the given datetime:
        URI:     mast:JWST/product/R2023040504-NRCA3_FP1-1.fits
        Date (MJD):      60039.4669
        Delta time:      0.3184 days
User requested choosing OPD time closest in time to 2023-04-05T03:33:47.591, which is R2023040504-NRCA3_FP1-1.fits, delta time 0.318 days
Importing and format-converting OPD from /Users/mperrin/software/stpsf-data/MAST_JWST_WSS_OPDs/R2023040504-NRCA3_FP1-1.fits
Backing out SI WFE and OTE field dependence at the WF sensing field point

Configured simulation instrument for:
    Instrument: NIRCam
    Filter: F335M
    Detector: NRCA5
    Apername: NRCA5_FULL
    Det. Pos.: (1024, 1024)
    Image plane mask: None
    Pupil plane mask: None

This function will:

• Create a stpsf instrument object for the relevant instrument

• Configure it to have the correct filter, detector, and other relevant instrument parameters for that science data file (e.g. coronagraph masks and so on).

• Load the measured telescope mirror alignment data from the closest-in-time wavefront sensing visit to that science data.

A reminder about output data products

Recall that the output data files from a simulated PSF are multi-extension FITS, containing both oversampled and detector-sampled outputs, and with and without added distortion effects. In general there will be 4 output FITS HDUs:

psf = inst.calc_psf()
psf.info()

Filename: (No file associated with this HDUList)
No.    Name      Ver    Type      Cards   Dimensions   Format
  0  OVERSAMP      1 PrimaryHDU     116   (316, 316)   float64
  1  DET_SAMP      1 ImageHDU       118   (79, 79)   float64
  2  OVERDIST      1 ImageHDU       170   (316, 316)   float64
  3  DET_DIST      1 ImageHDU       171   (79, 79)   float64

fig, axes = plt.subplots(figsize=(12,3), ncols=4)
for i in range(len(psf)):
    stpsf.display_psf(psf, ext=i, ax=axes[i], title=f'Ext {i}: {psf[i].header["EXTNAME"]}',
                       imagecrop=2, colorbar=False)

The difference between the oversampled and detector-sampled output products is readily apparent. The distortion effects are generally more subtle:

• In this example case, note the slightly blurred softer look of the DET_DIST output compared to DET_SAMP, or of OVERDIST compared to OVERSAMP. This aspect of the simulation is a model for charge transfer physics and inter-pixel capacitance within the detector which result in crosstalk between adjacent pixels.

• Also included as part of the distortion is a model for optical geometric distortions (including for instance slight differences between X and Y pixel scales, small rotations and skews of the detector pixel axes, the very-slightly-different position angles for each NIRCam detector, etc.). This attempts to forward-model the distortions which the “drizzle” pipeline algorithm corrects for, using the same astrometric calibration information for the instruments recorded in the science instrument aperture file.

In general it is this last extension 3, DET_DIST, that is expected to have the greatest fidelity to measured PSFs observed in detector coordinates, such as in ``rate.fits`` or ``cal.fits`` pipeline data products. Users are encouraged to use this part of the PSF sim output in most cases. Note, JWST has many instrument modes, and the precise level of fidelity achieved compared to observed PSFs varies depending on instrument and mode.

No STPSF output product precisely corresponds to PSFs as seen in pipeline-interpolated i2d.fits files; doing so would be very difficult because the “drizzle” algorithm is not strictly preserving of PSF morphology, and furthermore drizzled output products typically are generated from multiple dithered images and thus PSFs in such files are derived as averages of PSFs sampled at different positions per each dither. We recommend working with rate or cal files, not i2d files, for tasks requiring the most precise PSF calibration.

Understanding and matching positions on the detector

In general, to most precisely match a given observed PSF, you will need to set parameters for both ``detector_location`` and ``source_offset_{x/y}``. Set detector_location to the integer pixel coordinates for the center of some cutout subregion containing a PSF (similar to setting the center for e.g. astropy.nddata.Cutout2D to cut out that piece of the science data). Then set source_offset_x and source_offset_y to adjust the subpixel position of the target with respect to the center of that cutout, as described further below.

Adjusting the position of the simulated subregion

You can think of the simulation output as a model for a PSF on a particular subregion or cutout of the detector pixels in a given image. STPSF includes models and calibration data that allow simulating the small but nonzero variations in PSFs depending on the location of that cutout within the field of view.

By default, the PSF calculation is configured for the center of the selected detector, which will often give a reasonable match to PSFs across the detector. This may or may not be sufficient for any given use case, depending on the particular science needs.

To specify a particular detector region of interest, simply adjust the detector_position attribute. This attribute sets which pixel in the image (i.e. which pixel in your science data FITS file) corresponds to the center pixel location of the subregion of pixels output by stpsf.

(Note, for science images taken using a detector subarray instead of the full detector, after using setup_sim_to_match_data the detector_pixel coordinates will refer to the pixel coordinates in that subarray only, not the full detector. In other words, simply set detector_pixel to the desired (X,Y) coordinates within the specified science image file, and it will work correctly for either full-frame or subarray science data).

Here’s an example of setting the position to match a PSF cutout from a science image:

# Specify a location and box size. There's a source visible at these coordinates in the file we specified above.
psf_center = (1996, 1413)  # note this is in X, Y order
boxsize = 50

# Load that science data. Cut out the surface brightness and uncertainty for the desired location
obs_im = fits.open(filename)

obs_psf = astropy.nddata.Cutout2D(obs_im['SCI'].data, position=psf_center, size=boxsize).data
obs_psf_err = astropy.nddata.Cutout2D(obs_im['ERR'].data, position=psf_center, size=boxsize).data

obs_psf -= np.nanmedian(obs_psf)    # perform a simple background subtraction

# Configure the instrument model to that same center location

inst.detector_position = psf_center

# Compute PSF model
sim_psf = inst.calc_psf(fov_pixels=boxsize)

# Define a convenience function for comparisons. We'll use this several times in this notebook.

def plot_data_sim_comparison(obs_psf, obs_psf_err, sim_psf):
    """Display observed data, model, and difference; and compute a goodness of fit metric
    """
    fig, axes = plt.subplots(figsize=(13,3), ncols=3)

    vmax = np.nanmax(obs_psf)
    cmap = matplotlib.cm.gist_heat
    cmap.set_bad(cmap(0))
    axes[0].imshow(obs_psf,
                   norm = matplotlib.colors.LogNorm(vmax/1e4, vmax),
                   cmap=cmap, origin='lower')
    axes[0].set_title("Observed PSF from science data")

    stpsf.display_psf(sim_psf, ext='DET_DIST', vmax=0.1, vmin=1e-5, ax=axes[1], )

    scalefactor = np.nansum(obs_psf)
    difference = obs_psf - sim_psf["DET_DIST"].data *scalefactor
    axes[2].imshow(difference,
                   norm = matplotlib.colors.LogNorm(vmax/1e5, vmax),
                   cmap=cmap, origin='lower')
    # compute a simple goodness-of-fit metric
    chisqr = np.nansum(difference**2/obs_psf_err**2)/np.isfinite(difference).sum()
    axes[2].set_title("Difference")
    axes[2].text(3, 3, f"$\\chi^2$ = {chisqr:.2f}", color='white')
    for ax in [axes[0], axes[2]]:
        ax.set_xlabel("Pixels")

plot_data_sim_comparison(obs_psf, obs_psf_err, sim_psf)

Adjusting the position of the simulated point source

By default the PSF will be exactly centered in that output grid of pixels. Think of this as simulating a star that is precisely centered in that subarray of pixels on the detector.

To adjust the position of that point source, including for subpixel positions, set options['source_offset_{x/y}']. These allow you to specify an angular offset in arcseconds.

inst.options['source_offset_x'] = 0.010   # in units of arcseconds. So, this is about 1/5 of a NRC LW pixel
inst.options['source_offset_y'] = 0.005

sim_psf_offset = inst.calc_psf(fov_pixels=boxsize)

# Plot the data, showing that the PSF fit is slightly improved.
plot_data_sim_comparison(obs_psf, obs_psf_err, sim_psf_offset)

Other options to potentially improve fidelity

If further improvements in PSF fidelity are needed, there are additional knobs that can be turned. For example, you may also find it useful to adjust the .options['charge_diffusion_sigma'] paraneter as well. This provides a (highly simplified!) way to adjust the model to account for charge diffusion from the so-called “brighter fatter effect”. See the docs page about JWST detector effect models for more, or you may choose to consult the JWST help desk for advice on applications to particular science datasets. This level of precise tuning may not be needed for many science cases.

Here we show an example of adjusting the charge diffusion model, which is one way to slightly tweak the PSF FWHM to better match a given dataset.

inst.options['charge_diffusion_sigma'] = 0.028

sim_psf_offset_v2 = inst.calc_psf(fov_pixels=boxsize)

# Plot the data, showing that the PSF fit is again slightly improved.
plot_data_sim_comparison(obs_psf, obs_psf_err, sim_psf_offset_v2)

Checking telescope stability around a particular science observation.

JWST PSFs are in general quite stable over time, though not always. It may be useful to check whether there were significant wavefront variations around the time of a given observation, which can be done with the show_wfs_around_obs function. See the wavefront measurements docs page for further details.

In this case, we can verify that the telescope wavefront was stable for the two-day period including that science observation. The delta wavefront error is barely detectable (lower right plot below). The typical wavefront sensing measurement precision is about 6 nm rms, so the measured 6.9 nm rms shown here is just barely more than that sensing precision. This is an example of the excellent stability typical of most time periods.

stpsf.trending.show_wfs_around_obs(filename)

File jw02739010001_02103_00001_nrcalong_cal.fits observed at 2023-04-05T03:33:47.591
Retrieving WFS before that obs... WFS at 2023-04-03T00:55:44.076
Retrieving WFS after that obs... WFS at 2023-04-05T11:12:18.992