JWST Using Wavefronts Measured On Orbit

JWST Using Wavefronts Measured On Orbit#

STPSF, formerly knows as WebbPSF, includes code for using the results of in-flight wavefront sensing measurements, by retrieving Optical Path Difference (OPD) files. These can be used to create simulated PSFs appropriate for a given instrument, detector, and time.

This notebook serves as a reference for how to use STPSF to search for and retrieve an OPD file, and explains the contents of this file. It also shows how to load multiple OPD files and plot changes in the wavefront properties over time. It also contains an example case of using an OPD file to create and subtract a simulated PSF from an in-flight image.

Details on the use of STPSF’s in-flight OPDs are given here: https://stpsf.readthedocs.io/en/latest/jwst_measured_opds.html

Author: Marcio Melendez Hernandez
Last Updated: 13 Feb 2025

Index#

  • Imports and Setup

  • Finding the measured wavefront near a given date

  • OPD File description

  • Delta Wavefront Error Around Time of Observation

  • Setup a PSF simulation using a particular observation

  • Trending Wavefront Changes over Time

  • Wavefront time series and histogram plots

    • Table of all OPDs

    • Trending Plots

    • Plot over a time range

    • Load and plot a Single OPD

  • JWST Simulated PSF Subtraction from in-flight data

    • Extract a subarray around source of interest, and measure source location

    • Find the source

    • Use previous function to set up our simulation

    • Use photutils to create our model PSF

    • Subtract PSF and create residual image

    • PSF properties and differences

  • Improved IFU Sims

For other STPSF examples: https://stpsf.readthedocs.io/en/latest/more_examples.html

Imports and Setup#

import matplotlib.pylab as plt
import poppy
import numpy as np
import os
import datetime
import tarfile
from urllib.parse import urlparse
import requests

import stpsf

import astropy
from astropy.nddata import Cutout2D
from astropy.visualization import simple_norm
from astropy.io import fits
from astropy.stats import SigmaClip
from astropy.visualization.mpl_normalize import ImageNormalize
from astropy.visualization import LogStretch
from astropy.modeling import models, fitting

import photutils
from photutils.background import Background2D, MedianBackground
from photutils.detection import IRAFStarFinder
from photutils.psf import PSFPhotometry

# Files containing such information as the JWST pupil shape, instrument throughputs, and aperture positions are distributed separately from STPSF.
# To run STPSF, you must download these files and tell STPSF where to find them using the STPSF_PATH environment variable.

# Set environmental variables
os.environ["STPSF_PATH"] = "./data/stpsf-data"

# STPSF Data
stpsf_url = 'https://stsci.box.com/shared/static/kqfolg2bfzqc4mjkgmujo06d3iaymahv.gz'
stpsf_file = './stpsf-data-LATEST.tar.gz'
stpsf_folder = "./data"


def download_file(url, dest_path, timeout=60):
    parsed_url = urlparse(url)
    if parsed_url.scheme not in ["http", "https"]:
        raise ValueError(f"Unsupported URL scheme: {parsed_url.scheme}")

    response = requests.get(url, stream=True, timeout=timeout)
    response.raise_for_status()
    with open(dest_path, "wb") as f:
        for chunk in response.iter_content(chunk_size=8192):
            f.write(chunk)


# Gather stpsf files
stpsfExist = os.path.exists(stpsf_folder)
if not stpsfExist:
    os.makedirs(stpsf_folder)
    download_file(stpsf_url, stpsf_file)
    gzf = tarfile.open(stpsf_file)
    gzf.extractall(stpsf_folder, filter='data')

Finding the measured wavefront near a given date#

nrc = stpsf.NIRCam()
nrc.filter = 'F200W'

nrc.detector = 'NRCB2'
nrc.detector_position = (1024, 1024)

The interpixel capacitance effect is now included for oversampled outputs as well as for the detector-sampled outputs (in the geometric distortion extension). Remember that, in general, the last (“DET_DIST”) FITS extension of the output PSF FITS file are the output data product that most represents the PSF as actually observed on a detector. In any case, there are way to disable any detector effects or to adjust the empirical approximation of charge difussion. For example: nrc.options['charge_diffusion_sigma'] = 0

nrc.options['add_ipc'] = False

