Drizzled PSFs from HST/WFC3 Imaging

This example notebook shows how spike interacts with HST/WFC3 imaging to generate a drizzled PSF. This notebook uses observations from this UVIS program and this IR program. The easiest way to find these data in the MAST archive is by searching for observation ID idxf58020 (UVIS) and icor18020 (IR). spike uses calibrated, but not-yet co-added images (i.e., “Level 2” data products), so be sure that your download includes ‘flt’/’flc’ files. If you download all of the files associated with this program, you will have two working directories, one for the WFC3/UVIS observations which will include:

• idxf58a3q_flc.fits

• idxf58a6q_flc.fits

• idxf58aaq_flc.fits

• idxf58ztq_flc.fits

• idxf58zvq_flc.fits

• idxf58zyq_flc.fits

and one for the WFC3/IR observations which will include:

• icor18thq_flt.fits

• icor18tkq_flt.fits

• icor18toq_flt.fits

• icor18tsq_flt.fits

These data are also available here.

The principle of this example notebook is the same regardless of dataset, though, so it can be used as a direct guide for use with other data, too.

NOTE: Environment variables are not always read in properly within Jupyter notebooks/JupyterLab – see here for some additional details of how to access them.

from spike.psf import hst # import the relevant top-level module from spike

outputpath = 'psfs' #default is /psfs, defined from your working directory

Note that if you choose to create drizzled PSFs for the same data and object location using multiple PSF generation methods, you will want to define a new ‘savedir’ each time to avoid conflicts between files with the same name.

WFC3/UVIS

We will begin by working with WFC3/UVIS data.

datapath = '/path/to/uvis/data'

We’ll also define an object – in this case, a set of coordinates – for which the PSF will be generated.

obj = '10:00:31.9753 +02:11:52.964' #can also be a resolvable name or coordinates in decimal degrees

Now, we can call spike.psf.hst – as an initial example we’ll use the most minimal inputs. By default, in addition to saving drizzled PSFs and any intermediate products in the directory given by ‘savedir’ (here the argument will be the outputpath variable we defined), this function will return a dictionary that stores model PSFs indexed by filter and object.

If returnpsf = 'full' (default), the arrays in the dictionary will be the full drizzled image field of view. If returnpsf = 'crop', the arrays in the dictionary will be the cutout region around the object location with the size of the cutout determined by the ‘cutout_fov’ argument; the cropped region around the drizzled PSF will also be saved as a FITS file, including its WCS information (can be toggled off with savecutout = False).

spike assumes that the input files have not yet been “tweaked” and handles that step for you. If your input images have been tweaked or if you would like to skip this step for other reasons, specify pretweaked = True in calling spike.psf.hst.

Since the default PSF generation method for spike.psf.hst is TinyTim, which is not recommended for use with WFC3, we will select method = 'STDPSF'.

The outputs are very long due to all of the print statements from the tweak, PSF generation, and drizzling steps. Outputs have been cleared for cells where spike.psf.hst was run.

dpsf = hst(img_dir = datapath, obj = obj, img_type = 'flc', inst = 'WFC3', camera = 'UVIS',
           savedir = outputpath, pretweaked = True, method = 'STDPSF')

In this case dpsf will be very simple, since we are only looking at one set of coordinates and one filter:

dpsf

{'10:00:31.9753 +02:11:52.964': {'F275W': array([[nan, nan,  0., ..., nan, nan, nan],
         [nan, nan,  0., ..., nan, nan, nan],
         [ 0.,  0.,  0., ..., nan, nan, nan],
         ...,
         [nan, nan, nan, ..., nan, nan, nan],
         [nan, nan, nan, ..., nan, nan, nan],
         [nan, nan, nan, ..., nan, nan, nan]],
        shape=(2311, 4132), dtype='>f4')}}

spike can automatically plot model PSFs as they are generated, and generates PSFs regardless of whether they are plotted (or returned). To show the output drizzled PSFs here, we will plot them directly. Therefore, we will import matplotlib below. Since returnpsf = 'full' by default, the array stored in dpsf shows the drizzled PSF in the context of the entire co-added image footprint. This is equivalent to reading in the image stored in ‘150d07m59.6295s+2d11m52.964s_F275W_psf_drc.fits’.

import matplotlib.pyplot as plt
import numpy as np

dat = dpsf['10:00:31.9753 +02:11:52.964']['F275W']

fig = plt.figure()
plt.imshow(dat, vmin = np.nanpercentile(dat[dat != 0], 20), vmax = np.nanpercentile(dat[dat != 0], 75),
            origin = 'lower', cmap = 'Greys')
plt.xlabel('x (pix)')
plt.ylabel('y (pix)')

Text(0, 0.5, 'y (pix)')

We can index the array to just show the PSF below:

dat = dpsf['10:00:31.9753 +02:11:52.964']['F275W'][1182:1322, 2103:2243]

fig = plt.figure()
plt.imshow(dat, vmin = np.nanpercentile(dat[dat != 0], 20), vmax = np.nanpercentile(dat[dat != 0], 75),
            origin = 'lower', cmap = 'Greys', extent = [2103, 2243, 1182, 1322])
plt.xlabel('x (pix)')
plt.ylabel('y (pix)')

Text(0, 0.5, 'y (pix)')

Additional user choices

