This example was written to help users better understand the subtleties in improving image sampling for dithered data. One of the powers of the drizzling algorithm is that, given properly dithered images, it can restore much of the information lost due to undersampled images (Fruchter and Hook, 2002).
This work is based on ISR ACS 2015-01, which contains a more detailed discussion than presented here.
In practice, this requires the use of Astrodrizzle task within the Drizzlepac package. This example will cover the following topics:
1. Description of the Data
2. Run astrodrizzle several times using different settings for the 'final_pixfrac' and 'final_scale' parameters
3. Compare and evaluate results of using different 'final_pixfrac' and 'final_scale' valuesimport glob
import shutil
import os
from astropy.io import fits
from astroquery.mast import Observations
from astropy.visualization import LogStretch,ImageNormalize,LinearStretch
from astropy import wcs
import drizzlepac
import matplotlib.pyplot as plt
import numpy as np
%matplotlib inline
Four images of the spiral galaxy NGC 3370 (Program 11570) were acquired using the WFC3/IR F160W imaging filter. The WFC3-IR-DITHERBOX-MIN dither pattern designed to provide optimal sampling of the PSF was used.
Data are downloaded using the astroquery API to access the MAST archive. The astroquery.mast documentation has more examples for how to find and download data from MAST. For this example, the caibrated FLT exposures from visit 19 will be downloaded.
# Querying MAST and downloading files.
obs_ids = ['ib1f19010']
obs_table = Observations.query_criteria(obs_id=obs_ids, obstype='all')
download_tab = Observations.download_products(obs_table['obsid'], mrp_only=False, productSubGroupDescription=['FLT'])
# Cleaning up directories after dowloading from MAST
for file in download_tab['Local Path']:
os.rename(file, os.path.basename(file))
shutil.rmtree('mastDownload')
For the drizzle algorithm to work optimally, observations need to be dithered so that the PSF is optimally sampled. The code below creates plots to show the dithering of each of the images in the association. The plot on the left shows how the images were dithered on the sky (the POSTARG). The plot on the right shows how the offsets translate to pixel phase (e.g. subpixel position).
flt_files = glob.glob('*flt.fits')
plate_scale = fits.getval(flt_files[0], 'idcscale',ext=1)
print('Plate scale: {:>6.4f}'.format(plate_scale))
postarg1 = np.empty(len(flt_files), dtype=np.float)
postarg2 = np.empty(len(flt_files), dtype=np.float)
x_phase = np.empty(len(flt_files), dtype=np.float)
y_phase = np.empty(len(flt_files), dtype=np.float)
for i, im in enumerate(flt_files):
with fits.open(im) as hdu:
postarg1[i] = hdu[0].header['postarg1']
postarg2[i] = hdu[0].header['postarg2']
x_phase[i] = abs(np.modf(postarg1[i] / plate_scale)[0])
y_phase[i] = abs(np.modf(postarg2[i] / plate_scale)[0])
fig, ax = plt.subplots(1, 2, figsize=(12, 6))
ax[0].plot(postarg1, postarg2, 'k+-', ms=15, lw=1)
ax[0].set_xlabel('arcsec', fontsize='large')
ax[0].set_ylabel('arcsec', fontsize='large')
ax[1].plot(x_phase, y_phase, 'k+-', ms=15, lw=1)
ax[1].set_xlim(-0.05, 1)
ax[1].set_ylim(-0.05, 1)
ax[1].set_xlabel('X pixel phase', fontsize='large')
ax[1].set_ylabel('Y pixel phase', fontsize='large')
Before drizzling can be performed, a plate scale for the output image should be chosen. This will usually be dictated by the science needs. In theory, critical sampling of a PSF occurs with 2.355 pixels. The FWHM of the WFC3/IR detector at 1600 nm is ~0.151" (Table 7.5 of WFC Instrument Handbook). Using these values, a plate scale of 0.065"/pix will be selected for the drizzled output frame, approximately half of the native plate scale (0.1283"/pix).
The mechanism by which the drizzle algorithm improves sampling in the output image is by shrinking the input pixels before adding them to the output pixel grid (see Figure 2 in Fruchter and Hook (2002)). In practice this process is controlled in the AstroDrizzle package by adjusting the final_pixfrac parameter, which is the fractional size of the pixel to be used.
Below AstroDrizzle is called twice. The first time, no changes are made to the plate scale or the size of the drop. In the second call, the plate scale and pixfrac are reduced. For a four point dither stack, a final_pixfrac value of 0.8 is usually optimal.
drizzlepac.astrodrizzle.AstroDrizzle('*flt.fits',
output='f160w_noopt',
runfile='',
context=False,
build=True,
preserve=False,
clean=True,
driz_sep_bits='64,512',
final_bits='64,512',
final_wcs=True,
final_rot=0.)
drizzlepac.astrodrizzle.AstroDrizzle('*flt.fits',
output='f160w_opt',
runfile='',
context=False,
build=True,
preserve=False,
clean=True,
driz_sep_bits='64,512',
final_bits='64,512',
final_pixfrac=0.8,
final_wcs=True,
final_rot=0.,
final_scale=0.065)
The drizzled science and weight images produced from the first call to AstroDrizzle with no optimization of the plate scale and pixfrac are plotted below.
with fits.open('f160w_noopt_drz.fits') as hdu:
im1wcs = wcs.WCS(hdu[1].header)
sci1 = hdu[1].data
wht1 = hdu[2].data
norm1 = ImageNormalize(sci1, vmin=-0.2, vmax=25, stretch=LogStretch())
fig, ax = plt.subplots(1, 2, figsize=(16, 8), subplot_kw={'projection':im1wcs})
ax[0].imshow(sci1, norm=norm1, cmap='gray', origin='lower')
ax[1].imshow(wht1, cmap='gray', origin='lower')
The drizzled science image is on the left and the associated weight image is on the right, both without optimization of the plate scale and pixfrac.
To compare, the figure plotted below shows close ups of the same part of the sky from the two drizzled products.
radeclims = wcs.utils.pixel_to_skycoord([825, 875], [930, 980], im1wcs)
with fits.open('f160w_opt_drz.fits') as hdu:
im2wcs = wcs.WCS(hdu[1].header)
sci2 = hdu[1].data
wht2 = hdu[2].data
norm2 = ImageNormalize(sci2, vmin=-0.05, vmax=25, stretch=LogStretch())
fig = plt.figure(figsize=(16, 8))
ax = fig.add_subplot(1, 2, 1, projection=im1wcs)
ax.imshow(sci1, norm=norm2, cmap='gray')
ax.set_xlim(radeclims.to_pixel(im1wcs)[0])
ax.set_ylim(radeclims.to_pixel(im1wcs)[1])
ax = fig.add_subplot(1, 2, 2, projection=im2wcs)
ax.imshow(sci2, norm=norm2, cmap='gray')
ax.set_xlim(radeclims.to_pixel(im2wcs)[0])
ax.set_ylim(radeclims.to_pixel(im2wcs)[1])
The image on the left is without improved plate scale and pixfrac, and shows that the detector undersamples the PSF. The image on the right is the image with improved parameters where the resolution is greatly improved.
final_pixfrac parameter¶While the optimized final_pixfrac in the example above was chosen from experience with using a four-point dither, the value that should be used is not known a priori. The value could be different depending on several factors. For example, if the number of images is greater than four, the value used for final_pixfrac could be smaller since more images are available to fill in holes in the output grid. On the other hand, it is possible that no dithering was used during the observations. In that case, final_pixfrac should be left at 1.0, since shrinking the size of the drop could be detrimental.
Below, a series of experiments will be run to determine the best final_pixfrac value for the selected output plate scale, by varying pixfrac in steps of 0.1 over a range of values from 0.1 to 1.0.
pixfracs = np.arange(0.1, 1.1, 0.1)
for pixfrac in pixfracs:
outname = 'f160w_{:.1f}'.format(pixfrac)
drizzlepac.astrodrizzle.AstroDrizzle('*flt.fits',
output=outname,
runfile='',
context=False,
build=True,
preserve=False,
clean=True,
driz_sep_bits='64,512',
final_bits='64,512',
final_pixfrac=pixfrac,
final_wcs=True,
final_rot=0.,
final_scale=0.065)
When evaluating what value to use for final_pixfrac, THERE IS NO SINGLE METRIC THAT INDICATES WHAT VALUE TO USE, and several factors should be taken into account. The general philosophy is that the chosen value should improve the resolution of the image as much as possible, without causing any adverse effects.
The first thing to look for is an excessive number of holes in the science and weight images. The figure below shows the central region of the science and weight images produced by three different final_pixfrac values.
with fits.open('f160w_0.1_drz.fits') as hdu1:
sci1 = hdu1[1].data
wht1 = hdu1[2].data
with fits.open('f160w_0.8_drz.fits') as hdu2:
sci2 = hdu2[1].data
wht2 = hdu2[2].data
with fits.open('f160w_1.0_drz.fits') as hdu3:
sci3 = hdu3[1].data
wht3 = hdu3[2].data
fig,ax = plt.subplots(2, 3, figsize=(16, 10), sharex=True, sharey=True)
fig.subplots_adjust(hspace=0, wspace=0)
norm3 = ImageNormalize(wht2, vmin=400, vmax=600, stretch=LinearStretch())
ax[0, 0].imshow(sci1, norm=norm2, cmap='gray')
ax[0, 0].text(1640, 1920, 'final_pixfrac=0.1', fontsize='20', color='gold')
ax[0, 1].imshow(sci2, norm=norm2, cmap='gray')
ax[0, 1].text(1640, 1920, 'final_pixfrac=0.8', fontsize='20', color='gold')
ax[0, 2].imshow(sci3, norm=norm2, cmap='gray')
ax[0, 2].text(1640, 1920, 'final_pixfrac=1.0', fontsize='20', color='gold')
ax[1, 0].imshow(wht1, norm=norm3, cmap='gray')
ax[1, 1].imshow(wht2, norm=norm3, cmap='gray')
ax[1, 2].imshow(wht3, norm=norm3, cmap='gray')
ax[0, 0].set_xlim(radeclims.to_pixel(im2wcs)[0])
ax[0, 0].set_ylim(radeclims.to_pixel(im2wcs)[1])
This figure above shows the central region of the science and weight images produced by final_pixfrac values (left to right) of 0.1, 0.8 and 1.0. The top row is the science frame, the bottom is the weight image. The science image with the smallest final_pixfrac value of 0.1 shows a noisy background and holes in the image where no input pixels fall into the output grid because pixfrac is too small. Inspection of the weight map corresponding to the smallest final_pixfrac value shows many places with weights of zero, indicating that a final_pixfrac value of 0.1 is clearly too small and was only included in this example for illustrative purposes.
Another piece of information that can be useful is the amount of noise in the weight image. As suggested in the HST Dither Handbook Chapter 5 and in the DrizzlePac Handbook Section 6.3.3, statistics performed on the drizzled weight image should yield a RMS/median value less than ~0.2. This threshold controls the trade-off between improving image resolution versus increasing background noise due to pixel resampling.
The figure below shows the RMS/median as a function of final_pixfrac. One should take care to use the same region in the weight image as the region where the object of interest is located in the science image. If one is using the entire image for scientific analysis, then one should measure the statistics of the weight image where there is more variance.
whtlist = glob.glob('f160w_0.[0-9]*drz.fits')
std_med = np.empty(len(whtlist), dtype=float)
fraclist = np.empty(len(whtlist), dtype=float)
xlims = radeclims.to_pixel(im1wcs)[0].astype(int)
ylims = radeclims.to_pixel(im1wcs)[1].astype(int)
#Loop that measures statistics, also some information gathering
for i, im in enumerate(whtlist):
with fits.open(im) as hdu:
hdr = hdu[0].header
wht = hdu[2].data
if i == 0:
target = hdr['TARGNAME']
scale = str(hdr['D001SCAL'])
nimg = hdr['NDRIZIM']
wht_std = np.std(wht[ylims[0] : ylims[1], xlims[0] : xlims[1]])
wht_med = np.median(wht[ylims[0] : ylims[1], xlims[0] : xlims[1]])
std_med[i] = wht_std / wht_med
fraclist[i] = hdr['D001PIXF']
#Plotting commands
plt.clf()
plt.xlim(1.025, -0.025)
plt.ylim(0., 0.5)
plt.scatter(fraclist, std_med, s=50)
plt.axhline(0.2, ls='--', lw=3, c='r')
plt.xlabel('pixfrac', fontsize=18)
plt.ylabel('rms/median', fontsize=18)
plt.text(0.1, 0.45, target, fontsize=16, horizontalalignment='right')
plt.text(0.1, 0.40, 'Scale=' + str(scale) + '"', fontsize=16, horizontalalignment='right')
plt.text(0.1, 0.35, str(nimg) + ' images', fontsize=16, horizontalalignment='right')
The RMS/median increases steadily up to pixfrac of 0.6, where there is a change in the slope of the function, indicating that something odd might be happening beyond that point. Even though the RMS remains below 20% of the median all the way to final_pixfrac of 0.1, the higher slope in the curve combined with the information from inspection of the science images themselvse (the evidence of holes) indicates that the pixfrac should remain >0.6.
The final_pixfrac value has to be small enough to avoid degrading the combined output image, but large enough that when all images are "dropped" onto the final frame, coverage of the output frame is fairly uniform. In general, final_pixfrac should be slightly larger than the final output scale to allow some "spillover" to adjacent pixels. This will help avoid "holes" in the final product when a given pixel has been flagged as "bad" in several frames. As a rule of thumb, statistics performed on the drizzled weight image in the region of interest should yield an RMS value (standard deviation) that is less than 20% of the median value. This threshold is a balance between the benefits of improving the image resolution at the expense of increasing noise in the background.
Author: R. Avila, STScI ACS Team
Updated: December 12, 2018