output_path = os.getcwd()
nrc.load_wss_opd_by_date('2022-07-01T00:00:00', plot=True, output_path=output_path)
# choice : string . Default 'closest'
# Method to choose which OPD file to use, e.g. 'before', 'after'

  • Upper Left: This is the measured OPD as sensed in NIRCam at “field point 1” which is in the upper left corner of NRCA3, relatively close to the center of the NIRCam module. This observatory total OPD measurement includes both the telescope and NIRCam contributions to the WFE.

  • Upper Middle: This is the wavefront map for the NIRCam portion of the WFE at that field point. This is known from ground calibration test data, not measured in flight.

  • Upper Right: That NIRCam WFE contribution is subtracted from the total observatory WFE to yield this estimate of the OTE-only portion of the WFE.

  • Lower Left and Middle: These are models for the field dependence of the OTE OPD between the sensing field point in NRCA3 and the requested field ooint, in this case in NRCB2. This field dependence arises mostly from the figure of the tertiary mirror. These are used to transform the estimated OTE OPD from one field position to another.

  • Lower Right: This is the resulting estimate for the OTE OPD at the requested field point, in this case in NRCB2.

# What's inside an OPD?
opd_fn = 'R2022070403-NRCA3_FP1-0.fits'
opd = fits.open(os.path.join(output_path, opd_fn))

opd.info()

norm = ImageNormalize(stretch=LogStretch(), vmin=1e-10, vmax=1e-6)
plt.figure(figsize=[10, 18])
for i in range(1, len(opd)):
    plt.subplot(5, 3, i)
    if i == 2:
        plt.imshow(opd[i].data, norm=norm, origin='lower')
    else:
        plt.imshow(opd[i].data, origin='lower')

OPD File description#

  • RESULT_PHASE Image Extension. The FITS file contains a ‘RESULT’ image extension, which is the average of the optical path difference results from all images analyzed together in the phase retrieval

  • RESULT_PSF Image Extension The FITS file contains a ‘RESULT_PSF’ image extension, which is the PSF computed from the resultant phase by the WAS

  • EXPECTED Image Extension The FITS file contains an ‘EXPECTED’ image extension, which is the expected optical path difference if the WAS-recommended correction is applied

  • PUPIL_MASK Image Extension The FITS file contains a ‘PUPIL_MASK’ image extension, which is the pupil mask used to compute the PSF from the resultant phase

The FITS file contains five image extensions for each analyzed input image. In the case of the sensing program, there are 2 images +/-8WL

  • RAW_PSF Image Extensions For each input image, the Raw PSF extension will contain the raw extracted subimage taken from the calibrated science data

  • CALC_PSF Image Extensions For each input image, the Calculated PSF extension will contain an image which represents the estimated PSF as calculated by the phase retrieval process

  • CALC_AMP Image Extension For each input image, the Calculated Amplitude extension will contain an image which represents the estimated pupil amplitude as calculated by the phase retrieval process

  • HO_PHASE Image Extension For each input image, the High-Order Phase extension will contain an image which represents the retrieved phase information minus the Low-Order (controllable) phase as calculated by the phase retrieval process.

  • LO_PHASE Image Extension For each input image, the Low-Order Phase extension will contain an image which represents the retrieved controllable phase information as calculated by the phase retrieval process

# Let's simulate Webb's PSF from the OPD above
# The PSF below is calculated using the actual 
# as-measured-at-L2 state of the telescope WFE near 
# the requested date, in this case 2022 July 1.
psf = nrc.calc_psf(fov_pixels=101)

plt.figure(figsize=(16, 9))
stpsf.display_psf(psf, ext=1)

Delta Wavefront Error Around Time of Observation#

stpsf.trending.delta_wfe_around_time('2024-05-11 01:50:25.231')

Setup a PSF simulation using a particular observation#

def mast_retrieve_files(filename, output_path=None, verbose=False, redownload=False, token=None):
    """Download files from MAST.
    If file is already present locally, the download is skipped and the cached file is used.
    """
    import os
    from requests.exceptions import HTTPError
    from astroquery.mast import Mast
    if token:
        Mast.login(token=token)

    if output_path is None:
        output_path = '.'
    else:
        output_path = output_path

    output_filename = os.path.join(output_path, filename)

    if not os.path.exists(output_path):
        os.mkdir(output_path)

    if os.path.exists(output_filename) and not redownload:
        if verbose:
            print(f"Found file previously downloaded: {filename}")
        return output_filename

    data_uri = f'mast:JWST/product/{filename}'

    # Download the file
    url_path = Mast._portal_api_connection.MAST_DOWNLOAD_URL
    try:
        Mast._download_file(url_path + "?uri=" + data_uri, output_filename)
    except HTTPError as err:
        print(err)
    return output_filename