As noted above, by default, PSFs are generated for spike.psf.hst by TinyTim using a blackbody model of a G5V star at 6000 K, but users have complete control over the specifics of the PSF generation. To change the parameters within a method (e.g., updating model arguments, setting different detection thresholds for the empirical PSFs, …), spike.psf.hst directly takes keyword arguments for the specified ‘method’. Details of allowed arguments are available in the spike.psfgen documentation.

spike can also be used flexibly with a variety of different PSF generation methods – to change which one is used, simply use the ‘method’ argument. Built-in options are:

• ‘TinyTim’ (not recommended for WFC3)

• ‘STDPSF’

• ‘epsf’

• ‘PSFEx’

If method = 'USER', then the keyword argument ‘usermethod’ must also be specified, pointing to either a function or a directory that houses model PSFs that follow the [imgprefix]_[coords]_[band]_psf.fits, e.g., imgprefix_23.31+30.12_F814W_psf.fits or imgprefix_195.78-46.52_F555W_psf.fits naming scheme.

spike can also be used to generate PSFs for multiple objects/coordinates simultaneously. We’ll continue using STDPSFs below and have spike output cutouts around the drizzled PSF.

obj = ['10:00:31.9753 +02:11:52.964', '10:00:31.2334 +02:13:13.167']

dpsf = hst(img_dir = datapath, obj = obj, img_type = 'flc', inst = 'WFC3', camera = 'UVIS',
           savedir = outputpath, pretweaked = True, method = 'STDPSF', returnpsf = 'crop')

dpsf

{'10:00:31.9753 +02:11:52.964': {'F275W': array([[0., 0., 0., ..., 0., 0., 0.],
         [0., 0., 0., ..., 0., 0., 0.],
         [0., 0., 0., ..., 0., 0., 0.],
         ...,
         [0., 0., 0., ..., 0., 0., 0.],
         [0., 0., 0., ..., 0., 0., 0.],
         [0., 0., 0., ..., 0., 0., 0.]], shape=(151, 151), dtype='>f4')},
 '10:00:31.2334 +02:13:13.167': {'F275W': array([[0., 0., 0., ..., 0., 0., 0.],
         [0., 0., 0., ..., 0., 0., 0.],
         [0., 0., 0., ..., 0., 0., 0.],
         ...,
         [0., 0., 0., ..., 0., 0., 0.],
         [0., 0., 0., ..., 0., 0., 0.],
         [0., 0., 0., ..., 0., 0., 0.]], shape=(151, 151), dtype='>f4')}}

Since we’ve selected returnpsf = 'crop', the output array is already cropped to the immediate region around the PSF. Selecting both returnpsf = 'crop' and savecutout = True (default), will save a FITS file that stores the cropped drizzled PSF – the name is the same as the for the full saved image with ‘_crop’ appended.

dat = dpsf['10:00:31.9753 +02:11:52.964']['F275W']

fig, ax = plt.subplots(1, 2, gridspec_kw = {'wspace':0.4})
ax[0].imshow(dat, vmin = np.nanpercentile(dat[dat != 0], 20), vmax = np.nanpercentile(dat[dat != 0], 80),
            origin = 'lower', cmap = 'Greys')
ax[0].set_xlabel('x (pix)')
ax[0].set_ylabel('y (pix)')

dat = dpsf['10:00:31.2334 +02:13:13.167']['F275W']

ax[1].imshow(dat, vmin = np.nanpercentile(dat[dat != 0], 20), vmax = np.nanpercentile(dat[dat != 0], 80),
            origin = 'lower', cmap = 'Greys',)
ax[1].set_xlabel('x (pix)')
ax[1].set_ylabel('y (pix)')

Text(0, 0.5, 'y (pix)')

WFC3/IR

spike also works with WFC3/IR data.

datapath = '/path/to/ir/data'

obj = '10:00:36.9898 +02:09:14.229'

dpsf = hst(img_dir = datapath, obj = obj, img_type = 'flt', inst = 'WFC3', camera = 'IR',
           savedir = outputpath, pretweaked = True, method = 'STDPSF')

dpsf

{'10:00:36.9898 +02:09:14.229': {'F160W': array([[nan, nan, nan, ..., nan, nan, nan],
         [nan, nan, nan, ..., nan, nan, nan],
         [nan, nan, nan, ..., nan, nan, nan],
         ...,
         [ 0.,  0.,  0., ..., nan, nan, nan],
         [ 0.,  0.,  0., ..., nan, nan, nan],
         [nan, nan, nan, ..., nan, nan, nan]],
        shape=(963, 1089), dtype='>f4')}}

dat = dpsf['10:00:36.9898 +02:09:14.229']['F160W']

fig = plt.figure()
plt.imshow(dat, vmin = np.nanpercentile(dat[dat != 0], 20), vmax = np.nanpercentile(dat[dat != 0], 95),
            origin = 'lower', cmap = 'Greys')
plt.xlabel('x (pix)')
plt.ylabel('y (pix)')

Text(0, 0.5, 'y (pix)')

dat = dpsf['10:00:36.9898 +02:09:14.229']['F160W'][399:539, 387:527]

fig = plt.figure()
plt.imshow(dat, vmin = np.nanpercentile(dat[dat != 0], 20), vmax = np.nanpercentile(dat[dat != 0], 95),
            origin = 'lower', cmap = 'Greys', extent = [387, 527, 399, 539])
plt.xlabel('x (pix)')
plt.ylabel('y (pix)')

Text(0, 0.5, 'y (pix)')