Drizzled PSFs from HST/ACS Imaging
This example notebook shows how spike interacts with HST/ACS imaging to generate a drizzled PSF. This notebook uses observations from this program. The easiest way to find these data in the MAST archive is by searching for observation ID j8pu42010. spike uses calibrated, but not-yet co-added images (i.e., “Level 2” data products), so be sure that your download includes ‘flc’ files. If you download all of
the files associated with this program, your working directory will include:
j8pu42ecq_flc.fits
j8pu42egq_flc.fits
j8pu42esq_flc.fits
j8pu42evq_flc.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.
Initial setup and visualization
[ ]:
from spike.psf import hst # import the relevant top-level module from spike
datapath = '/path/to/acs/data'
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.
We’ll also define an object – in this case, a set of coordinates – for which the PSF will be generated.
[ ]:
obj = '10:00:33.0178 +02:09:52.304' #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.
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 = 'ACS', camera = 'WFC',
savedir = outputpath, pretweaked = True)
In this case dpsf will be very simple, since we are only looking at one set of coordinates and one filter:
[4]:
dpsf
[4]:
{'10:00:33.0178 +02:09:52.304': {'F475W': array([[ 0., 0., 0., ..., nan, nan, nan],
[ 0., 0., 0., ..., nan, nan, nan],
[ 0., 0., 0., ..., nan, nan, nan],
...,
[nan, nan, nan, ..., 0., 0., 0.],
[nan, nan, nan, ..., 0., 0., 0.],
[nan, nan, nan, ..., 0., 0., 0.]],
shape=(2359, 4162), 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
‘150d08m15.267s+2d09m52.304s_F475W_psf_drc.fits’.
[5]:
import matplotlib.pyplot as plt
import numpy as np
dat = dpsf['10:00:33.0178 +02:09:52.304']['F475W']
fig = plt.figure()
plt.imshow(dat, vmin = np.nanpercentile(dat[dat != 0], 20), vmax = np.nanpercentile(dat[dat != 0], 97),
origin = 'lower', cmap = 'Greys')
plt.xlabel('x (pix)')
plt.ylabel('y (pix)')
[5]:
Text(0, 0.5, 'y (pix)')
We can index the array to just show the PSF below:
[6]:
dat = dpsf['10:00:33.0178 +02:09:52.304']['F475W'][325:465, 940:1080]
fig = plt.figure()
plt.imshow(dat, vmin = np.nanpercentile(dat[dat != 0], 20), vmax = np.nanpercentile(dat[dat != 0], 97),
origin = 'lower', cmap = 'Greys', extent = [940, 1080, 325, 465])
plt.xlabel('x (pix)')
plt.ylabel('y (pix)')
[6]:
Text(0, 0.5, 'y (pix)')
Specifying user preferences
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. As an example, we will change the TinyTim model inputs to use an O6 star at 45000 K:
[ ]:
dpsf = hst(img_dir = datapath, obj = obj, img_type = 'flc', inst = 'ACS', camera = 'WFC',
savedir = outputpath, pretweaked = True, listchoice = 'O6', temp = 45000)
As before, we’ll plot the outputs directly:
[8]:
dat = dpsf['10:00:33.0178 +02:09:52.304']['F475W']
datcrop = dpsf['10:00:33.0178 +02:09:52.304']['F475W'][325:465, 940:1080]
fig, ax = plt.subplots(1, 2, gridspec_kw = {'width_ratios':[1.75, 1], 'wspace':0.3})
ax[0].imshow(dat, vmin = np.nanpercentile(dat[dat != 0], 20), vmax = np.nanpercentile(dat[dat != 0], 97),
origin = 'lower', cmap = 'Greys')
ax[0].set_xlabel('x (pix)')
ax[0].set_ylabel('y (pix)')
ax[1].imshow(datcrop, origin = 'lower', cmap = 'Greys', extent = [940, 1080, 325, 465],
vmin = np.nanpercentile(datcrop[datcrop != 0], 20),
vmax = np.nanpercentile(datcrop[datcrop != 0], 97))
ax[1].set_xlabel('x (pix)')
ax[1].set_ylabel('y (pix)')
[8]:
Text(0, 0.5, 'y (pix)')
Generating multiple PSFs
spike can also be used to generate PSFs for multiple objects/coordinates simultaneously.
[ ]:
obj = ['10:00:33.0178 +02:09:52.304', '10:00:28.1798 +02:08:40.507']
dpsf = hst(img_dir = datapath, obj = obj, img_type = 'flc', inst = 'ACS', camera = 'WFC',
savedir = outputpath, pretweaked = True)
[10]:
dpsf
[10]:
{'10:00:33.0178 +02:09:52.304': {'F475W': array([[ 0., 0., 0., ..., nan, nan, nan],
[ 0., 0., 0., ..., nan, nan, nan],
[ 0., 0., 0., ..., nan, nan, nan],
...,
[nan, nan, nan, ..., 0., 0., 0.],
[nan, nan, nan, ..., 0., 0., 0.],
[nan, nan, nan, ..., 0., 0., 0.]],
shape=(2359, 4162), dtype='>f4')},
'10:00:28.1798 +02:08:40.507': {'F475W': array([[ 0., 0., 0., ..., nan, nan, nan],
[ 0., 0., 0., ..., nan, nan, nan],
[ 0., 0., 0., ..., nan, nan, nan],
...,
[nan, nan, nan, ..., 0., 0., 0.],
[nan, nan, nan, ..., 0., 0., 0.],
[nan, nan, nan, ..., 0., 0., 0.]],
shape=(2359, 4162), dtype='>f4')}}
[11]:
dat = dpsf['10:00:33.0178 +02:09:52.304']['F475W']
datcrop = dpsf['10:00:33.0178 +02:09:52.304']['F475W'][325:465, 940:1080]
fig, ax = plt.subplots(2, 2, gridspec_kw = {'width_ratios':[1.75, 1], 'wspace':0.3})
ax[0, 0].imshow(dat, vmin = np.nanpercentile(dat[dat != 0], 20), vmax = np.nanpercentile(dat[dat != 0], 97),
origin = 'lower', cmap = 'Greys')
ax[0, 0].set_xlabel('x (pix)')
ax[0, 0].set_ylabel('y (pix)')
ax[0, 1].imshow(datcrop, origin = 'lower', cmap = 'Greys', extent = [940, 1080, 325, 465],
vmin = np.nanpercentile(datcrop[datcrop != 0], 20),
vmax = np.nanpercentile(datcrop[datcrop != 0], 97))
ax[0, 1].set_xlabel('x (pix)')
ax[0, 1].set_ylabel('y (pix)')
dat = dpsf['10:00:28.1798 +02:08:40.507']['F475W']
datcrop = dpsf['10:00:28.1798 +02:08:40.507']['F475W'][1504:1644, 2605:2745]
ax[1, 0].imshow(dat, vmin = np.nanpercentile(dat[dat != 0], 20), vmax = np.nanpercentile(dat[dat != 0], 97),
origin = 'lower', cmap = 'Greys')
ax[1, 0].set_xlabel('x (pix)')
ax[1, 0].set_ylabel('y (pix)')
ax[1, 1].imshow(datcrop, origin = 'lower', cmap = 'Greys', extent = [2605, 2745, 1504, 1644],
vmin = np.nanpercentile(datcrop[datcrop != 0], 20),
vmax = np.nanpercentile(datcrop[datcrop != 0], 97))
ax[1, 1].set_xlabel('x (pix)')
ax[1, 1].set_ylabel('y (pix)')
[11]:
Text(0, 0.5, 'y (pix)')
PSF generation methods, automatic cropping
One advantage of spike is that it can 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’
‘TinyTim_Gillis’
‘STDPSF’
‘epsf’
‘PSFEx’
‘ACS_ePSF’ – added in v1.2.4
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.
We’ll use STDPSFs below and have spike output cutouts around the drizzled PSF.
[ ]:
obj = '10:00:33.0178 +02:09:52.304'
dpsf = hst(img_dir = datapath, obj = obj, img_type = 'flc', inst = 'ACS', camera = 'WFC',
savedir = outputpath, pretweaked = True, method = 'STDPSF', returnpsf = 'crop')
[14]:
dpsf
[14]:
{'10:00:33.0178 +02:09:52.304': {'F475W': 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.
[15]:
dat = dpsf['10:00:33.0178 +02:09:52.304']['F475W']
fig = plt.figure()
plt.imshow(dat, vmin = np.nanpercentile(dat[dat != 0], 20), vmax = np.nanpercentile(dat[dat != 0], 97),
origin = 'lower', cmap = 'Greys')
plt.xlabel('x (pix)')
plt.ylabel('y (pix)')
[15]:
Text(0, 0.5, 'y (pix)')