file = 'jw04500-o056_t012_nircam_f212n-wlm8-nrca3_wfscmb-04.fits'
mast_retrieve_files(file)

inst = stpsf.setup_sim_to_match_file(file)
position = (512, 1024) # source position
size_pixels = 128 # size in pixels
inst.detector_position = position

single_stpsf_nrc = inst.calc_psf(fov_pixels=size_pixels, display=True)

Wavefront time series and histogram plots#

Table of all OPDs#

Retrieve a table of all available OPDs and plot the measurements over time

opdtable = stpsf.mast_wss.retrieve_mast_opd_table()

opdtable = stpsf.mast_wss.deduplicate_opd_table(opdtable)

opdtable

Download all the OPDs from the opdtable object. Note that some functions need to have all the OPDs available. 

stpsf.mast_wss.download_all_opds(opdtable)

Plot over a time range#

start_date = datetime.datetime.fromisoformat('2022-11-22')
end_date = datetime.datetime.fromisoformat('2022-12-01')
stpsf.trending.wavefront_time_series_plot(opdtable, start_date=start_date,
                                          end_date=end_date, ymax=85, ymin=60,
                                          label_visits=True, label_events=True)

# We can also plot all measured wavefront drifts over specified time periods
start_time = astropy.time.Time('2024-01-01T00:00:00')
end_time = astropy.time.Time.now()

stpsf.trending.wavefront_drift_plots(opdtable, start_time=start_time, end_time=end_time, n_per_row=10)

# We can plot histograms of wavefront error levels over time
stpsf.trending.wfe_histogram_plot(opdtable, thresh=70, ote_only=True)

# Or select a particular time period
start_date = astropy.time.Time('2024-04-01')
end_date = astropy.time.Time.now()
stpsf.trending.wfe_histogram_plot(opdtable, start_date=start_date, thresh=70, ote_only=True)

Load and plot a Single OPD#

Inspect the OPD table above for an OPD around your observation

nrc = stpsf.NIRCam()
nrc.filter = 'F212N'
opd_fn = 'R2022120204-NRCA3_FP1-1.fits'
nrc.load_wss_opd(opd_fn, output_path=output_path)
fov_pixels = 511
psf = nrc.calc_psf(oversample=4, fov_pixels=fov_pixels)

psf.info()

# display the PSF and plot the encircled energy
plt.figure(figsize=(12, 6))
plt.subplot(1, 2, 1)
stpsf.display_psf(psf, colorbar_orientation='horizontal')
axis2 = plt.subplot(1, 2, 2)
stpsf.display_ee(psf, ext=1, ax=axis2)

JWST Simulated PSF Subtraction from in-flight data#

https://stpsf.readthedocs.io/en/latest/jwst_psf_subtraction.html

file = 'jw02725-o484_t097_nircam_clear-f356w-nrcalong_wfscmb-04.fits'
mast_retrieve_files(file)

# Open and inspect observation
psf = fits.open(file)
data = psf['SCI'].data
mask = (psf['DQ'].data % 2).astype('bool')
sigma_clip = SigmaClip(sigma=3.0)
bkg_estimator = MedianBackground()
# norm = simple_norm(data1, 'log',min_cut = min_cut, max_cut = max_cut)
bkg = Background2D(data, (50, 50), filter_size=(3, 3),
                   sigma_clip=sigma_clip, bkg_estimator=bkg_estimator, mask=mask)

plt.figure(figsize=(6, 6))
norm = simple_norm(data - bkg.background, 'asinh', min_percent=20, max_percent=98.99)
plt.imshow(data - bkg.background, norm=norm, origin='lower', cmap='viridis')

Extract a subarray around source of interest, and measure source location#

position = (550, 1950) # approximate location. This will be improved below.
size_pixels = 200

data_source = Cutout2D(data - bkg.background, position, size_pixels).data
plt.figure(figsize=(6, 6))
norm = simple_norm(data_source, 'asinh', min_percent=20, max_percent=99.6)
plt.imshow(data_source, norm=norm, origin='lower', cmap='viridis')
plt.title("Cutout around a galaxy merger and a star")

Find the source#

starfind = IRAFStarFinder(fwhm=3.0, threshold=100. * bkg.background_median)
sources = starfind(data_source)
plt.figure(figsize=(6, 6))
plt.imshow(data_source, norm=norm, origin='lower', cmap='viridis')
plt.title("Cutout around a galaxy merger and a star")
plt.scatter(sources['xcentroid'], sources['ycentroid'], color='red', marker='x')

# Let's re-center the image 
positions_original = Cutout2D(data - bkg.background, position, size_pixels).to_original_position((sources['xcentroid'], sources['ycentroid']))
position = (positions_original[0].value[0], positions_original[1].value[0])
size = 200
data_source = Cutout2D(data - bkg.background, position, size_pixels).data

plt.figure(figsize=(6, 6))
norm = simple_norm(data_source, 'asinh', min_percent=20, max_percent=99.6)
plt.imshow(data_source, norm=norm, origin='lower', cmap='viridis')
plt.title("Cutout centered on the star")
print(position)

Use previous function to set up our simulation#

inst = stpsf.setup_sim_to_match_file(file)
inst.detector_position = position

single_stpsf_nrc = inst.calc_psf(fov_pixels=size_pixels)

single_stpsf_nrc[3].header

Use photutils to create our model PSF#

We’ll use photutils to fit and subtract the PSF from the data. First, we take the output simulated PSF and convert it into a photutils fittable model:

stpsf_model = photutils.psf.FittableImageModel(single_stpsf_nrc['DET_DIST'].data, normalize=True, oversampling=1)

fit_shape = (5, 5)
psfphot = PSFPhotometry(stpsf_model, fit_shape,
                        finder=starfind, aperture_radius=4)

# fit model PSF to background subtracted data
phot = psfphot(data_source)

Check the location of the source to be subtracted#

print(phot) # single source

Subtract PSF and create residual image#

residual = psfphot.make_residual_image(data_source, psf_shape=(size_pixels, size_pixels))

plt.figure(figsize=(14, 16))

plt.subplot(1, 2, 1)
plt.imshow(data_source, norm=norm, origin='lower', cmap='viridis')
plt.title('Original - zoom')

plt.subplot(1, 2, 2)
plt.imshow(residual, norm=norm, origin='lower', cmap='viridis')
plt.title('Clean - PSFPhotometry zoom')

PSF properties and differences#

pixelscale = nrc.pixelscale 
ee_pixel_radius = 2.5
ee_arcsec_radius = ee_pixel_radius * pixelscale
ee_psf = poppy.measure_ee(single_stpsf_nrc, ext=3, normalize='total')
ee_val = ee_psf(ee_arcsec_radius)
print("ee ({}px, {:.3f}arsec)  = {:.4f}".format(ee_pixel_radius, ee_arcsec_radius, ee_val.item(0)))

def measure_fwhm(array):
    """Fit a Gaussian2D model to a PSF and return the fitted PSF
    the FWHM is x and y can be found with fitted_psf.x_fwhm, fitted_psf.y_fwhm

    Parameters
    ----------
    array : numpy.ndarray
        Array containing PSF

    Returns
    -------
    x_fwhm : float
        FWHM in x direction in units of pixels

    y_fwhm : float
        FWHM in y direction in units of pixels
        
    x_mean : float
        x centroid position in units of pixels
    
    y_mean : float
        y centroid position in units of pixels
    
    """
    yp, xp = array.shape
    y, x = np.mgrid[:yp, :xp]
    p_init = models.Gaussian2D(amplitude=array.max(), x_mean=xp * 0.5, y_mean=yp * 0.5)
    fit_p = fitting.LevMarLSQFitter()
    fitted_psf = fit_p(p_init, x, y, array)
    return fitted_psf

fitted_psf = measure_fwhm(single_stpsf_nrc["DET_SAMP"].data)
print("FWHM X-direction: {:.3f}, FWHM y-direction: {:.2f}".format(fitted_psf.x_fwhm, fitted_psf.y_fwhm))

# Note that the OPD is loaded in STPSF instrument object
# norm = ImageNormalize(stretch=LinearStretch(), vmin = 1e-9 , vmax = 1e-7)
plt.figure(figsize=[10, 10])
plt.imshow(nrc.pupilopd[0].data, origin='lower')

# check OPD header information
nrc.pupilopd[0].header

# Let's calculate the difference between two different times
# We can use the same instrument setup from before
opd_fn = 'R2022120404-NRCA3_FP1-1.fits'
inst.load_wss_opd(opd_fn)
psf2 = inst.calc_psf(fov_pixels=size_pixels)

pixelscale = nrc.pixelscale
ee_pixel_radius = 2.5
ee_arcsec_radius = ee_pixel_radius * pixelscale
ee_psf = poppy.measure_ee(psf2, ext=1, normalize='total')
ee_val = ee_psf(ee_arcsec_radius)
print("ee ({}px, {:.3f}arsec)  = {:.4f}".format(ee_pixel_radius, ee_arcsec_radius, ee_val.item(0)))

stpsf.display_psf_difference(single_stpsf_nrc, psf2, imagecrop=2, title='Difference between two OPDs', cmap='gist_heat')

def miri_psfs_for_ee():
    miri = stpsf.MIRI()

    # opd_fn = 'R2022120404-NRCA3_FP1-1.fits'
    # miri.load_wss_opd(opd_fn, output_path = output_path)
    miri.load_wss_opd_by_date('2022-07-12T00:00:00', plot=True, output_path=output_path)

    for wave in [5.0, 7.5, 10, 14]:
        fov = 18

        outname = "PSF_MIRI_%.1fum_wfed.fits" % (wave)
        psf = miri.calc_psf(outname, monochromatic=wave * 1e-6,
                            oversample=4, fov_arcsec=fov, display=True)
    return psf


def plot_ee_curves():
    plt.clf()
    for iw, wave in enumerate([5.0, 7.5, 10, 14]):

        ees60 = []
        ees51 = []
        ax = plt.subplot(2, 2, iw+1)
       
        name = "PSF_MIRI_%.1fum_wfed.fits" % (wave)
        stpsf.display_ee(name, ax=ax, mark_levels=False)

        eefn = stpsf.measure_ee(name)
        ees60.append(eefn(0.60))
        ees51.append(eefn(0.51))

        ax.text(1, 0.6, 'Mean EE inside 0.60": %.3f' % np.asarray(ees60).mean())
        ax.text(1, 0.5, 'Mean EE inside 0.51": %.3f' % np.asarray(ees51).mean())

        ax.set_title(f"Wavelength = {wave:.1f} $\\mu$m")

        ax.axvline(0.6, ls=":", color='k')
        ax.axvline(0.51, ls=":", color='k')

    plt.tight_layout()   

miri_psfs_for_ee()

plot_ee_curves()

nsp = stpsf.NIRSpec()
# or you can specify a full path name.
# please make an output PSF with its center
# aligned to the center of a single pixel
nsp.options['parity'] = 'odd'

opd_fn = 'R2022120404-NRCA3_FP1-1.fits'
nsp.load_wss_opd(opd_fn, output_path=output_path)
    
waves = np.linspace(0.8, 5, 50) * 1e-6 # iterate over wavelengths in meters

for iw, wavelength in enumerate(waves):
    psffile = 'psf_NIRSPec_mono_%.1fum_opd1.fits' % (wavelength * 1e6)
    psf = nsp.calc_psf(fov_arcsec=3, oversample=4, 
                       monochromatic=wavelength, display=False,
                       outfile=psffile)
    ax = plt.subplot(8, 8, iw + 1)
    stpsf.display_psf(psffile, ext='DET_SAMP', colorbar=False, imagecrop=8)
    ax.set_title('')
    ax.xaxis.set_visible(False)
    ax.yaxis.set_visible(False)
    ax.text(-3.5, 0, '{0:.1f}'.format(wavelength * 1e6))

plt.figure(figsize=(8, 12))
nsp.image_mask = 'MSA all open'
nsp.display()
msapsf = nsp.calc_psf(monochromatic=2e-6, oversample=8)
stpsf.display_psf(msapsf, ext='DET_SAMP')

miri_psfs_for_ee()

plot_ee_curves()

Improved IFU sims#

miri = stpsf.MIRI()
miri.mode = 'IFU'
miri.band = '2A'
waves = miri.get_IFU_wavelengths()
cube = miri.calc_datacube_fast(waves)

cube.info()

nrs = stpsf.NIRSpec()
nrs.mode = 'IFU'
nrs.disperser = 'PRISM'
nrs.filter = 'CLEAR'
waves = nrs.get_IFU_wavelengths()
cube = nrs.calc_datacube_fast(waves)

cube.info()
Contents